text:
export CircleGraph, RandomCircleGraph, RandomCircleRepresentation
export isCircleGraph

"""
`list2locs(list)` takes a list of values as input. Every item on the list
should appear exactly twice. This function returns a dictionary mapping
items on the list to a pair of positive integers giving the two locations
of that element in the list. It's a helper function used by
`CircleGraph`.
"""
function list2locs(list::Array{T,1}) where {T}
    # build a map giving the locations of each symbol on the list
    err_msg = "Invalid input list"
    # check that the list length is even
    n = length(list)
    if n % 2 == 1
        error(err_msg)
    end

    positions = Dict{T,Vector{Int}}()

    for i = 1:n
        elt::T = list[i]
        if haskey(positions, elt)  # already have 1st position
            if positions[elt][2] > 0  # and already have 2nd; that's bad
                error(err_msg)
            else
                positions[elt][2] = i
            end
        else # never seen this element before
            positions[elt] = [i, 0]
        end
    end
    return positions
end

"""
`CircleGraph(list)` takes a list of elements and creates a circle
graph from them. Every element on the list should appear exactly
twice on the list, or else an error is thrown.
"""
function CircleGraph(list::Array{T,1}) where {T}
    positions = list2locs(list)

    elts = collect(keys(positions))
    n = length(elts)

    G = SimpleGraph{T}()
    for v in elts
        add!(G, v)
    end

    for i = 1:n-1
        v = elts[i]
        a = positions[v]
        for j = i+1:n
            w = elts[j]
            b = positions[w]
            if a[1] < b[1] < a[2] < b[2] || b[1] < a[1] < b[2] < a[2]
                add!(G, v, w)
            end
        end
    end

    cache_save(G, :CircleRepresentation, list)
    cache_save(G, :name, "Circle graph")
    return G
end

"""
`CircleGraph(str::String)` creates a circle graph from
a string of ASCII characters; each character must appear exactly
twice in the string.
"""
function CircleGraph(str::String)
    list = [ch for ch in str]
    return CircleGraph(list)
end

"""
`RandomCircleRepresentation(n)` returns a `2n` long list
consisting of the numbers `1:n` each exacty twice in some
random order.
"""
function RandomCircleRepresentation(n::Int)
    list = [collect(1:n); collect(1:n)]
    p = randperm(2n)
    return list[p]
end

"""
`RandomCircleGraph(n)` creates a random circle graph with `n`
vertices.
"""
function RandomCircleGraph(n::Int)
    return CircleGraph(RandomCircleRepresentation(n))
end

# This code by Tara

# returns true if G is a circle graph, false otherwise

using SimpleGF2

"""
`isCircleGraph(G)` returns `true` if `G` is a circle graph
and `false` otherwise.
"""
function isCircleGraph(G::SimpleGraphs.SimpleGraph)
    V = eltype(G)
    vertNum = Dict{V,Int}()
    counter = 1
    verts = length(vlist(G))
    z = zeros(GF2, verts * verts + 1)
    C = deepcopy(z)
    vertices = vlist(G)
    for v in vertices
        vertNum[v] = counter
        counter = counter + 1
    end
    for e in elist(G)
        x = deepcopy(z)
        x[(vertNum[e[1]]-1)*verts+vertNum[e[2]]] = 1
        x[(vertNum[e[2]]-1)*verts+vertNum[e[1]]] = 1
        x[verts*verts+1] = 1
        C = hcat(C, x)
    end
    n = length(vertices)
    for v in vertices
        for i = 1:n
            for j = 1+i:n
                v1 = vertices[i]
                v2 = vertices[j]
                if (has(G, v, v1) && !has(G, v, v2)) ||
                   (has(G, v, v2) && !has(G, v, v1)) ||
                   v == v1 ||
                   v == v2 ||
                   v1 == v2
                    continue
                elseif has(G, v, v1) && has(G, v, v2)
                    if !has(G, v1, v2)
                        x = deepcopy(z)
                        x[(vertNum[v]-1)*verts+vertNum[v1]] = 1
                        x[(vertNum[v]-1)*verts+vertNum[v2]] = 1
                        x[(vertNum[v1]-1)*verts+vertNum[v2]] = 1
                        x[(vertNum[v2]-1)*verts+vertNum[v1]] = 1
                        x[verts*verts+1] = 1
                        C = hcat(C, x)
                    end
                elseif !has(G, v, v1) && !has(G, v, v2)
                    if has(G, v1, v2)
                        x = deepcopy(z)
                        x[(vertNum[v]-1)*verts+vertNum[v1]] = 1
                        x[(vertNum[v]-1)*verts+vertNum[v2]] = 1
                        C = hcat(C, x)
                    end
                end
            end
        end
    end
    C = C'
    try
        solve_augmented(C)
    catch
        return false
    end
    return true
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Random, Distributions, StatsBase

@testset "bat_sample" begin
    @testset "IIDSampling" begin
        dist = MvNormal([0.4, 0.6], [2.0 1.2; 1.2 3.0])

        @test length(@inferred(bat_sample(dist, IIDSampling(nsamples = 10^3))).result) == 10^3

        @test @inferred(bat_sample(Random.GLOBAL_RNG, dist)).result isa DensitySampleVector
        @test @inferred(bat_sample(dist, BAT.IIDSampling())).result isa DensitySampleVector
        @test @inferred(bat_sample(Random.GLOBAL_RNG, dist, BAT.IIDSampling())).result isa DensitySampleVector

        samples = @inferred(bat_sample(dist, IIDSampling(nsamples = 10^5))).result
        @test isapprox(mean(samples.v), [0.4, 0.6]; rtol = 0.05)
        @test isapprox(cov(samples.v), [2.0 1.2; 1.2 3.0]; rtol = 0.05)
        @test all(isequal(1), samples.weight)

        resamples = @inferred(bat_sample(samples, OrderedResampling(nsamples = length(samples)))).result
        @test samples == resamples

        dist_bmode = @inferred(bat_findmode(dist)).result
        @test @inferred(length(dist_bmode)) == 2

        dist_sample_vector_bmode = @inferred(bat_findmode(samples)).result
        @test @inferred(length(dist_sample_vector_bmode)) == 2

        isapprox(var(bat_sample(Normal(), BAT.IIDSampling(nsamples = 10^3)).result), [1], rtol = 10^-1)
    end
end

ext: jl


text:
import FastGroupBy: grouptwo!
using BenchmarkTools

const N = 200_000_000
const K = 100

srand(1);
x = rand(1:N÷K, N);
y = rand(N);

@time grouptwo!(x,y); # 2.95


text:
"""
    StartingPtPrior{T}

Types inheriting from the abstract type `StartingPtPrior` indicate the prior
that is put on the starting point of the observed path of some stochastic
process. `T` denotes the DataType of the starting point.
"""
abstract type StartingPtPrior{T} end

@doc raw"""
    Base.rand(G::StartingPtPrior, [z, ρ=0.0])

Sample a new starting point according to its prior distribution. An
implementation with arguments `z`, `ρ` implements a preconditioned
Crank-Nicolson scheme with memory parameter `ρ` and a current **non-centered
variable** `z`. `z` is also referred to as the **driving noise**.
"""
Base.rand(G::StartingPtPrior, z, ρ=0.0) = error("Not Implemented")
Base.rand(G::StartingPtPrior) = error("Not Implemented")

"""
    start_pt(z, G::StartingPtPrior, P)

Compute a new starting point from the white noise for a given posterior
distribution obtained from combining prior `G` and the likelihood encoded by the
object `P`.
"""
start_pt(z, G::StartingPtPrior, P) = error("Not Implemented")

"""
    start_pt(z, G::StartingPtPrior)

Compute a new starting point from the white noise for a given prior
distribution `G`
"""
start_pt(z, G::StartingPtPrior) = error("Not Implemented")

"""
    inv_start_pt(y, G::StartingPtPrior, P)

Compute the driving noise that is needed to obtain starting point `y` under
prior `G` and the likelihood in `P`
"""
inv_start_pt(y, G::StartingPtPrior, P) = error("Not Implemented")

"""
    logpdf(G::StartingPtPrior, y)

log-probability density function evaluated at `y` of a prior distribution `G`
"""
Distributions.logpdf(G::StartingPtPrior, y) = error("Not Implemented")

#===============================================================================
                            Struct definitions
===============================================================================#

"""
    struct KnownStartingPt{T} <: StartingPtPrior{T}
        y::T
    end

Indicates that the starting point is known and stores its value in `y`

    KnownStartingPt(y::T) where T

Base constructor.
"""
struct KnownStartingPt{T} <: StartingPtPrior{T}
    y::T
    KnownStartingPt(y::T) where T = new{T}(y)
end

"""
    struct GsnStartingPt{T,S,TM} <: StartingPtPrior{T} where {S}
        μ::T
        Σ::S
        Λ::S
        μ₀::T
        Σ₀::UniformScaling
    end

Indicates that the starting point is equipped with a Gaussian prior with
mean `μ` and covariance matrix `Σ` (and pre-computed precision `Λ`:=`Σ`⁻¹).
Sampling is always done via non-centred parametrization, by sampling white noise
`z` according to Gaussian with zero mean and identity covariance: `μ₀` and `Σ₀`,
and then transforming `z` to a variable with mean and covariance `μ` and `Σ`.

    GsnStartingPt(μ::T, Σ::S)

Base constructor. It initialises the mean `μ` and covariance `Σ` parameters and
`Λ` is set according to `Λ`:=`Σ`⁻¹.
"""
struct GsnStartingPt{T,S,TM} <: StartingPtPrior{T} where {S}
    μ::T
    Σ::S
    Λ::S
    μ₀::T
    Σ₀::UniformScaling

    GsnStartingPt{T,S,TM}(
        μ::T,
        Σ::S,
        Λ::S,
        μ₀::T,
        Σ₀::UniformScaling
    ) where {T,S,TM} = new{T,S,TM}(μ, Σ, Λ, μ₀, Σ₀)
end

GsnStartingPt(μ, Σ) = GsnStartingPt(μ, Σ, ismutable(μ))

function GsnStartingPt(μ::T, Σ, ::Val{true}) where T
    typeof(Σ) <: Vector && length(Σ) == 1 && (Σ = reshape(Σ, 1, 1))
    ismutable(Σ)::Val{true}
    @assert size(Σ,1) == size(Σ,2) == length(μ)
    GsnStartingPt{T,typeof(Σ),:mutable}(
        μ, Σ, inv(Σ), zeros(eltype(T), length(μ)), I
    )
end

function GsnStartingPt(μ::T, Σ, ::Val{false}) where T
    typeof(Σ) <: SVector && length(Σ) == 1 && (Σ = SMatrix{1,1}(Σ))
    ismutable(Σ)::Val{false}
    @assert size(Σ,1) == size(Σ,2) == length(μ)
    GsnStartingPt{T,typeof(Σ),:immutable}(μ, Σ, inv(Σ), zero(T), I)
end

#===============================================================================
                Sampling a starting point according to a prior
===============================================================================#

"""
    rand([rng::Random.AbstractRNG], G::GsnStartingPt, z, ρ=0.0)

Sample new white noise using Crank-Nicolson scheme with memory parameter `ρ` and
a previous value of the white noise stored inside object `G`
"""
function Base.rand(rng::Random.AbstractRNG, G::GsnStartingPt, z, ρ=0.0)
    zᵒ = rand(rng, Gaussian(G.μ₀, G.Σ₀))
    sqrt(1-ρ)*zᵒ + sqrt(ρ)*z # preconditioned Crank-Nicolson
end

Base.rand(G::GsnStartingPt, z, ρ=0.0) = rand(Random.GLOBAL_RNG, G, z, ρ)

"""
    rand([rng::Random.AbstractRNG], G::GsnStartingPt)

Sample new starting point according to its prior distribution.
"""
function Base.rand(rng::Random.AbstractRNG, G::GsnStartingPt)
    rand(rng, Gaussian(G.μ, G.Σ))
end

Base.rand(G::GsnStartingPt) = rand(Random.GLOBAL_RNG, Gaussian(G.μ, G.Σ))

"""
    rand([rng::Random.AbstractRNG], G::KnownStartingPt, args...)

Starting point is known. Nothing can be sampled. Returning known starting point.
"""
function Base.rand(rng::Random.AbstractRNG, G::KnownStartingPt, args...)
    rand(Random.GLOBAL_RNG, G, args...)
end

Base.rand(G::KnownStartingPt, args...) = G.y

#===============================================================================
                Compute the starting point from random seed
===============================================================================#

"""
    start_pt(z, G::GsnStartingPt, P)

Compute a new starting point from the white noise for a given posterior
distribution obtained from combining prior `G` and the likelihood encoded by the
object `P`.
"""
function start_pt(z, G::GsnStartingPt, P)
    μ_post = (P.H[1] + G.Λ) \ (P.Hν[1] + G.Λ * G.μ)
    Σ_post = inv(P.H[1] + G.Λ)
    Σ_post = 0.5 * (Σ_post + Σ_post') # remove numerical inaccuracies
    unwhiten(Σ_post, z) + μ_post
end

"""
    start_pt(z, G::GsnStartingPt)

Compute a new starting point from the white noise for a given prior
distribution `G`
"""
start_pt(z, G::GsnStartingPt) = unwhiten(G.Σ, z) + G.μ

"""
    start_pt(z, G::KnownStartingPt, P)

Return starting point
"""
start_pt(z, G::KnownStartingPt, P) = G.y

#TODO check if this is needed
"""
    start_pt(G::KnownStartingPt, P)

Return starting point
"""
start_pt(G::KnownStartingPt, P) = G.y

#===============================================================================
                Compute the random seed from the starting point
===============================================================================#

"""
    inv_start_pt(y, G::GsnStartingPt, P)

Compute the driving noise that is needed to obtain starting point `y` under
prior `G` and the likelihood in `P`
"""
function inv_start_pt(y, G::GsnStartingPt, P)
    μ_post = (P.H[1] + G.Λ) \ (P.Hν[1] + G.Λ * G.μ)
    Σ_post = inv(P.H[1] + G.Λ)
    Σ_post = 0.5 * (Σ_post + Σ_post')
    try
        whiten(Σ_post, y-μ_post)
    catch e
        print("incorrect matrix: ", Σ_post, "\n\n")
    end
    whiten(Σ_post, y-μ_post)
end

"""
    inv_start_pt(y, G::KnownStartingPt, P)

Starting point known, no need for dealing with white noise, use convention of
returning `y`
"""
inv_start_pt(y, G::KnownStartingPt, P) = y

#===============================================================================
                            compute logpdf
===============================================================================#

"""
    logpdf(::KnownStartingPt, y)

Nothing to do so long as `y` is equal to the known starting point inside `G`
"""
function Distributions.logpdf(G::KnownStartingPt, y)
    (G.y == y) ? 0.0 : error("Starting point is known, but a different value ",
                             "was passed to logpdf.")
end

"""
    logpdf(G::GsnStartingPt, y)

log-probability density function evaluated at `y` of a prior distribution `G`
"""
Distributions.logpdf(G::GsnStartingPt, y) = logpdf(Gaussian(G.μ, G.Σ), y)


text:
using GeometricEquations
using GeometricEquations: arrtype, datatype
using GeometricEquations: check_initial_conditions, check_parameters
using GeometricEquations: _functions, _solutions
using Test

@test AbstractEquationODE  <: GeometricEquation
@test AbstractEquationDAE  <: GeometricEquation
@test AbstractEquationSDE  <: GeometricEquation
@test AbstractEquationPODE <: GeometricEquation
@test AbstractEquationPDAE <: GeometricEquation
@test AbstractEquationPSDE <: GeometricEquation

struct TestEquation <: GeometricEquation{Nothing,Nothing,Nothing} end

testeq = TestEquation()

@test functions(testeq) == NamedTuple()
@test solutions(testeq) == NamedTuple()

@test_throws ErrorException functions(testeq, NullParameters())
@test_throws ErrorException solutions(testeq, NullParameters())

@test_throws ErrorException _functions(testeq)
@test_throws ErrorException _solutions(testeq)

@test_throws ErrorException _functions(testeq, NullParameters())
@test_throws ErrorException _solutions(testeq, NullParameters())

@test_throws ErrorException invariants(testeq)
@test_throws ErrorException parameters(testeq)
@test_throws ErrorException periodicity(testeq)

@test_throws ErrorException check_initial_conditions(testeq, NamedTuple())

@test_throws ErrorException check_parameters(testeq, NullParameters())
@test_throws ErrorException check_parameters(testeq, NamedTuple())
@test_throws ErrorException check_parameters(testeq, 1)

@test_throws ErrorException datatype(testeq, NamedTuple())
@test_throws ErrorException arrtype(testeq, NamedTuple())

@test hassolution(testeq) == false
@test hasvectorfield(testeq) == false
@test hasprimary(testeq) == false
@test hassecondary(testeq) == false

@test hasinvariants(testeq) == false
@test hasparameters(testeq) == false
@test hasperiodicity(testeq) == false

@test hashamiltonian(testeq) == false
@test haslagrangian(testeq) == false


text:
P = [
    1 2 3 1e3 1e3 1e3 1e3
    3 1e3 2 1 1e3 1e3 1e3
    1e3 3 1e3 1e3 1 2 1e3
    1 1e3 1e3 1e3 3 1e3 2
    1 3 1e3 1e3 2 1e3 1e3
    2 1e3 3 1e3 1 1e3 1e3
    1e3 1e3 2 1 1e3 1e3 3
    1 3 1e3 1e3 2 1e3 1e3
    1 3 1e3 1e3 2 1e3 1e3
    1e3 1e3 2 3 1e3 1e3 1
    1 2 3 1e3 1e3 1e3 1e3
    1 2 3 1e3 1e3 1e3 1e3
    1e3 1e3 1 2 1e3 1e3 3
    1e3 1e3 3 2 1e3 1e3 1e3
    2 1 1e3 3 1e3 1e3 1e3
    1 2 1e3 3 1e3 1e3 1e3
    1e3 1 1e3 2 1e3 1e3 3
    1e3 3 1 1e3 1e3 1e3 2
    1e3 1e3 1e3 3 1 2 1e3
    1e3 2 1e3 3 1e3 1e3 1
    3 2 1e3 1 1e3 1e3 1e3
    3 2 1 1e3 1e3 1e3 1e3
    1 2 3 1e3 1e3 1e3 1e3
    1e3 1e3 2 3 1e3 1e3 1
    1e3 1 3 2 1e3 1e3 1e3
    1 1e3 2 1e3 3 1e3 1e3
    1 2 1e3 1e3 1e3 1e3 3
    1 2 1e3 3 1e3 1e3 1e3
    1 2 1e3 3 1e3 1e3 1e3
    2 1 3 1e3 1e3 1e3 1e3
    1e3 1e3 3 1e3 1e3 1 2
    3 2 1e3 1e3 1 1e3 1e3
    1 1e3 2 1e3 3 1e3 1e3
    1e3 2 1e3 1e3 1 1e3 3
    3 1e3 2 1 1e3 1e3 1e3
    1e3 1 3 1e3 1e3 1e3 2
    1e3 1e3 2 1 1e3 3 1e3
    1e3 3 1 2 1e3 1e3 1e3
    1 1e3 2 3 1e3 1e3 1e3
    1 1e3 3 2 1e3 1e3 1e3
    3 1 1e3 1e3 1e3 1e3 2
    1 1e3 2 1e3 3 1e3 1e3
    1 1e3 1e3 1e3 2 3 1e3
    3 1e3 2 1e3 1e3 1e3 1
    2 1 3 1e3 1e3 1e3 1e3
    1e3 3 2 1e3 1e3 1 1e3
    1e3 2 1e3 3 1 1e3 1e3
    1 3 1e3 1e3 1e3 1e3 2
    2 3 1e3 1 1e3 1e3 1e3
    1e3 1 1e3 2 1e3 1e3 3
]


text:
"""
    set_Branch_Scheme!(x::BnBSolver,BM::String)

Sets the search scheme to "best", "breadth", or "depth" first schemes.
"""
function set_Branch_Scheme!(x::BnBSolver,BM::String)
  if (BM == "best")
    x.Node_Select = NS_best
    x.Branch_Sto = BM_depth_best!
  elseif (BM == "breadth")
    x.Node_Select = NS_depth_breadth
    x.Branch_Sto = BM_breadth!
  elseif (BM == "depth")
    x.Node_Select = NS_depth_breadth
    x.Branch_Sto = BM_depth_best!
  else
    error("Invalid branching method")
  end
end

"""
    set_Bisect_Func!(x::BnBSolver,BF::String,nx::Int64)

Sets the bisection function to BF = "relative midpoint" or
BF = "absolute midpoint" and disregards the first nx components of the interval
box storage.
"""
function set_Bisect_Func!(x::BnBSolver,BF::String,nx::Int64)
  if (nx<0)
    if (BF == "relative midpoint")
      x.Bisect_Func = Bisect_Rel
    elseif (BF == "absolute midpoint")
      x.Bisect_Func = Bisect_Abs
    else
      error("Invalid bisection method")
    end
  else
    if (BF == "relative midpoint")
      x.Bisect_Func = (S::BnBSolver,B::BnBModel,N) -> Bisect_Rel_Imp(S,B,N,nx)
    elseif (BF == "absolute midpoint")
      x.Bisect_Func = (S::BnBSolver,B::BnBModel,N) -> Bisect_Abs_Imp(S,B,N,nx)
    else
      error("Invalid bisection method")
    end
  end
end

"""
    set_Verbosity!(x::BnBSolver,VB::String)

Sets the verbosity (console output) to either "None", "Normal", or "Full".
"""
function set_Verbosity!(x::BnBSolver,VB::String)
  if (VB == "None"||VB == "Normal"||VB == "Full")
    x.Verbosity = VB
  else
    error("Invalid Verbosity Option")
  end
end

"""
    set_to_default!(x::BnBSolver)

Returns the B&B solver to the default settings (does not include problems or processing routines.
"""
function set_to_default!(x::BnBSolver)
  x.Term_Check = Term_Check
  x.Branch_Sto = BM_depth_best!
  x.Node_Select = NS_best
  x.Bisect_Func = Bisect_Abs
  x.Verbosity = "Normal"
  x.max_iter = Inf
  x.iter_lim = false
  x.max_nodes = 1E6
  x.BnB_atol = 1E-4
  x.BnB_rtol = 1E-4
  x.itr_intv = 1
  x.hdr_intv = 20
  x.converged = Conv_Check
  x.BnB_digits = 3
  x.hist_return = false
  x.opt = []
  x.exhaust = false
end


text:
using WTExperiments
using Test

@testset "WTExperiments.jl" begin
    # Write your tests here.
end


text:
module Geometry

safe = 1000

function setSafe(s)
    safe = s
end

# returns a function that describes a parabola defined by a focus p and a directrix of vertical component ly
function parabola(p::Tuple{Real, Real}, ly::Real)
   if p[2] == ly
      return nothing
   else
      return function (x)
         return (x^2 - 2p[1]*x + p[1]^2 + p[2]^2 -ly^2)/2(p[2] - ly)
      end
   end
end

function parabolaIntersection(p::Tuple{Real, Real}, q::Tuple{Real, Real}, ly::Real)
   f = parabola(p, ly)
   g = parabola(q, ly)

   if p[2] == ly && q[2] == ly # both degenerate parabolas
       midpoint = (p[1] + q[1])/2
       return (midpoint, ly + safe)
   elseif p[2] == ly
      return (p[1], g(p[1]))
   elseif q[2] == ly
      return (q[1], f(q[1]))
   end

   r = p[2] - ly
   s = q[2] - ly

   a = 1/r - 1/s
   b = 2(q[1]/s - p[1]/r)
   c = (p[1]^2 + p[2]^2 - ly^2)/r - (q[1]^2 + q[2]^2 - ly^2)/s

   delta = b^2 - 4a*c # discriminant

   if a == 0 # generating points of both parabolas have the same y
      x = (p[1] + q[1])/2
      return (x, f(x))
   end

   if delta < 0 # no intersection, should only happen when the first couple of points have the same y coordinate
      return nothing
   else
      x1 = (-b - sqrt(delta))/2a
      x2 = (-b + sqrt(delta))/2a

      if p[2] > q[2] # use left intersection
         if x1 < x2
            return (x1, f(x1))
         else
            return (x2, f(x2))
         end
      else # use right intersection
         if x1 > x2
            return (x1, f(x1))
         else
            return (x2, f(x2))
         end
      end
   end
end

function circumcircle(A::Tuple{Real, Real}, B::Tuple{Real, Real}, C::Tuple{Real, Real})
   # P is the midpoint between A and B
   P = ((A[1] + B[1])/2, (A[2] + B[2])/2)

   # v is a vector orthogonal to B-A; P + λv is the perpendicular bisector of AB, therefore the center of the circumcircle is in it
   v = (A[2] - B[2], B[1] - A[1])

   # Q is the midpoint between B and C
   Q = ((B[1] + C[1])/2, (B[2] + C[2])/2)

   # w is a vector orthogonal to C-B; Q + μw is the perpendicular bisector of BC, therefore the center of the circumcircle is in it
   w = (B[2] - C[2], C[1] - B[1])

   λ = (w[2]*(P[1] - Q[1]) + w[1]*(Q[2] - P[2])) / (w[1]*v[2] - w[2]*v[1])

   O = (P[1] + λ*v[1], P[2] + λ*v[2])

   r = sqrt((O[1] - A[1])^2 + (O[2] - A[2])^2)

   return O, r
end

function distance(A::Tuple{Real, Real}, B::Tuple{Real, Real})
   return sqrt(distanceSquared(A, B))
end

function distanceSquared(A::Tuple{Real, Real}, B::Tuple{Real, Real})
   return (A[1]-B[1])^2 + (A[2]-B[2])^2
end

# tests if middle arc of focus b will converge between foci a and c when the sweep line reaches ly
function arcWillConverge(a::Tuple{Real, Real}, b::Tuple{Real, Real}, c::Tuple{Real, Real})
   det = (b[1] - a[1])*(c[2] - a[2]) - (c[1] - a[1])*(b[2] - a[2])

   return (det < 0)
end

function addVector(a::Tuple{Real, Real}, b::Tuple{Real, Real})
   return (a[1] + b[1], a[2] + b[2])
end

function subVector(a::Tuple{Real, Real}, b::Tuple{Real, Real})
   return (a[1] - b[1], a[2] - b[2])
end

function rotateVectorCCW(a::Tuple{Real, Real})
   return (-a[2], a[1])
end

function multVector(c::Real, a::Tuple{Real, Real})
   return (c*a[1], c*a[2])
end

#function someParabolaIsWrong() # for debugging
#   NLINES = 1000
#   NPOINTSFOREACHLINE = 1000
#
#   function parabolaIsWrong(p::Tuple{Real, Real}, ly::Real)
#      MAX_ERROR = 10^(-5.0)
#
#      X = linspace(0, 100, 1000)
#      f = parabola(p, ly)
#
#      distanceFromLine = map(x -> f(x) - ly, X)
#      distanceFromPoint = map(x -> distance(p, (x, f(x))), X)
#
#      return any(i -> i > MAX_ERROR, distanceFromLine - distanceFromPoint)
#   end
#
#   lines = rand(0:20.0, NLINES)
#
#   for ly in lines
#      PX = rand(-100.0:100.0, NPOINTSFOREACHLINE)
#      PY = rand(ly:ly+200.0, NPOINTSFOREACHLINE)
#
#      for i in 1:NPOINTSFOREACHLINE
#         if parabolaIsWrong((PX[i], PY[i]), ly)
#            return true
#         end
#      end
#   end
#
#   return false
#end

end # module

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
using QuantumESPRESSOParser
using Test

@testset "QuantumESPRESSOParser.jl" begin
    include("Inputs/Inputs.jl")
    # include("Outputs/PWscf.jl")
end


text:
function relu(x)
  return max(x, 0.)
end

function sigmoid(x)
  return 1./(1 + exp(-x))
end

function softmax(x)
  return exp(x)./sum(exp(x),2)
end


text:
println("Starting runtests.jl $(join(ARGS, " ")) ...")
using FunctionalDataUtils, FactCheck
FactCheck.setstyle(:compact)

shouldtest(f, a) = length(ARGS) == 0 || in(a, ARGS) ? facts(f, a) : nothing
shouldtestcontext(f, a) = length(ARGS) < 2 || a == ARGS[2] ? facts(f, a) : nothing

shouldtest("machinelearning") do
    shouldtestcontext("randbase") do
        # @test_equal size(randbase(2,10)) (2,10)
        # @test_almostequal sqrt(sum(randbase(2,10).^2, 2)) [1. 1.]'
        # @test_almostequal (@p randbase 2 10  | (.^) _ 2 |  sum _ 2) [1. 1.]'
    end
    shouldtestcontext("randproj") do
    end
    shouldtestcontext("loocv") do
    end

end

shouldtest("computing") do
    shouldtestcontext("fasthash") do
        @fact fasthash(1) --> "20eae029a26a15420f9ed6d7a9999e823f251b9d44adea5a504e1d8a4f7c5512"
        @fact fasthash(1.) --> "cb5415fc791f65b70b303678e3601d60ac527cc4cdd3a10d1f85b7423b270112"
        @fact fasthash('a') --> "94923c7a4899a5baf82f29551f5e4652e3b7ee50699bbd0b4ddba1538e5963b2"
        @fact fasthash("a") --> "ca978112ca1bbdcafac231b39a23dc4da786eff8147c4e72b9807785afee48bb"
        @fact fasthash(["a"]) --> "da3811154d59c4267077ddd8bb768fa9b06399c486e1fc00485116b57c9872f5"
        @fact fasthash([1]) --> "7c9fa136d4413fa6173637e883b6998d32e1d675f88cddff9dcbcf331820f4b8"
        @fact fasthash(UInt8[1]) --> "4bf5122f344554c53bde2ebb8cd2b7e3d1600ad631c385a5d7cce23c7785459a"
        @fact fasthash(Int8[1]) --> "4bf5122f344554c53bde2ebb8cd2b7e3d1600ad631c385a5d7cce23c7785459a"
        @fact fasthash(UInt16[1]) --> "47dc540c94ceb704a23875c11273e16bb0b8a87aed84de911f2133568115f254"
        @fact fasthash(Int16[1]) --> "47dc540c94ceb704a23875c11273e16bb0b8a87aed84de911f2133568115f254"
        @fact fasthash(UInt32[1]) --> "67abdd721024f0ff4e0b3f4c2fc13bc5bad42d0b7851d456d88d203d15aaa450"
        @fact fasthash(Int32[1]) --> "67abdd721024f0ff4e0b3f4c2fc13bc5bad42d0b7851d456d88d203d15aaa450"
        @fact fasthash(UInt64[1]) --> "7c9fa136d4413fa6173637e883b6998d32e1d675f88cddff9dcbcf331820f4b8"
        @fact fasthash(Int64[1]) --> "7c9fa136d4413fa6173637e883b6998d32e1d675f88cddff9dcbcf331820f4b8"
        @fact fasthash(Any[1,"test"]) --> "b59f20a7e3513ef702971ec5cbd288f19d760df46d49d5d27edb82fe73fc257b"
        data = Any[1,1., 'a', "a", rand(10,200), zeros(1000,10000), Dict("a" => 1, 2 => "b")]
        r = [fasthash(x) for x in data]
        @fact all(r .!= nothing) --> true   # just to avoid showing "0 facts verified". we did survive the above, after all.
        @fact fasthash(1:2) --> "48adce367d9daa9d76033a4693b20d7fab23d357bf97e49a036f2dd32b1afff4"
        @fact fasthash(1:3:10) --> "dc29cadf1701b9023b1170a8101b40928770dbbd967e2740afa4986499522ca3"
    end
end

shouldtest("computervision") do
    shouldtestcontext("warp") do
        img = rand(10,20,3)
        a = rand(3:9,2,4)
        b = copy(a)
        b[:,end] = 10
        @fact size(warp(img,a,b))  -->  (10,20,3)
        @fact size(warp(img,a,b,(8,9)))  -->  (8,9,3)

    end
    shouldtestcontext("iimg") do
        a = [1 2 3]
        @fact size(iimg(a)) --> (1,3)
        @fact typeof(iimg(a)) --> Array{Float64,2}
        @fact iimg(a) --> [1 3 6]
        @fact iimg([1 2 3; 4 5 6]) --> [1.0 3.0 6.0; 5.0 12.0 21.0]
        @fact iimg(ones(1,2,3)) --> float(cat(3,[1 2],[2 4],[3 6]))
    end

    shouldtestcontext("interp3") do
        a = [o for m in 1:2, n in 1:2, o in 1:2]
        @fact interp3(a,1.,1.,1.5) --> 1.5
        @fact interp3(a,2.,2.,1.5) --> 1.5
        a = [n for m in 1:2, n in 1:2, o in 1:2]
        @fact interp3(a,1.,1.5,1.) --> 1.5
        @fact interp3(a,2.,1.5,2.) --> 1.5
        a = [m for m in 1:2, n in 1:2, o in 1:2]
        @fact interp3(a,1.5,1.,1.) --> 1.5
        @fact interp3(a,1.5,2.,2.) --> 1.5
    end

    shouldtestcontext("interp3with01coords") do
        a = [o for m in 1:2, n in 1:2, o in 1:2]
        @fact interp3with01coords(a,0.,0.,0.5) --> 1.5
        @fact interp3with01coords(a,1.,1.,0.5) --> 1.5
        a = [n for m in 1:2, n in 1:2, o in 1:2]
        @fact interp3with01coords(a,0.,0.5,0.) --> 1.5
        @fact interp3with01coords(a,1.,0.5,1.) --> 1.5
        a = [m for m in 1:2, n in 1:2, o in 1:2]
        @fact interp3with01coords(a,0.5,0.,0.) --> 1.5
        @fact interp3with01coords(a,0.5,1.,1.) --> 1.5
    end

    shouldtestcontext("resize") do
        @fact resize(ones(Int,1,1), [1, 1]) --> ones(Int, 1, 1)
        @fact resize(ones(Float64,1,1), 1) --> ones(Float64, 1, 1)
        @fact resize(ones(Float64,1,1,1), [1, 1, 1]) --> ones(Float64, 1, 1, 1)
        @fact resize(ones(Int,1,1), [1, 10]) --> ones(Int, 1, 10)
        @fact maximum(abs(resize([1. 2.], [1, 4])-[1. 1.33333 1.66666 2.]))<0.01 --> true
        @fact resize([1. 3.], [1, 3]) --> [1. 2. 3.]
    end

    shouldtestcontext("bwlabel") do
        img = [1 2 3]
        @fact bwlabel(img) --> img

        img = [1 2 3 0 5 6]
        @fact bwlabel(img) --> [1 2 3 1 5 6]
        @fact bwlabel(img,4) --> [1 2 3 4 5 6]

        img = [1 1 0 5 2 2; 3 3 0 0 4 4]
        bwlabel(img)
    end

    shouldtestcontext("border") do
        img = zeros(5,7)
        img[3,4] = 1
        @fact border(img)  -->  img
        img[2:4,3:6] = 1
        b = copy(img)
        b[3,4:5] = 0
        @fact border(img)  -->  b
    end

    shouldtestcontext("bwdist") do
        a = [0 0 0 1 0 0]
        @fact bwdist(a)  -->  [3 2 1 0 1 2]
        @fact bwdist([1 10;2 12], [1 10; 1 10])  -->  [1,2]
    end

    shouldtestcontext("rle") do
        data = [1 1 2 2 2 3 3 3 3]
        @fact unrle(rle(data)) --> data
        data = rand(1:5, 1, 1000)
        @fact unrle(rle(data)) --> data
        for t = [Int16, Int32, Int64, Float32, Float64]
            for s1 = 1:10, s2 = 1:10
                data = rand(t, s1, s2)
                @test_equal unrle(rle(data)) data
            end
        end
    end

    shouldtestcontext("resizeminmax") do
        resize(a) = size(resizeminmax(a, [30,40], [130, 140]))
        @fact resize(ones(1,1)) --> (30,40)
        @fact resize(ones(100,200)) --> (85,140)
        @fact resize(ones(10,200)) --> (30,140)
        @fact resize(ones(50,10)) --> (125,40)
        @fact resize(ones(50,60)) --> (50,60)
        @fact resize(ones(1,600)) --> (30,140)
        @fact resize(ones(600,1)) --> (130,40)
    end

    shouldtestcontext("blocks") do
        a = ones(16,32)
        ind = stridedblocksub(a, 8; keepshape = true)
        @fact size(ind) --> (2,4)
        ind = stridedblocksub(a, 8)
        @fact len(ind) --> 8
        ind =  stridedblocksub(a, 8, 1)
        @fact len(ind) --> 9*25
        a = flatten(Any[ones(2,2) 2*ones(2,2); 3*ones(2,2) 4*ones(2,2)])
        subs = stridedblocksub(a, 2)
    end

    shouldtestcontext("meshgrid") do
        @fact meshgrid(2,3)  -->  [1 2 1 2 1 2; 1 1 2 2 3 3]
        @fact meshgrid(1:2,1:3)  -->  [1 2 1 2 1 2; 1 1 2 2 3 3]
        @fact meshgrid(2:3,2:4)  -->  [2 3 2 3 2 3; 2 2 3 3 4 4]
        @fact meshgrid(2:3,[2, 0, 4])  -->  [2 3 2 3 2 3; 2 2 0 0 4 4]
        @fact centeredmeshgrid(2:3,2:4)  -->  [0 1 0 1 0 1; -1 -1 0 0 1 1]
    end

    shouldtestcontext("inpointcloud") do
        l = row(linspace(0,1,100))
        o = row(ones(100))
        shape2d = [l o flipdim(l,2) 0*o; 0*o l o flipdim(l,2)]*10.+[2.5;5.5]
        shouldbe2d = zeros(30,30)
        shouldbe2d[3:12,6:15] = 1
        coords = meshgrid(30,30)
        result2d = reshape(map(coords, x->inpointcloud(x, shape2d)), 30,30);
        @fact result2d  -->  shouldbe2d

        plate1 = [meshgrid(1:10,1:10); zeros(1,100)]
        plate2 = [meshgrid(1:10,1:10); 10*ones(1,100)]
        plate3 = plate1[[1, 3, 2],:]
        plate4 = plate2[[1, 3, 2],:]
        plate5 = plate1[[3, 1, 2],:]
        plate6 = plate2[[3, 1, 2],:]
        shape3d = [plate1 plate2 plate3 plate4 plate5 plate6].+[2.5;5.5;5.5]

        shouldbe3d = zeros(30,30,30)
        shouldbe3d[3:12,6:15,6:15] = 1
        coords = meshgrid(30,30,30)
        result3d = reshape(map(coords, x->inpointcloud(x, shape3d)), 30,30,30);
        result3d == shouldbe3d
        @fact sum(result3d-shouldbe3d)  -->  38
    end

    shouldtestcontext("sampler") do
        image = rand(200,100)
        sampler = makeSampler(image, 32)
        coords = asint(50 + 50*randn(2,100))
        if VERSION.minor >= 4
            @fact size(@p map coords sampler) --> (32,32,100)
        else
            @fact size(@p map coords x->sampler(x)) --> (32,32,100)
        end
        out = zeros(1024,1)
        pos = col([50,50])
        @fact size(sample!(out, pos, sampler)) --> (1024,1)

        sampler = makeSampler(image, 32, col = true)
        @fact size(sample(pos, sampler)) --> (1024,1)

        sampler = makeSampler(image, 32, col = true, normmeanstd = true)
        r = sample(pos, sampler)
        @fact abs(mean(r))  < 1e-5 --> true
        @fact abs(std(r)-1) < 1e-5 --> true

        image = [1 2 3 4 5; 10 20 30 40 50; 11 22 33 44 55]
        sampler = makeSampler(image, 3, centered = false)
        @fact sample(col([1,2]), sampler) --> [2 3 4; 20 30 40; 22 33 44]
        sampler = makeSampler(image, 3, centered = true)
        @fact sample(col([1,2]), sampler) --> [1 2 3; 10 20 30; 11 22 33]
        @fact sample(col([1,3]), sampler) --> [2 3 4; 20 30 40; 22 33 44]

        image = ones(20,10)
        sampler = makeSampler(image, 3, centered = false)
        sample(sampler,ones(siz(image)))
        sample(sampler,siz(image))

        sampler = makeSampler(image, 3, centered = true)
        sample(sampler,ones(siz(image)))
        sample(sampler,siz(image))
    end
end

shouldtest("graphics") do
    shouldtestcontext("jetcolormap") do
        r = jetcolormap(10)
        @fact all(r .>= 0) --> true
        @fact all(r .<= 1) --> true
    end

    shouldtestcontext("jetcolors") do
        @fact jetcolors(1:3) --> roughly([0.0 0.5198019801980198 0.5396039603960396; 0.0 1.0 0.0; 0.5396039603960396 0.4801980198019802 0.0])
        @fact jetcolors(1:3,2,4) --> roughly([0.0 0.0 0.5198019801980198; 0.0 0.0 1.0; 0.5396039603960396 0.5396039603960396 0.4801980198019802])
        @fact jetcolorants(1:3) --> Any[RGB{Float64}(0.0,0.0,0.5396039603960396),RGB{Float64}(0.5198019801980198,1.0,0.4801980198019802),RGB{Float64}(0.5396039603960396,0.0,0.0)]
    end

     shouldtestcontext("asimagesc") do
        img = [1 2 3 4 5 6 7 8 9]
        r = asimagescrgb(img)
        @fact size(r,3) --> 3
    end

    shouldtestcontext("blocksvisu") do
        @fact blocksvisu([1 2 3]) --> [1 0 3; 0 0 0; 2 0 0]
        @fact blocksvisu([1 2 3; 1 2 3; 1 2 3; 1 2 3]) --> [1 1 0 3 3; 1 1 0 3 3; 0 0 0 0 0; 2 2 0 0 0; 2 2 0 0 0]
    end

    shouldtestcontext("overlaygradient") do
        img = rand(100,200)
        sp = zeros(size(img))
        sp[25:75, 50:150] = 1
        @fact size(overlaygradient(img, sp)) --> (100,200,3)
    end
 end

shouldtest("numerical") do
    shouldtestcontext("nanfunctions") do
        @fact nanmedian([1 2 3], 1)  -->  [1.0 2.0 3.0]
        @fact isequal(nanmedian([1 2 NaN], 1),  [1.0 2.0 NaN]) --> true
        @fact nanmedian([1 2 3], 2)  -->  row([2.0])
        @fact nanmedian([1 2 NaN], 2)  -->  row([1.5])
        @fact nanmedian([1, 2, NaN])  -->  1.5
    end
    shouldtestcontext("distance") do
        @fact distance(1,1)  -->  0
        @fact distance(1.,1.)  -->  0.
        @fact distance([1,2,3],[2,2,3])  -->  1
        @fact distance([0 1; 0 1],[0 1; 0 1])  -->  [0. sqrt(2); sqrt(2) 0.]
    end
    shouldtestcontext("norms") do
        @fact normsum([1 2 1])  -->  [1/4 1/2 1/4]
        @fact norm01([1 2 3])  -->  [0 1/2 1]
        @fact normeuclid([1 1])  -->  1./[sqrt(2) sqrt(2)]
        @fact normmeanstd([1 2 3])  -->  [-1 0 1]
        @fact normquantile(collect(0:10))  -->  roughly(asfloat64(collect(0:10))/10)
        @fact normquantile(collect(0:2:20))  -->  roughly(asfloat64(collect(0:10))/10)
    end
    shouldtestcontext("extrema") do
        @fact (@p maximum [1 2 3; 4 5 0] (x->x[2]))  -->  [2 5]'
        @fact (@p minimum [1 2 3; 4 5 0] (x->x[2]))  -->  [3 0]'
    end
    shouldtestcontext("valuemap") do
        @fact valuemap([0 1 2 3], [10 20 30])  -->  [0. 10. 20. 30.]
    end
    shouldtestcontext("clamp") do
        @fact clamp(1,1,1)  -->  1
        @fact clamp(0,1,2)  -->  1
        @fact clamp(2,1,2)  -->  2
        @fact clamp([-1 2 3 10; -5 10 11  100],[1 10],[3 15])  -->  [1 2 3 3; 10 10 11 15]
        @fact clamp(Float32[-1 2 3 10; -5 10 11  100],[1 10],[3 15])  -->  [1 2 3 3; 10 10 11 15]
        @fact clamp([1,2],[3,4],[5,6])   -->  [3,4]
        @fact clamp([0 1 11; -5 3 20], rand(10,15))  -->  [1 1 10; 1 3 15]
        @fact clamp([0 1 11; -5 3 20; -4 3 5], rand(10,15,4))  -->  [1 1 10; 1 3 15; 1 3 4]
    end
    shouldtestcontext("rand") do
        @fact size(randsiz(siz(rand(2,3)))) --> (2,3)
        @fact eltype(randsiz([2 3]', Float32)) --> Float32
    end
    shouldtestcontext("plus") do
        @fact (@p plus zeros(2,3) ones(2,1))  -->  ones(2,3)
        @fact (@p minus zeros(2,3) ones(2,1))  -->  -ones(2,3)
        @fact (@p times ones(2,3) zeros(2,1))  -->  zeros(2,3)
        @fact (@p divby ones(2,3) 2*ones(2,1))  -->  0.5*ones(2,3)
    end
end

shouldtest("utils") do
    shouldtestcontext("histdict") do
        d = [Dict(:a => 1), Dict(:a => 1), Dict(:a => 10)]
        r = Dict(1 => 2, 10 => 1)
        @fact histdict(d, :a)  -->  r
        @fact histdict([1,1,10])  -->  r
    end
end

exitwithteststatus()


text:
using Random
import StatsBase
export Daisy, Land, DaisyWorld

mutable struct Daisy <: AbstractAgent
    id::Int
    pos::Dims{2}
    breed::Symbol
    age::Int
    albedo::Float64 # 0-1 fraction
end

mutable struct Land <: AbstractAgent
    id::Int
    pos::Dims{2}
    temperature::Float64
end

const DaisyWorld = ABM{<:GridSpace,Union{Daisy,Land}};

"""
``` julia
daisyworld(;
    griddims = (30, 30),
    max_age = 25,
    init_white = 0.2,
    init_black = 0.2,
    albedo_white = 0.75,
    albedo_black = 0.25,
    surface_albedo = 0.4,
    solar_change = 0.005,
    solar_luminosity = 1.0,
    scenario = :default,
    seed = 165
)
```
Same as in [Daisyworld](@ref).

To access the `Daisy` and `Land` types, simply call
``` julia
using Agents.Models: Daisy, Land
```
"""
function daisyworld(;
    griddims = (30, 30),
    max_age = 25,
    init_white = 0.2, # % cover of the world surface of white breed
    init_black = 0.2, # % cover of the world surface of black breed
    albedo_white = 0.75,
    albedo_black = 0.25,
    surface_albedo = 0.4,
    solar_change = 0.005,
    solar_luminosity = 1.0, # initial luminosity
    scenario = :default,
    seed = 165,
)

    rng = MersenneTwister(seed)
    space = GridSpace(griddims)
    properties = Dict(
        :max_age => max_age,
        :surface_albedo => surface_albedo,
        :solar_luminosity => solar_luminosity,
        :solar_change => solar_change,
        :scenario => scenario,
        :tick => 0,
    )
    ## create a scheduler that only schedules Daisies
    daisysched(model) = [a.id for a in allagents(model) if a isa Daisy]
    model = ABM(
        Union{Daisy,Land},
        space;
        scheduler = daisysched,
        properties,
        rng,
        warn = false,
    )

    ## fill model with `Land`: every grid position has 1 land instance
    fill_space!(Land, model, 0.0) # zero starting temperature

    ## Populate with daisies: each position has only one daisy (black or white)
    grid = collect(positions(model))
    num_positions = prod(griddims)
    white_positions =
        StatsBase.sample(grid, Int(init_white * num_positions); replace = false)
    for wp in white_positions
        wd = Daisy(nextid(model), wp, :white, rand(model.rng, 0:max_age), albedo_white)
        add_agent_pos!(wd, model)
    end
    allowed = setdiff(grid, white_positions)
    black_positions =
        StatsBase.sample(allowed, Int(init_black * num_positions); replace = false)
    for bp in black_positions
        wd = Daisy(nextid(model), bp, :black, rand(model.rng, 0:max_age), albedo_black)
        add_agent_pos!(wd, model)
    end

    return model, daisyworld_agent_step!, daisyworld_model_step!
end

function update_surface_temperature!(pos::Dims{2}, model::DaisyWorld)
    ids = ids_in_position(pos, model)
    ## All grid positions have at least one agent (the land)
    absorbed_luminosity = if length(ids) == 1
        ## Set luminosity via surface albedo
        (1 - model.surface_albedo) * model.solar_luminosity
    else
        ## more than 1 agents: daisy exists
        ## Set luminosity via daisy albedo
        (1 - model[ids[2]].albedo) * model.solar_luminosity
    end
    ## We expect local heating to be 80 ᵒC for an absorbed luminosity of 1,
    ## approximately 30 for 0.5 and approximately -273 for 0.01.
    local_heating = absorbed_luminosity > 0 ? 72 * log(absorbed_luminosity) + 80 : 80
    ## Surface temperature is the average of the current temperature and local heating.
    T0 = model[ids[1]].temperature
    model[ids[1]].temperature = (T0 + local_heating) / 2
end

function diffuse_temperature!(pos::Dims{2}, model::DaisyWorld)
    ratio = get(model.properties, :ratio, 0.5) # diffusion ratio
    ids = nearby_ids(pos, model)
    meantemp = sum(model[i].temperature for i in ids if model[i] isa Land) / 8
    land = model[ids_in_position(pos, model)[1]] # land at current position
    ## Each neighbor land patch is giving up 1/8 of the diffused
    ## amount to each of *its* neighbors
    land.temperature = (1 - ratio) * land.temperature + ratio * meantemp
end

function propagate!(pos::Dims{2}, model::DaisyWorld)
    ids = ids_in_position(pos, model)
    if length(ids) > 1
        daisy = model[ids[2]]
        temperature = model[ids[1]].temperature
        ## Set optimum growth rate to 22.5 ᵒC, with bounds of [5, 40]
        seed_threshold = (0.1457 * temperature - 0.0032 * temperature^2) - 0.6443
        if rand(model.rng) < seed_threshold
            ## Collect all adjacent position that have no daisies
            empty_neighbors = Tuple{Int,Int}[]
            neighbors = nearby_positions(pos, model)
            for n in neighbors
                if length(ids_in_position(n, model)) == 1
                    push!(empty_neighbors, n)
                end
            end
            if !isempty(empty_neighbors)
                ## Seed a new daisy in one of those position
                seeding_place = rand(model.rng, empty_neighbors)
                a = Daisy(nextid(model), seeding_place, daisy.breed, 0, daisy.albedo)
                add_agent_pos!(a, model)
            end
        end
    end
end

function daisyworld_agent_step!(agent::Daisy, model::DaisyWorld)
    agent.age += 1
    agent.age >= model.max_age && kill_agent!(agent, model)
end

daisyworld_agent_step!(agent::Land, model::DaisyWorld) = nothing

function daisyworld_model_step!(model)
    for p in positions(model)
        update_surface_temperature!(p, model)
        diffuse_temperature!(p, model)
        propagate!(p, model)
    end
    model.tick += 1
    solar_activity!(model)
end

function solar_activity!(model::DaisyWorld)
    if model.scenario == :ramp
        if model.tick > 200 && model.tick <= 400
            model.solar_luminosity += model.solar_change
        end
        if model.tick > 500 && model.tick <= 750
            model.solar_luminosity -= model.solar_change / 2
        end
    elseif model.scenario == :change
        model.solar_luminosity += model.solar_change
    end
end


text:
###############################################################
# This macro automatically generates a function body that
# throws an error if the interface function is not implemented.
###############################################################

macro pomdp_func(signature)
    if signature.head == :(=) # in this case a default implementation has already been supplied
        return esc(signature)
    end
    @assert signature.head == :call

    # get the names of all the arguments
    args = [strip_arg(expr) for expr in signature.args[2:end]]

    # get the name of the function
    fname = signature.args[1]

    error_string = "POMDPs.jl: No implementation of $fname for "

    # add each of the arguments to error string
    for (i,a) in enumerate(args)
        error_string *= "$a::\$(typeof($a))"
        if i == length(args)-1
            error_string *= ", and "
        elseif i != length(args)
            error_string *= ", "
        else
            error_string *= "."
        end
    end

    # if you are modifying this and want to debug, it might be helpful to print
    # println(error_string)

    body = Expr(:call, :error, parse("\"$error_string\""))

    return Expr(:function, esc(signature), body)
end

# strip_arg strips anything extra (type annotations, default values, etc) from an argument
# for now this cannot handle keyword arguments (it will throw an error)

strip_arg(arg::Symbol) = arg # once we have a symbol, we have stripped everything, so we can just return it

function strip_arg(arg_expr::Expr)
    if arg_expr.head == :parameters # keywork argument
        error("extract_arg_names can't handle keyword args yet (parsing arg expression $(arg_expr))")
    elseif arg_expr.head == :(::) # argument is type annotated, remove the annotation
        return strip_arg(arg_expr.args[1])
    elseif arg_expr.head == :kw # argument has a default value, remove the default
        return strip_arg(arg_expr.args[1])
    else
        error("strip_arg encountered something unexpected. arg_expr was $(arg_expr)")
    end
end

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
"""
Collection of the methods used for generating the graph representing the
flow model.
"""
module Enumeration

using DocStringExtensions
export Node, Edge, gen_all_edges

"""
A node of the GuillotineModels.Flow model.

$(TYPEDFIELDS)

"""
struct Node{S, N}
	"Unique identifier of the node in the model."
	idx :: N
	"The size of the plate dimension parallel to the cuts."
	par :: S
	"The size of the plate dimension perpendicular to the cuts."
	per :: S
	"The cut orientation: one means vertical cuts; two means horizontal cuts."
	ori :: UInt8
end

"""
An edge of the GuillotineModels.Flow Model.

$(TYPEDFIELDS)

"""
struct Edge{N, E}
  "Unique identifier of the edge in the model."
	indx :: E
	"The node from where the flow comes. Source."
	head :: N
	"The node from where the flow ends. Sink."
	tail :: N
	"The backward edge global index, if it exists, otherwise zero."
	back :: E
end

"""
    gen_rr_fow_edges!(::Type{N}, d :: [D], l :: [S], L :: S) :: ([N], [(S, S)])

!!! **Internal use.**

# Arguments

An integer type capable of representing the nodes `N`. The pieces demand `d`,
their length (or width) `l`, and a plate length (or width) `L`. This method
assumes the selected pieces fit the selected plate (on both dimensions, this
is, in the one that was given and the one that wasn't).

# Return

A tuple of two vectors:
1) a vector `v` of length `L`; `v[y] == 1` iff there is a linear combination
   of the pieces (respecting demand) that gives `y`; otherwise `v[y] == 0`.
2) a vector containing all ordered pairs of linear combinations in which the
   difference between the two is exactly the size of a piece.
"""
function gen_rr_fow_edges!(
	::Type{N}, # eltype of the marked vector
	d :: Vector{D}, l :: Vector{S}, L :: S
) :: Tuple{Vector{N}, Vector{Tuple{S, S}}} where {D, S, N}
	@assert length(d) == length(l)
	n = length(d)
	marked :: Vector{N} = fill(zero(N), L)
	edges = Vector{Tuple{S, S}}() # start point and end point
	iszero(n) && return marked, edges
	# Mark the cuts of the first piece.
	marked[l[1]] = one(N)
	push!(edges, (zero(S), l[1]))
	y = l[1] # y: used to iterate capacity, inherited from knapsack papers
	for _ = 2:min(d[1], L ÷ l[1])
		y_ = y
		push!(edges, (y_, y += l[1]))
		marked[y] = one(N)
	end
	# Mark the cuts of all other pieces.
	for pii = 2:n # PIece Index
		li = l[pii] # length of i
		di = d[pii] # demand of i
		for y = (L - li):-1:1
			if !iszero(marked[y])
				yrli = y + li # yrli: y + repeated lengths of i
				push!(edges, (y, yrli))
				marked[yrli] = one(N)
				for _ = 2:di
					yrli_ = yrli
					yrli += li
					yrli > L && break
					push!(edges, (yrli_, yrli))
					marked[yrli] = one(N)
				end
			end
		end
		# Mark cuts using just multiple copies of the same piece type (could be
		# done using a dummy cut at index zero with value true, but the array has
		# no index zero).
		y = li
		push!(edges, (zero(S), y))
		marked[y] = one(N)
		for _ = 2:min(di, L ÷ li)
			y_ = y
			push!(edges, (y_, y += li))
			marked[y] = one(N)
		end
	end

	marked, edges
end

"""
    reduce_dlw(d::[D], l::[S], w::[S], L::S, W::S) :: ([D], [S], [S], D)

!!! **Internal use.**

Create a copy of `d`, `l`, and `w` with only the pieces that fit the plate
`L` x `W` (i.e., pieces with both dimensions smaller than the corresponding
plate dimension). The last value returned is the number of pieces that pass
this criteria (and, consequently, the common size of the returned vectors).
"""
function reduce_dlw(
	d :: Vector{D}, l :: Vector{S}, w :: Vector{S}, L :: S, W :: S
) :: Tuple{Vector{D}, Vector{S}, Vector{S}, D} where {D, S}
	@assert length(d) == length(l)
	@assert length(d) == length(w)
	N = convert(D, length(d))
	# Remove items that cannot fit into L anyway.
	d_, l_, w_ = similar.((d, l, w))
	j = zero(D)
	for i in one(D):N
		if l[i] <= L && w[i] <= W
			j += one(D)
			d_[j] = d[i]
			l_[j] = l[i]
			w_[j] = w[i]
		end
	end
	resize!(d_, j)
	resize!(l_, j)
	resize!(w_, j)
	d_, l_, w_, j
end

"""
    merge_duplicates(d :: [D], per :: [S]) :: ([D], [S])

!!! **Internal use.**

If `per` has no repeated values just return the parameters. Otherwise,
for each value in `per` with duplicates, keep just one copy of such value
inside `per` and delete the other copies, and the same positions in `d`;
the position in `d` corresponding to the kept copy will have the sum of
the `d` values for all copies deleted and the copy kept.
Sort `per` by increasing value and swap positions in `d` to keep the
correspondence.
"""
function merge_duplicates(
	d :: Vector{D}, per :: Vector{S}
) :: Tuple{Vector{D}, Vector{S}} where {D, S}
	@assert length(d) == length(per)
	isempty(d) && return d, per

	y2d = fill(zero(D), maximum(per))
	has_duplicates = false
	for i in eachindex(d)
		!iszero(y2d[per[i]]) && (has_duplicates = true)
		y2d[per[i]] += d[i]
	end

	!has_duplicates && return d, per

	per_ = Vector{S}()
	d_ = Vector{D}()
	for (y, dy) in enumerate(y2d)
		if !iszero(dy)
			push!(d_, dy)
			push!(per_, y)
		end
	end

	d_, per_
end

"""
    globalize!(see description and source for parameters and return)

!!! **Internal use.**

The model may be seen as constituted of many independent graphs, but as there
is interaction between them (i.e., subplate edges increase the flow of a
backward edge in other graph) their nodes and edges must follow a global
numbering and the backward edges must be known by all of graphs (as subplate
edges need to refer to them as child plates).

This method takes the "local/raw" lists of node (`y2node_idx`) and edge
(`raw_fow_edges`) numbering (that is just the normal cut positions); the index
of the last globalized node (`last_gnode_idx`) and edge (`last_gedge_idx`); a
matrix that gives the piece index given the piece dimensions (`lw2pii`); a
tridimensional array that gives the backward edge index given its plate
dimensions and allowed cut orientation (`ppo2gbedge_idx`); the allowed cut
`ori`entation of the given graph; and the size of the plate dimension parallel
to the allowed cuts (`PAR`).

The method return the "globalized" list of nodes and edges, as the index of
the last globalized node and edge.

The method changes `ppo2gbedge_idx`: if a subplate (forward) edge is created
and the corresponding backward edge was not yet numbered, it is given a
"globalized" index (all backward edges are created at another step at the end).
"""
function globalize!(
	y2node_idx     :: Vector{N},
	raw_fow_edges  :: Vector{Tuple{S, S}}, # raw forward edges
	last_gnode_idx :: N,
	last_gedge_idx :: E,
	lw2pii         :: AbstractArray{D, 2},
	ppo2gbedge_idx :: Array{E, 3},
	ori            :: UInt8,
	PAR            :: S
) :: Tuple{Vector{Node{S, N}}, Vector{Edge{N, E}}, N, E} where {D, S, N, E}
	# The first node of a flow is the point zero, that does not appear in a
	# discretization because the array has no position zero.
	lgni = origin = last_gnode_idx + one(N)
	lgei = last_gedge_idx

	glo_nodes = Vector{Node{S, N}}()
	push!(glo_nodes, Node{S, N}(origin, PAR, zero(S), ori))

	for y in one(S):length(y2node_idx)
		if !iszero(y2node_idx[y])
			y2node_idx[y] = (lgni += one(N))
			push!(glo_nodes, Node{S, N}(lgni, PAR, y, ori))
		end
	end

	# glo_fow_edges: vector in which the values are the globalized forward edges
	# (we call it globalized because indx is the global name of that arc, and
	# both head and tail are global names of nodes).
	glo_fow_edges = Vector{Edge{N, E}}()

	for (y1, y2) in raw_fow_edges
		@assert y1 >= zero(y1) # there are no negative positions
		@assert y2 > y1 # i.e., all edges are forward
		indx = (lgei += one(E))
		if iszero(y1)
			head = origin
		else
			head = y2node_idx[y1]
		end
		tail = y2node_idx[y2]
		# NOTE: the backward edges are not created here. The backward edges are
		# created after all forward edges in the whole model are created.
		# NOTE: the ppo2gbedge_idx is indexed by parallel-perpendicular-orientation
		# but as we check a child plate where the orientation is inversed, then
		# the ori value is inverted, and parallel and perpendicular for the current
		# flow/plate are reversed for the child plate.
		bedge_ppo = (y2 - y1, PAR, 0x03 - ori)
		if iszero(ppo2gbedge_idx[bedge_ppo...])
			ppo2gbedge_idx[bedge_ppo...] = (lgei += one(E))
		end
		back = ppo2gbedge_idx[bedge_ppo...]
		edge = Edge{N, E}(indx, head, tail, back)
		#(iszero(head) || iszero(tail)) && @show edge, @__LINE__
		push!(glo_fow_edges, edge)
		# If this foward edge cuts a plate with the exact same size as some piece
		# then we add the "final plate"/piece edge too.
		yd = y2 - y1
		max_par = size(lw2pii, 1)
		max_per = size(lw2pii, 2)
		#@show PAR, yd, max_par, max_per, ori
		if PAR <= max_par && yd <= max_per && !iszero(lw2pii[PAR, yd])
			indx = (lgei += one(E))
			edge = Edge{N, E}(indx, head, tail, lw2pii[PAR, yd])
			push!(glo_fow_edges, edge)
		end
	end

	glo_nodes, glo_fow_edges, lgni, lgei
end

"""
    gen_u_fow_edges(glo_nodes, y2node_idx, last_gedge_idx, ppo2gbedge_idx)

!!! **Internal use.**

This method generate the unrestricted subplate edges, i.e., the subplate edges
representing a cut that creates a subplate in which the size of the dimension
that is perpendicular to the cut is not the size of a single piece, but of a
linear combination of them. Such edges are necessary to solve the unrestricted
case but not to solve restricted case.

# Arguments

* `glo_nodes::Vector{Node{S, N}}`: A list of globalized nodes from a
  single graph (asserts will fail if nodes are from different graphs).
* `y2node_idx::Vector{N}`: A vector in which `!iszero(y2node_idx[y])`
  only if `y` is the `per` value of some node in `glo_nodes` and zero
  otherwise.
* `last_gedge_idx::E`: The highest global edge identifier already in use.
* `ppo2gbedge_idx::Array{E, 3}`: A table that translates the
  parallel-perpendicular-orientation triple to the global backward
  edge index, if it exists. If there is the need to refer to a circulation
  edge but its global identifier does not yet exists, it is created and
  saved to the table.

# Procedure

Loop through all pairs `i` and `j` (i < j) of node indexes (not
counting a possibly dummy sink node, but counting the source node),
create a edge between the two nodes if `glo_nodes[j].per - glo_nodes[i].per
<= glo_nodes[i].per`.

# Returns

1. The list of new edges, already globalized.
2. The last edge index attributed.
"""
function gen_u_fow_edges(
	glo_nodes :: Vector{Node{S, N}},
	y2node_idx :: Vector{N},
	last_gedge_idx :: E,
	ppo2gbedge_idx :: Array{E, 3}
) :: Tuple{Vector{Edge{N, E}}, E} where {S, N, E}
	u_fow_edges = Vector{Edge{N, E}}()
	isempty(glo_nodes) && return u_fow_edges, last_gedge_idx
	@assert iszero(glo_nodes[1].per)
	@assert isone(length(unique!(map(gn -> gn.par, glo_nodes))))
	@assert isone(length(unique!(map(gn -> gn.ori, glo_nodes))))
	lgei = last_gedge_idx
	for node in @view glo_nodes[2:end-1]
		bedge_ppo = (node.per, node.par, 0x03 - node.ori)
		if iszero(ppo2gbedge_idx[bedge_ppo...])
			ppo2gbedge_idx[bedge_ppo...] = (lgei += one(E))
		end
		back = ppo2gbedge_idx[bedge_ppo...]

		edge = Edge(lgei += one(E), glo_nodes[1].idx, node.idx, back)
		#(iszero(edge.head) || iszero(edge.tail)) && @show edge, @__LINE__
		push!(u_fow_edges, edge)
	end
	for i = 2:(length(glo_nodes)-1)
		for j = (i+1):length(glo_nodes)
			@assert glo_nodes[j].per > glo_nodes[i].per
			per_dist = glo_nodes[j].per - glo_nodes[i].per
			per_dist > glo_nodes[i].per && break
			iszero(y2node_idx[per_dist]) && continue
			bedge_ppo = (
				per_dist,
				glo_nodes[1].par, # all par are the same
				0x03 - glo_nodes[1].ori # all ori are the same
			)
			if iszero(ppo2gbedge_idx[bedge_ppo...])
				ppo2gbedge_idx[bedge_ppo...] = (lgei += one(E))
			end
			back = ppo2gbedge_idx[bedge_ppo...]
			edge = Edge(
				lgei += one(E), glo_nodes[i].idx, glo_nodes[j].idx, back
			)
			#(iszero(edge.head) || iszero(edge.tail)) && @show edge, @__LINE__
			push!(u_fow_edges, edge)
		end
	end
	u_fow_edges, lgei
end

"""
    gen_w_fow_edges(glo_nodes, last_gedge_idx) :: ([Edge], E)

!!! **Internal use.**

Given a set of "globalized" nodes, and the index of the last "globalized"
edge, create the set of corresponding waste edges (already "globalized").

For now, it is not very smart and just connect the second node to the third,
third to the fourth, and so on. Ideally, it should connect every node
(except the first) to the next node *that has a backward edge associated*
(not counting the node in question, clearly, no self-loops allowed).
"""
function gen_w_fow_edges(
	glo_nodes :: Vector{Node{S, N}},
	last_gedge_idx :: E
) :: Tuple{Vector{Edge{N, E}}, E} where {S, N, E}
	w_fow_edges = Vector{Edge{N, E}}()
	lgei = last_gedge_idx
	isempty(glo_nodes) && return w_fow_edges, lgei
	@assert iszero(glo_nodes[1].per)
	@assert isone(length(unique!(map(gn -> gn.par, glo_nodes))))
	@assert isone(length(unique!(map(gn -> gn.ori, glo_nodes))))
	#sink_idx = glo_nodes[end].idx
	# Start from the second position, to avoid wasting the whole plate (if
	# the whole plate is to be wasted, this should happen at higher level).
	# Add waste arcs from every node to the next. Maybe it could be smarter
	# than this, but at least is O(n).
	for i = 2:(length(glo_nodes)-1)
		#for j = (i+1):length(glo_nodes)
		edge = Edge(lgei += 1, glo_nodes[i].idx, glo_nodes[i+1].idx, zero(E))
		#(iszero(edge.head) || iszero(edge.tail)) && @show edge, @__LINE__
		push!(w_fow_edges, edge)
		#end
	end
	w_fow_edges, lgei
end

"""
    gen_closed_flow(...)

!!! **Internal use.**

Given some "global" structures/counters (`last_gnode_idx`, `last_gedge_idx`,
`lw2pii`, `ppo2gbedge_idx`), the piece set (`d`, `par`, `per`) and the plate
dimensions and its allowed cut orientation (`ori`, `PAR`, `PER`), this method
builds the graph that represents all allowed cuts over such plate. The only
changed parameter is `ppo2gbedge_idx`.

# Arguments

* `last_gnode_idx::N`: The highest global node identifier already in use.
* `last_gedge_idx::E`: The highest global edge identifier already in use.
* `lw2pii::AbstractArray{D, 2}`: A convenient table that translates the
  dimensions of a piece to the global piece code if there exists a piece
  that fits the description. It is not changed but a transposition of it is
  made frequently.
* `ppo2gbedge_idx::Array{E, 3}`: A table that translates the
  parallel-perpendicular-orientation triple to the global backward
  edge index, if it exists. If there is the need to refer to a circulation
  edge but its global identifier does not yet exists, it is created and
  saved to the table.
* `d::Vector{D}`: The demand of the pieces.
* `par::Vector{S}`: The size of the pieces in the dimension that is
  parallel to the cuts made in this flow.
* `per::Vector{S}`: The size of the pieces in the dimension that is
  perpendicular to the cuts made in this flow.
* `ori::UInt8`: If it is one, the flow is making vertical cuts, and
  therefore `par` is `l` and `per` is `w`. If it is two, the flow is making
  horizontal cuts and therefore `par` is `w` and `per` is `l`.
* `PAR::S`: The size of the flow/plate in the dimension that is parallel to
  the cuts made.
* `PER::S`: The size of the flow/plate in the dimension that is perpendicular
  to the cuts made.

# Returns

* A dense list of globalized nodes (all nodes in sequence).
* A sparse list of globalized nodes (if position `y` has a corresponding node,
  then `v[y] == globalized_node_id`, otherwise `v[y] == 0`).
* A dense list of the globalized edges.
* The index of the last globalized node (after the procedure).
* The index of the last globalized edge (after the procedure).
"""
function gen_closed_flow(
	last_gnode_idx :: N, last_gedge_idx :: E,
	lw2pii :: AbstractArray{D, 2}, ppo2gbedge_idx :: Array{E, 3},
	d :: Vector{D}, par :: Vector{S}, per :: Vector{S},
	ori :: UInt8, PAR :: S, PER :: S
) :: Tuple{Vector{Node{S, N}}, Vector{N}, Vector{Edge{N, E}}, N, E} where {D, S, N, E}
	# Abbreviate.
	lgni = last_gnode_idx
	lgei = last_gedge_idx

	# Reduce the pieces considered to just the ones that fit the flow/plate.
	d, par, per, n = reduce_dlw(d, par, per, PAR, PER)
	d, per = merge_duplicates(d, per)
	par = nothing # not used anymore after here
	# y2node_idx: the sparse discretization vector.
	#   If y2node_idx[y] is zero, then there is no way to reach that position y
	#   with a piece combination, otherwise, there is a way. A reachable position
	#   is first marked with one, and at the end, with the universal node index.
	# rr_fow_edges: raw restriced forward arcs. The edges (y1, y2)
	#   where y1 < y2, and both are positions inside 0:PER, and y2 - y1 is
	#   always the size of a single piece in the dimension perpendicular to
	#   the cuts.
	y2node_idx, rr_fow_edges = gen_rr_fow_edges!(
		N, d, per, PER
	)
	#@show rr_fow_edges

	dummy_sink = false
	if iszero(y2node_idx[PER])
		dummy_sink = true
		y2node_idx[PER] = one(n)
	end

	# rr_fow_edges: not raw anymore, but yet restricted
	glo_nodes, r_fow_edges, lgni, lgei = globalize!(
		y2node_idx, rr_fow_edges, lgni, lgei, lw2pii, ppo2gbedge_idx, ori, PAR
	)

	w_fow_edges, lgei = gen_w_fow_edges(glo_nodes, lgei)

	# NOTE: symmetries may be broken with extra constraints. If some edge e1 is
	# used, then no edge e2, that is both "larger" than that e1 and that may be
	# reached after traversing e1, may be used.
	u_fow_edges, lgei = gen_u_fow_edges(
		glo_nodes, y2node_idx, lgei, ppo2gbedge_idx
	)
	glo_nodes, y2node_idx, vcat(r_fow_edges, w_fow_edges, u_fow_edges), lgni, lgei
	#=
	glo_nodes, y2node_idx, append!(r_fow_edges, w_fow_edges), lgni, lgei
	=#
end

"""
    gen_nodes_and_edges(N, E, d, l, w, L, W) :: ([Node], [Edge], N, E)

!!! **Internal use.**

Given the integer types used to number nodes and edges, the piece set, and the
original plate length, create all nodes and edges except by some dummy ones
(that are handled specially and, so, should not even be mixed here).

# Notes

* This method is aware the edge identifiers `1` to `n+4` are special.
* This method creates the special edges n+2 and n+4, but not any other in
  the `1:n+4` special identifier interval.

# High-level procedure description

* Compute the discretized lengths (widths).
* For each discretized length (width), create a flow with max width (length).
* The backward edges in not-yet-reached flow B may be referred in
  currently-being-built flow A (a subplate edge from A may enable more flow to
  a backward edge in B). To solve this problem, a structure is passed around,
  it is a matrix of three dimensions, which may indexed by the values that are
  unique for a backward edge (par_dim, cut_ori, per_dim), and stores the
  identifier for the unique backward edge with such attributes. If some forward
  edge created need to refer to a backward edge it is receives its identifier
  at that moment and is always referred by it, the corresponging `Edge` struct
  however, is only created at the final of the procedure, when all backward
  edges are known.
* Backward Edges are similar to waste edges in the fact both of they have the
  `back` field as zero (they never enable another backward edge), however only
  backward edges have tails pointing to the source nodes of flows (and,
  consequently, the only ones with tails before heads).
"""
function gen_all_edges(
	::Type{N}, ::Type{E},
	d :: Vector{D}, l :: Vector{S}, w :: Vector{S},
	L :: S, W :: S
) :: Tuple{Vector{Node{S, N}}, Vector{Edge{N, E}}, N, E} where {D, S, N, E}
	@assert length(d) == length(l)
	@assert length(d) == length(w)
	n = length(d)
	lgni = zero(N)
	lgei = convert(E, n + 4) # the first n + 4 edge indexes are reserved
	edges = Vector{Edge{N, E}}()
	nodes = Vector{Node{S, N}}()
	lw2pii = fill(zero(E), L, W)
	for pii = 1:n
		(l[pii] > L || w[pii] > W) && continue
		if !iszero(lw2pii[l[pii], w[pii]])
			@warn(
				"CAUTION: the instance has two pieces with the exact same length" *
				" and width, the first is piece nº $(lw2pii[l[pii], w[pii]]) and" *
				" the second is nº $(pii), both have length $(l[pii]) and width " *
				"$(w[pii])"
			)
		end
		lw2pii[l[pii], w[pii]] = pii
	end
	max_LW = max(L, W)
	ppo2gbedge_idx = fill(zero(E), max_LW, max_LW, 2)
	origin_vf_by_L = fill(zero(N), L)
	origin_hf_by_W = fill(zero(N), W)
	vflows_by_L = [Vector{N}() for _ = 1:L]
	hflows_by_W = [Vector{N}() for _ = 1:W]

	# Mark the two root backward edges to be created at the end
	ppo2gbedge_idx[L, W, 0x01] = n + 2
	ppo2gbedge_idx[W, L, 0x02] = n + 4

	# Create the cutting of the root plate with the default orientation
	# (vertical).
	vroot_nodes, vflows_by_L[L], vroot_edges, lgni, lgei = gen_closed_flow(
		lgni, lgei, lw2pii, ppo2gbedge_idx, d, l, w, 0x01, L, W
	)
	@assert !isempty(vroot_nodes)
	@assert iszero(vroot_nodes[1].per)
	@assert vroot_nodes[1].par == L
	origin_vf_by_L[L] = vroot_nodes[1].idx
	disc_W = map(n -> n.per, vroot_nodes)
	#@show vroot_edges
	append!(edges, vroot_edges)
	append!(nodes, vroot_nodes)

	#for ppo in CartesianIndices(ppo2gbedge_idx) # DEBUG PRINT FOR
	#  iszero(ppo2gbedge_idx[ppo]) && continue
	#  @show ppo
	#  @show ppo2gbedge_idx[ppo]
	#end

	# Create the forward edge that inverts the default orientation
	# of the root plate.
	push!(edges, Edge{N, E}(
		n + 3, origin_vf_by_L[L], lgni, n + 4
	))
	#@show last(edges), @__LINE__

	hroot_nodes, hflows_by_W[W], hroot_edges, lgni, lgei = gen_closed_flow(
		lgni, lgei, lw2pii', ppo2gbedge_idx, d, w, l, 0x02, W, L
	)
	origin_hf_by_W[W] = hroot_nodes[1].idx
	disc_L = map(n -> n.per, hroot_nodes)
	#@show hroot_edges
	append!(edges, hroot_edges)
	append!(nodes, hroot_nodes)

	#for ppo in CartesianIndices(ppo2gbedge_idx) # DEBUG PRINT FOR
	#  iszero(ppo2gbedge_idx[ppo]) && continue
	#  @show ppo
	#  @show ppo2gbedge_idx[ppo]
	#end

	#unique_lw
	#@show disc_L
	#@show disc_W
	for L_ in @view disc_L[2:end-1]#setdiff(unique(l), [L])#
		vnodes, vflows_by_L[L_], vedges, lgni, lgei = gen_closed_flow(
			lgni, lgei, lw2pii, ppo2gbedge_idx, d, l, w, 0x01, L_, W
		)
		origin_vf_by_L[L_] = vnodes[1].idx
		append!(edges, vedges)
		append!(nodes, vnodes)
	end
	for W_ in @view disc_W[2:end-1]#setdiff(unique(w), [W])#
		hnodes, hflows_by_W[W_], hedges, lgni, lgei = gen_closed_flow(
			lgni, lgei, lw2pii', ppo2gbedge_idx, d, w, l, 0x02, W_, L
		)
		origin_hf_by_W[W_] = hnodes[1].idx
		append!(edges, hedges)
		append!(nodes, hnodes)
	end
	#@show vflows_by_L
	#@show hflows_by_W

	for ppo in CartesianIndices(ppo2gbedge_idx)
		pari, peri, orii = ppo[1], ppo[2], ppo[3]
		iszero(ppo2gbedge_idx[ppo]) && continue
		#@show ppo
		if orii == 0x01
			while !iszero(peri) && iszero(vflows_by_L[pari][peri]); peri -= 1; end
			@assert !iszero(peri)
			push!(edges, Edge{N, E}(
				ppo2gbedge_idx[ppo], vflows_by_L[pari][peri], origin_vf_by_L[pari],
				zero(N)
			))
			edge = last(edges)
			#(iszero(edge.head) || iszero(edge.tail)) && @show edge, ppo, peri, @__LINE__
		else
			@assert orii == 0x02
			while !iszero(peri) && iszero(hflows_by_W[pari][peri]); peri -= 1; end
			push!(edges, Edge{N, E}(
				ppo2gbedge_idx[ppo], hflows_by_W[pari][peri], origin_hf_by_W[pari],
				zero(N)
			))
			edge = last(edges)
			#(iszero(edge.head) || iszero(edge.tail)) && @show edge, ppo, peri, @__LINE__
		end
	end

	nodes, edges, lgni, lgei
end

end # module


text:
using StatsBase
using Base.Test

a = rand(5, 6)
b = rand(5, 6)

@test_approx_eq sqL2dist(a, b) sum(abs2(a - b))
@test_approx_eq L2dist(a, b)   sqrt(sqL2dist(a, b))
@test_approx_eq L1dist(a, b)   sum(abs(a - b))
@test_approx_eq Linfdist(a, b) maximum(abs(a - b))

@test_approx_eq gkldiv(a, b) sum(a .* log(a ./ b) - a + b)

@test_approx_eq meanad(a, b) mean(abs(a - b))
@test_approx_eq maxad(a, b)  maximum(abs(a - b))
@test_approx_eq msd(a, b)    mean(abs2(a - b))
@test_approx_eq rmsd(a, b)   sqrt(msd(a, b))
@test_approx_eq rmsd(a, b; normalize=true)  rmsd(a, b) / (maximum(a) - minimum(a))
@test_approx_eq psnr(a, b, 2)  10 * log10(4 / msd(a, b))


text:
using Compat: eachcol
using Crystallography: Cell, ReciprocalPoint, MonkhorstPackGrid
using Functors: fmap
using StaticArrays: SVector, SMatrix

import ..Inputs: optionpool

export AtomicSpecies,
    AtomicSpeciesCard,
    AtomicPosition,
    AtomicPositionsCard,
    CellParametersCard,
    AtomicForce,
    AtomicForcesCard,
    KPointsCard,
    KMeshCard,
    GammaPointCard,
    SpecialPointsCard
export optconvert, optionpool, eachatom, getpotentials

"""
    AtomicSpecies(atom::Union{AbstractChar,String}, mass::Float64, pseudopot::String)
    AtomicSpecies(x::AtomicPosition, mass, pseudopot)

Represent each line of the `ATOMIC_SPECIES` card in QE.

The `atom` field accepts at most 3 characters.

# Examples
```julia
julia> AtomicSpecies("C1", 12, "C.pbe-n-kjpaw_psl.1.0.0.UPF")
AtomicSpecies("C1", 12.0, "C.pbe-n-kjpaw_psl.1.0.0.UPF")

julia> AtomicSpecies(
           AtomicPosition('S', [0.500000000, 0.288675130, 1.974192764]),
           32.066,
           "S.pz-n-rrkjus_psl.0.1.UPF",
       )
AtomicSpecies("S", 32.066, "S.pz-n-rrkjus_psl.0.1.UPF")
```
"""
struct AtomicSpecies
    "Label of the atom. Max total length cannot exceed 3 characters."
    atom::String
    """
    Mass of the atomic species in atomic unit.

    Used only when performing molecular dynamics (MD) run
    or structural optimization runs using damped MD.
    Not actually used in all other cases (but stored
    in data files, so phonon calculations will use
    these values unless other values are provided).
    """
    mass::Float64
    """
    File containing pseudopotential for this species.

    See also: [`pseudoformat`](@ref)
    """
    pseudopot::String
    function AtomicSpecies(atom::Union{AbstractChar,AbstractString}, mass, pseudopot)
        @assert length(atom) <= 3 "`atom` can have at most 3 characters!"
        return new(string(atom), mass, pseudopot)
    end
end

"""
    pseudoformat(data::AtomicSpecies)

Return the pseudopotential format of the `AtomicSpecies`.
"""
# pseudoformat(data::AtomicSpecies) = pseudoformat(data.pseudopot)

"""
    AtomicSpeciesCard <: Card

Represent the `ATOMIC_SPECIES` card in QE. It does not have an "option".
"""
struct AtomicSpeciesCard <: Card
    data::Vector{AtomicSpecies}
end
@batteries AtomicSpeciesCard eq = true hash = true

"""
    getpotentials(card::AtomicSpeciesCard)

Get the pseudopotential names from an `AtomicSpeciesCard`.
"""
function getpotentials(card::AtomicSpeciesCard)
    return map(card.data) do atomic_species
        atomic_species.pseudopot
    end
end

"""
    AtomicPosition(atom::Union{AbstractChar,String}, pos::Vector{Float64}[, if_pos::Vector{Int}])
    AtomicPosition(x::AtomicSpecies, pos, if_pos)

Represent each line of the `ATOMIC_POSITIONS` card in QE.

The `atom` field accepts at most 3 characters.

# Examples
```julia
julia> AtomicPosition('O', [0, 0, 0])
AtomicPosition("O", [0.0, 0.0, 0.0], Bool[1, 1, 1])

julia> AtomicPosition(
           AtomicSpecies('S', 32.066, "S.pz-n-rrkjus_psl.0.1.UPF"),
           [0.500000000, 0.288675130, 1.974192764],
       )
AtomicPosition("S", [0.5, 0.28867513, 1.974192764], Bool[1, 1, 1])
```
"""
struct AtomicPosition
    "Label of the atom as specified in `AtomicSpecies`."
    atom::String
    "Atomic positions. A three-element vector of floats."
    pos::SVector{3,Float64}
    """
    Component `i` of the force for this atom is multiplied by `if_pos(i)`,
    which must be either `0` or `1`.  Used to keep selected atoms and/or
    selected components fixed in MD dynamics or structural optimization run.

    With `crystal_sg` atomic coordinates the constraints are copied in all equivalent
    atoms.
    """
    if_pos::SVector{3,Bool}
    function AtomicPosition(atom::Union{AbstractChar,AbstractString}, pos, if_pos)
        @assert length(atom) <= 3 "`atom` can have at most 3 characters!"
        return new(string(atom), pos, if_pos)
    end
end
AtomicPosition(atom, pos) = AtomicPosition(atom, pos, trues(3))
AtomicPosition(x::AtomicSpecies, pos, if_pos) = AtomicPosition(x.atom, pos, if_pos)
# Introudce mutual constructors since they share the same atoms.
AtomicSpecies(x::AtomicPosition, mass, pseudopot) = AtomicSpecies(x.atom, mass, pseudopot)

"""
    AtomicPositionsCard <: Card

Represent the `ATOMIC_POSITIONS` card in QE.

# Arguments
- `data::AbstractVector{AtomicPosition}`: A vector containing `AtomicPosition`s.
- `option::String="alat"`: allowed values are: "alat", "bohr", "angstrom", "crystal", and "crystal_sg".
"""
struct AtomicPositionsCard <: Card
    data::Vector{AtomicPosition}
    option::String
    function AtomicPositionsCard(data, option = "alat")
        @assert option in optionpool(AtomicPositionsCard)
        return new(data, option)
    end
end
AtomicPositionsCard(cell::Cell, option) =
    AtomicPositionsCard(map(cell.types, eachcol(cell.positions)) do atom, position
            AtomicPosition(string(atom), position)
        end, option)
@batteries AtomicPositionsCard eq = true hash = true

"Represent the abstraction of `CELL_PARAMETERS` and `REF_CELL_PARAMETERS` cards in QE."
abstract type AbstractCellParametersCard <: Card end

"""
    CellParametersCard{T<:Real} <: AbstractCellParametersCard
    CellParametersCard(data::AbstractMatrix, option::String)

Represent the `CELL_PARAMETERS` cards in `PWscf` and `CP` packages.
"""
struct CellParametersCard <: AbstractCellParametersCard
    data::SMatrix{3,3,Float64}
    option::String
    function CellParametersCard(data, option = "alat")
        @assert option in optionpool(CellParametersCard)
        return new(data, option)
    end
end
CellParametersCard(lattice::Lattice{<:Real}, option) =
    CellParametersCard(transpose(lattice.data), option)
CellParametersCard(lattice::Lattice{<:Length}) =
    CellParametersCard(fmap(x -> ustrip(u"bohr", x), lattice), "bohr")
CellParametersCard(cell::Cell, option) = CellParametersCard(transpose(cell.lattice), option)

struct AtomicForce
    atom::String
    force::SVector{3,Float64}
    function AtomicForce(atom::Union{AbstractChar,AbstractString}, force)
        @assert length(atom) <= 3 "`atom` can have at most 3 characters!"
        return new(string(atom), force)
    end
end

struct AtomicForcesCard <: Card
    data::Vector{AtomicForce}
end
@batteries AtomicForcesCard eq = true hash = true

# See https://github.com/JuliaCollections/IterTools.jl/blob/0ecaa88/src/IterTools.jl#L1008-L1032 & https://github.com/JuliaLang/julia/blob/de3a70a/base/io.jl#L971-L1054
struct EachAtom{T}
    card::T
end

eachatom(card::Union{AtomicSpeciesCard,AtomicPositionsCard,AtomicForcesCard}) =
    EachAtom(card)

Base.length(iter::EachAtom) = length(iter.card.data)

Base.iterate(iter::EachAtom, state = 1) =
    state > length(iter) ? nothing : (iter.card.data[state], state + 1)

Base.eltype(iter::EachAtom) = eltype(iter.card.data)

"""
    optconvert(new_option::AbstractString, card::AbstractCellParametersCard)

Convert the option of an `AbstractCellParametersCard` from "bohr" to "angstrom", or its reverse.

!!! warning
    It does not support conversion between `"alat"` and the others.
"""
function optconvert(new_option::AbstractString, card::AbstractCellParametersCard)
    old_option = optionof(card)
    if new_option == old_option
        return card  # No conversion is needed
    else
        typeof(card)(if (old_option => new_option) == ("bohr" => "angstrom")
            return @. ustrip(u"angstrom", card.data * u"bohr")
        elseif (old_option => new_option) == ("angstrom" => "bohr")
            return @. ustrip(u"bohr", card.data * u"angstrom")
        else
            error("unknown conversion rule $(old_option => new_option)!")
        end)
    end
end # function optconvert

abstract type KPointsCard <: Card end

struct KMeshCard <: KPointsCard
    data::MonkhorstPackGrid
end

struct GammaPointCard <: KPointsCard end

"""
    SpecialKPointsCard(data, option)

Represent the `K_POINTS` card in QE.

# Arguments
- `data::Union{MonkhorstPackGrid,GammaPoint,AbstractVector{SpecialKPoint}}`: A Γ point, a Monkhorst--Pack grid or a vector containing `SpecialKPoint`s.
- `option::String="tpiba"`: allowed values are: "tpiba", "automatic", "crystal", "gamma", "tpiba_b", "crystal_b", "tpiba_c" and "crystal_c".
"""
struct SpecialPointsCard <: KPointsCard
    data::Vector{ReciprocalPoint}
    option::String
    function SpecialPointsCard(data, option = "tpiba")
        @assert option in optionpool(SpecialPointsCard)
        return new(data, option)
    end
end
function SpecialPointsCard(data::AbstractMatrix, option = "tpiba")
    @assert size(data, 2) == 4
    return SpecialPointsCard(map(x -> ReciprocalPoint(x...), eachrow(data)), option)
end
@batteries SpecialPointsCard eq = true hash = true

optionof(::KMeshCard) = "automatic"
optionof(::GammaPointCard) = "gamma"

optionpool(::Type{AtomicPositionsCard}) =
    ("alat", "bohr", "angstrom", "crystal", "crystal_sg")
optionpool(::Type{CellParametersCard}) = ("alat", "bohr", "angstrom")
optionpool(::Type{KMeshCard}) = ("automatic",)
optionpool(::Type{GammaPointCard}) = ("gamma",)
optionpool(::Type{SpecialPointsCard}) =
    ("tpiba", "crystal", "tpiba_b", "crystal_b", "tpiba_c", "crystal_c")

groupname(::Type{AtomicSpeciesCard}) = "ATOMIC_SPECIES"
groupname(::Type{AtomicPositionsCard}) = "ATOMIC_POSITIONS"
groupname(::Type{CellParametersCard}) = "CELL_PARAMETERS"
groupname(::Type{<:KPointsCard}) = "K_POINTS"


text:
function meanvar_seq(xs)
    m = mean(xs)
    return m, var(xs, mean = m)
end

function meanvar_threads(xs)
    m = ThreadsFolds.mean(xs)
    s = ThreadsFolds.sum(x -> (x - m)^2, xs)
    return m, s / (length(xs) - 1)
end

function meanvar_tapir(xs)
    m = TapirFolds.mean(xs)
    s = TapirFolds.sum(x -> (x - m)^2, xs)
    return m, s / (length(xs) - 1)
end


text:
# t3r3: fixed
Fixed(body1::AbstractBody{T}, body2; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T =
    Translational3{T}(body1, body2, p1 = p1, p2 = p2), Rotational3{T}(body1, body2, qoffset = qoffset)

# t2r3: prismatic
Prismatic(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T =
    Translational2{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational3{T}(body1, body2, qoffset = qoffset)

# t1r3: planar
Planar(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T =
    Translational1{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational3{T}(body1, body2, qoffset = qoffset)

# t0r3: fixed orientation (chicken's head)
FixedOrientation(body1::AbstractBody{T}, body2; qoffset = one(UnitQuaternion{T})) where T =
    Translational0{T}(body1, body2), Rotational3{T}(body1, body2, qoffset = qoffset)

# t3r2: Revolute
Revolute(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T =
    Translational3{T}(body1, body2, p1 = p1, p2 = p2), Rotational2{T}(body1, body2, axis = axis, qoffset = qoffset)

# t2r2: Cylindrical
Cylindrical(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T =
    Translational2{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational2{T}(body1, body2, axis = axis, qoffset = qoffset)

# t1r2: Planar with rotation along axis (turtle bot)
PlanarAxis(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3), qoffset = one(UnitQuaternion{T})) where T =
    Translational1{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational2{T}(body1, body2, axis = axis, qoffset = qoffset)

# t0r2: something like a gyro?
# ?

# What is r1?
# t3r1: universal?
# ?
# t2r1: ?
# ?
# t1r1: ?
# ?
# t0r1: ?
# ?

# t3r0: Spherical
Spherical(body1::AbstractBody{T}, body2; p1 = szeros(T, 3), p2 = szeros(T, 3)) where T =
    Translational3{T}(body1, body2, p1 = p1, p2 = p2), Rotational0{T}(body1, body2)

# t2r0: Point-on-Line
CylindricalFree(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3)) where T =
    Translational2{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational0{T}(body1, body2)

# t1r0: Planar with arbitrary rotation
PlanarFree(body1::AbstractBody{T}, body2, axis; p1 = szeros(T, 3), p2 = szeros(T, 3)) where T =
    Translational1{T}(body1, body2, p1 = p1, p2 = p2, axis = axis), Rotational0{T}(body1, body2)

# t0r0: Free
Floating(body1::AbstractBody{T}, body2) where T =
    Translational0{T}(body1, body2), Rotational0{T}(body1, body2)


text:
"""
# Coping with NP-completeness

- Intelligent exhaustive search
    - Backtracking
    - Branch-and-bound
- Approximation algorithms
    - Vertex cover
    - Clustering
    - TSP
    - General TSP
    - Knapsack
    - The approximability hierarchy
- Local search heuristics
    - Traveling salesman, once more
    - Graph partitioning
    - Dealing with local optima
    - Randomization and restarts
    - Simulated annealing
"""
module Chapter9

end # module


text:
#Adding Julia Packages to Jupyter

Pkg.add("DataFrames")
Pkg.add("RDatasets")
quit();

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
@testset "Data Access" begin

# ISO 4217
@testset "ISO 4217" begin
    @testset "Currency Name" begin
        @test currencyinfo(:USD) == "US Dollar"
        @test currencyinfo(20EUR) == "Euro"
        @test currencyinfo(Monetary{:BND, BigInt}) == "Brunei Dollar"
    end

    @testset "Numeric Code" begin
        @test iso4217num(:BZD) == 84
        @test iso4217num(:USD) == 840
        @test iso4217num(:XPF) == 953
        @test iso4217num(Monetary{:XAF}) == 950
        @test iso4217num(Monetary{:XOF, BigInt}) == 952
    end

    @testset "Alphabetic Code" begin
        @test iso4217alpha(:DJF) == "DJF"
        @test iso4217alpha(Monetary{:BND, BigInt}) == "BND"
    end
end

@testset "Symbols" begin
    @test shortsymbol(:USD) == "\$"
    @test shortsymbol(:CAD) == "\$"
    @test shortsymbol(100EUR) == "€"
    @test shortsymbol(Monetary{:GBP}) == "£"
    @test shortsymbol(Monetary{:AUD, BigInt, 4}) == "\$"

    @test longsymbol(:USD) == "US\$"
    @test longsymbol(:CAD) == "C\$"
    @test longsymbol(100EUR) == "€"
    @test longsymbol(Monetary{:GBP}) == "GB£"
    @test longsymbol(Monetary{:AUD, BigInt, 4}) == "A\$"
end

# Currency
@testset "currency()" begin
    @test currency(2GBP) == :GBP
    @test currency(zero(Monetary{:CNY})) == :CNY
end

# Default precision
@testset "decimals()" begin
    @test decimals(:USD) == 2
    @test decimals(:JPY) == 0
    @test decimals(Monetary{:USD, BigInt}) == 2
    @test decimals(Monetary{:JPY, Int, 4}) == 4
    @test decimals(:XAU) == -1
    @test decimals(Monetary(:XAU; precision=4)) == 4
    @test decimals(Monetary{:XAU, Int, 3}) == 3
end

# Decimals
@testset ".val" begin
    @test (-USD).val == -1
    @test (100JPY).val == 100
    @test (100USD).val == 100
    @test Monetary(:USD, 25000; precision=3).val == 25
    @test Monetary(:USD; precision=6).val == 1
    @test zero(Monetary{:USD, Int, 8}).val == 0
end

end  # testset Data Access


text:
module Debugger

using Markdown
using Base.Meta: isexpr
using REPL
using REPL.LineEdit
using REPL.REPLCompletions

using JuliaInterpreter: JuliaInterpreter, Frame, @lookup, FrameCode, BreakpointRef, debug_command, leaf, root

using JuliaInterpreter: pc_expr, moduleof, linenumber, extract_args,
                        root, caller, whereis, get_return, nstatements

const SEARCH_PATH = []
function __init__()
    append!(SEARCH_PATH,[joinpath(Sys.BINDIR,"../share/julia/base/"),
            joinpath(Sys.BINDIR,"../include/")])
    return nothing
end

export @enter, @run, break_on_error

include("LineNumbers.jl")
using .LineNumbers: SourceFile, compute_line

mutable struct DebuggerState
    frame::Union{Nothing, Frame}
    level::Int
    broke_on_error::Bool
    repl
    terminal
    main_mode
    julia_prompt::Ref{LineEdit.Prompt}
    standard_keymap
    overall_result
end
DebuggerState(stack, repl, terminal) = DebuggerState(stack, 1, false, repl, terminal, nothing, Ref{LineEdit.Prompt}(), nothing, nothing)
DebuggerState(stack, repl) = DebuggerState(stack, repl, nothing)

function active_frame(state)
    frame = state.frame
    for i in 1:(state.level - 1)
        frame = caller(frame)
    end
    @assert frame !== nothing
    return frame
end

break_on_error(v::Bool) = JuliaInterpreter.break_on_error[] = v

include("locationinfo.jl")
include("repl.jl")
include("commands.jl")
include("printing.jl")

function _make_frame(mod, arg)
    args = try
        extract_args(mod, arg)
    catch e
        return :(throw($e))
    end
    quote
        theargs = $(esc(args))
        frame = JuliaInterpreter.enter_call_expr(Expr(:call,theargs...))
        frame = JuliaInterpreter.maybe_step_through_wrapper!(frame)
        JuliaInterpreter.maybe_next_call!(frame)
        frame
    end
end

macro make_frame(arg)
    _make_frame(__module__, arg)
end

macro enter(arg)
    quote
        let frame = $(_make_frame(__module__,arg))
            RunDebugger(frame)
        end
    end
end

macro run(arg)
    quote
        let frame = $(_make_frame(__module__,arg))
            RunDebugger(frame; initial_continue=true)
        end
    end
end

end # module