text:
using Test
using PlayingCards
using TexasHoldem
TH = TexasHoldem

@testset "Game: show" begin
    players = ntuple(2) do i
        Player(Bot5050(), i)
    end
    game = Game(players)
    sprint(show, game)

    game.table.state = PreFlop()
    sprint(show, game)
end

@testset "Game: pre-dealt deck" begin
    deck = ordered_deck()
    shuffle!(deck)
    players = ntuple(3) do i
        Player(Bot5050(), i, pop!(deck, 2))
    end
    table = Table(;players=players,deck=deck,cards=pop!(deck, 5))
    game = Game(players;deck=deck,table=table)
end

@testset "Game: contrived game" begin
    players = ntuple(3) do i
        Player(Bot5050(), i)
    end
    game = Game(players)
    players = TH.players_at_table(game)
    TH.deal!(game.table, TH.blinds(game.table))
    # Round 1
    check!(game, players[1])
    check!(game, players[2])
    fold!(game, players[3])

    # Round 2
    check!(game, players[1])
    check!(game, players[2])

    # Round 3
    raise_to!(game, players[1], 10)
    @test TH.checked(players[1]) == false
    call!(game, players[2])

    # Round 4
    raise_to!(game, players[1], 20)
    @test TH.checked(players[1]) == false
    fold!(game, players[2])

    # All-in cases
    players = ntuple(3) do i
        Player(Bot5050(), i)
    end
    game = Game(players)
    players = TH.players_at_table(game)
    TH.deal!(game.table, TH.blinds(game.table))
    # Round 1
    check!(game, players[1])
    check!(game, players[2])
    fold!(game, players[3])

    # Round 2
    raise_to!(game, players[1], players[1].bank_roll)
    @test TH.checked(players[1]) == false
    call!(game, players[2])

    @test_throws AssertionError raise_to!(game, players[1], 10000) # raise exceeds bank roll!
end


text:
# exported hyperbolic functions

over_sch(::Type{Float64}) = 710.0
over_sch(::Type{Float32}) = 89f0

"""
    sinh(x)

Compute hyperbolic sine of `x`.
"""
function sinh(x::T) where {T<:Union{Float32,Float64}}
    u = abs(x)
    d = expk2(Double(u))
    d = dsub(d, drec(d))
    u = T(d) * T(0.5)
    u = abs(x) > over_sch(T) ? T(Inf) : u
    u = isnan(u) ? T(Inf) : u
    u = flipsign(u, x)
    u = isnan(x) ? T(NaN) : u
    return u
end

"""
    cosh(x)

Compute hyperbolic cosine of `x`.
"""
function cosh(x::T) where {T<:Union{Float32,Float64}}
    u = abs(x)
    d = expk2(Double(u))
    d = dadd(d, drec(d))
    u = T(d) * T(0.5)
    u = abs(x) > over_sch(T) ? T(Inf) : u
    u = isnan(u) ? T(Inf) : u
    u = isnan(x) ? T(NaN) : u
    return u
end

over_th(::Type{Float64}) = 18.714973875
over_th(::Type{Float32}) = 18.714973875f0

"""
    tanh(x)

Compute hyperbolic tangent of `x`.
"""
function tanh(x::T) where {T<:Union{Float32,Float64}}
    u = abs(x)
    d = expk2(Double(u))
    e = drec(d)
    d = ddiv(dsub(d, e), dadd(d, e))
    u = T(d)
    u = abs(x) > over_th(T) ? T(1.0) : u
    u = isnan(u) ? T(1) : u
    u = flipsign(u, x)
    u = isnan(x) ? T(NaN) : u
    return u
end

"""
    asinh(x)

Compute the inverse hyperbolic sine of `x`.
"""
function asinh(x::T) where {T<:Union{Float32,Float64}}
    y = abs(x)

    d = y > 1 ? drec(x) : Double(y, T(0.0))
    d = dsqrt(dadd2(dsqu(d), T(1.0)))
    d = y > 1 ? dmul(d, y) : d

    d = logk2(dnormalize(dadd(d, x)))
    y = T(d)

    y = (abs(x) > SQRT_MAX(T) || isnan(y)) ? flipsign(T(Inf), x) : y
    y = isnan(x) ? T(NaN) : y
    y = isnegzero(x) ? T(-0.0) : y

    return y
end

"""
    acosh(x)

Compute the inverse hyperbolic cosine of `x`.
"""
function acosh(x::T) where {T<:Union{Float32,Float64}}
    d = logk2(dadd2(dmul(dsqrt(dadd2(x, T(1.0))), dsqrt(dsub2(x, T(1.0)))), x))
    y = T(d)

    y = (x > SQRT_MAX(T) || isnan(y)) ? T(Inf) : y
    y = x == T(1.0) ? T(0.0) : y
    y = x < T(1.0) ? T(NaN) : y
    y = isnan(x) ? T(NaN) : y

    return y
end

"""
    atanh(x)

Compute the inverse hyperbolic tangent of `x`.
"""
function atanh(x::T) where {T<:Union{Float32,Float64}}
    u = abs(x)
    d = logk2(ddiv(dadd2(T(1.0), u), dsub2(T(1.0), u)))
    u = u > T(1.0) ? T(NaN) : (u == T(1.0) ? T(Inf) : T(d) * T(0.5))

    u = isinf(x) || isnan(u) ? T(NaN) : u
    u = flipsign(u, x)
    u = isnan(x) ? T(NaN) : u

    return u
end

ext: jl


text:
using Test, PseudoArcLengthContinuation, LinearAlgebra, Setfield
const PALC = PseudoArcLengthContinuation

source_term(x; a = 0.5, b = 0.01) = 1 + (x + a*x^2)/(1 + b*x^2)
dsource_term(x; a = 0.5, b = 0.01) = (1-b*x^2+2*a*x)/(1+b*x^2)^2
d2source_term(x; a = 0.5, b = 0.01) = -(2*(-a+3*a*b*x^2+3*b*x-b^2*x^3))/(1+b*x^2)^3

println("\n\n\n--> Test Fold continuation")

function F_chan(x, α, β = 0.)
	f = similar(x)
	N = length(x)
	f[1] = x[1] - β
	f[N] = x[N] - β
	for i=2:N-1
		f[i] = (x[i-1] - 2 * x[i] + x[i+1]) * (N-1)^2 + α * source_term(x[i], b = β)
	end
	return f
end

function Jac_mat(u, α, β = 0.)
	N = length(u)
	J = zeros(N,N)
	J[1,1] = 1.0
	J[n,n] = 1.0
	for i = 2:N-1
		J[i,i-1] = (N-1)^2
		J[i,i+1] = J[i,i-1]
		J[i,i] = -2 * J[i,i-1] + α * dsource_term(u[i],b = β)
	end
	return J
end

Jac_fd(u0, α, β) = PALC.finiteDifferences(u->F_chan(u,α, β),u0)

# not really precise Finite Differences
n = 101
sol = rand(n)
sol[end] = sol[1]
J_fold_fd  = Jac_fd(sol,3,0.01)
J_fold_exp = Jac_mat(sol,3,0.01)
@test (J_fold_exp - J_fold_fd) |> x->norm(x,Inf64) < 1e-2

(J_fold_exp - J_fold_fd)[10,10]

n = 101
	a = 3.3
	sol = [(i-1)*(n-i)/n^2+0.1 for i=1:n]
	opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = false)
	# ca fait dans les 69.95k Allocations
	out, hist, flag = @time PALC.newton(
							x -> F_chan(x, a, 0.01),
							x -> Jac_mat(x, a, 0.01),
							sol,
							opt_newton, normN = x->norm(x,Inf64))

# test with secant continuation
opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, pMax = 4.1, maxSteps = 150, newtonOptions = opt_newton)
	br, u1 = @time PALC.continuation(
				(x,p) -> F_chan(x,p, 0.01),
				(x,p) -> (Jac_mat(x,p, 0.01)),
				printSolution = (x,p)->norm(x,Inf64),
				out, a, opts_br0, plot = false, verbosity = 0)

# test with Bordered tangent continuation
br_tg, u1 = @time PALC.continuation(
				(x,p) -> F_chan(x,p, 0.01),
				(x,p) -> (Jac_mat(x,p, 0.01)),
				printSolution = (x,p)->norm(x,Inf64),
				out,a,opts_br0,plot = false, verbosity = 0, tangentAlgo = PALC.BorderedPred())

# test with natural continuation
opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.05, ds= 0.01, pMax = 4.1, newtonOptions = opt_newton)
	br_nat, u1 = @time PALC.continuation(
				(x,p) -> F_chan(x,p, 0.01),
				(x,p) -> (Jac_mat(x,p, 0.01)),
				printSolution = (x,p)->norm(x,Inf64),
				out,0.,opts_br0,plot = false, verbosity = 0, tangentAlgo = PALC.NaturalPred())

# idem with Matrix-Free solver
function dF_chan(x, dx, α, β = 0.)
	out = similar(x)
	n = length(x)
	out[1] = dx[1]
	out[n] = dx[n]
	for i=2:n-1
		out[i] = (dx[i-1] - 2 * dx[i] + dx[i+1]) * (n-1)^2 + α * dsource_term(x[i], b = β) * dx[i]
	end
	return out
end

ls = PALC.GMRESKrylovKit(dim = 100)
	opt_newton_mf = PALC.NewtonPar(tol = 1e-11, verbose = true, linsolver = ls, eigsolver = DefaultEig())
	out_mf, hist, flag = @time PALC.newton(
		x -> F_chan(x, a, 0.01),
		x -> (dx -> dF_chan(x, dx, a, 0.01)),
		sol,
		opt_newton_mf)

opts_cont_mf  = PALC.ContinuationPar(dsmin = 0.01, dsmax = 0.1, ds= 0.01, pMax = 4.1, nev = 5, newtonOptions = setproperties(opt_newton_mf;maxIter = 70, verbose = false, tol = 1e-8), maxSteps = 150)

brmf, u1 = @time PALC.continuation(
		(x, p) -> F_chan(x, p, 0.01),
		(x, p) -> (dx -> dF_chan(x, dx, p, 0.01)),
		out, a, opts_cont_mf, verbosity = 0)

brmf, u1 = @time PALC.continuation(
	(x, p) -> F_chan(x, p, 0.01),
	(x, p) -> (dx -> dF_chan(x, dx, p, 0.01)),
	out, a, opts_cont_mf,
	tangentAlgo = BorderedPred(), verbosity = 0)

brmf, u1 = @time PALC.continuation(
	(x, p) -> F_chan(x, p, 0.01),
	(x, p) -> (dx -> dF_chan(x, dx, p, 0.01)),
	out, a, opts_cont_mf,
	tangentAlgo = SecantPred(),
	linearAlgo = PALC.MatrixFreeBLS())

brmf, u1 = @time PALC.continuation(
	(x, p) -> F_chan(x, p, 0.01),
	(x, p) -> (dx -> dF_chan(x, dx, p, 0.01)),
	out, a, opts_cont_mf,
	tangentAlgo = BorderedPred(),
	linearAlgo = PALC.MatrixFreeBLS())
####################################################################################################
# deflation newton solver, test of jacobian expression
deflationOp = DeflationOperator(2.0, (x,y) -> dot(x,y), 1.0, [out])

# quick test of scalardM deflation
# PALC.scalardM(deflationOp, sol, 5sol)

chanDefPb   = DeflatedProblem(x -> F_chan(x,a, 0.01),x -> Jac_mat(x,a, 0.01),deflationOp)

opt_def = setproperties(opt_newton; tol = 1e-10, maxIter = 1000)
outdef1, _,_ = @time PALC.newton(
						u -> chanDefPb(u),
						out .* (1 .+0.01*rand(n)),
						opt_def)

# we now compare the jacobians for the deflated problem either using finite differences or the explicit jacobian
rhs = rand(n)
J_def_fd = PALC.finiteDifferences(u->chanDefPb(u),1.5*out)
res_fd =  J_def_fd \ rhs

Jacdf = (u0, pb::DeflatedProblem,ls = opt_def.linsolve ) -> (return (u0, pb, ls))
Jacdfsolver = DeflatedLinearSolver()

res_explicit = Jacdfsolver(Jacdf(1.5out, chanDefPb, opt_def.linsolver),rhs)[1]

println("--> Test jacobian expression for deflated problem")
@test norm(res_fd - res_explicit,Inf64) < 1e-4

opt_def = setproperties(opt_newton; tol = 1e-10, maxIter = 1000)
outdef1, _, _ = @time PALC.newton(
		x -> F_chan(x, a, 0.01),
		x -> Jac_mat(x, a, 0.01),
		out.*(1 .+ 0.01*rand(n)),
		opt_def, deflationOp)
####################################################################################################
# Fold continuation, test of Jacobian expression
indfold = 1
foldpt = FoldPoint(br, indfold)
foldpb = FoldProblemMinimallyAugmented(
		(x, α) ->   F_chan(x, α, 0.01),
		(x, α) -> (Jac_mat(x, α, 0.01)),
		(x, α) -> transpose(Jac_mat(x, α, 0.01)),
		br.foldpoint[indfold][6],
		br.foldpoint[indfold][6],
		opts_br0.newtonOptions.linsolver)
foldpb(foldpt) |> norm

outfold, _ = PALC.newtonFold(
		(x, p) -> F_chan(x, p, 0.01),
		(x, p) -> Jac_mat(x, p, 0.01),
		foldpt,
		br.foldpoint[indfold][6],
		NewtonPar(verbose=true) )
	println("--> Fold found at α = ", outfold.p, " from ", br.foldpoint[indfold][3])

# example with KrylovKit
P = Jac_mat(sol.*0,0,0)
optils = NewtonPar(verbose=true, linsolver = GMRESKrylovKit(atol=1e-9, Pl=lu(P)), tol=1e-7)
outfold, _ = PALC.newtonFold(
		(x, p) -> F_chan(x, p, 0.01),
		(x, p) -> Jac_mat(x, p, 0.01),
		foldpt,
		br.foldpoint[indfold][6],
		optils )
	println("--> Fold found at α = ", outfold.p, " from ", br.foldpoint[indfold][3])

# continuation of the fold point
optcontfold = PALC.ContinuationPar(dsmin = 0.001, dsmax = 0.15, ds= 0.01, pMax = 4.1, pMin = 0., a = 2., theta = 0.3, newtonOptions = NewtonPar(verbose=true, tol = 1e-8), maxSteps = 5)
	outfoldco, hist, flag = @time PALC.continuationFold(
						(x, α, β) ->  F_chan(x, α, β),
						(x, α, β) -> Jac_mat(x, α, β),
						br, indfold,
						0.01,
						optcontfold, plot = false)

# user defined Fold Problem
indfold = 1

# we define the following wrappers to be able to use finite differences
Bd2Vec(x) = vcat(x.u, x.p)
Vec2Bd(x) = BorderedArray(x[1:end-1], x[end])
foldpbVec(x) = Bd2Vec(foldpb(Vec2Bd(x)))

outfold, _ = PALC.newton(x -> foldpbVec(x),
			# (x, α) -> Jac_mat(x, α, 0.01),
			Bd2Vec(foldpt),
			NewtonPar(verbose=true) )
	println("--> Fold found at α = ", outfold[end], " from ", br.foldpoint[indfold][3])

rhs = rand(n+1)
Jac_fold_fdMA(u0) = PALC.finiteDifferences( u-> foldpbVec(u), u0)
J_fold_fd = Jac_fold_fdMA(Bd2Vec(foldpt))
res_fd =  J_fold_fd \ rhs

Jac_fold_MA(u0, β, pb::FoldProblemMinimallyAugmented) = (return (u0, pb, x->x))
jacFoldSolver = FoldLinearSolveMinAug()
res_explicit = jacFoldSolver(Jac_fold_MA(foldpt, 0.01, foldpb), Vec2Bd(rhs), true)

Jac_fold_MA(foldpt, 0.01, foldpb)[2]

# test whether the Jacobian Matrix for the Fold problem is correct
println("--> FD σp = ", J_fold_fd[end,end])
println("--> MA σp = ", res_explicit[end][end,end])

@test J_fold_fd[end,:] -  J_fold_fd[end,:] |> x->norm(x,Inf64) <1e-7
@test J_fold_fd[:,end] -  J_fold_fd[:,end] |> x->norm(x,Inf64) <1e-7
# this part is where we compare the FD Jacobian with Jac_chan, not really good
@test (res_explicit[end] - J_fold_fd)[1:end-1,1:end-1] |> x->norm(x,Inf64) < 1e-1

# check our solution of the bordered problem
res_exp = res_explicit[end] \ rhs
@test norm(res_exp - Bd2Vec(res_explicit[1]), Inf64) < 1e-8
####################################################################################################
# Use of different eigensolvers
opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = false, eigsolver = EigKrylovKit())
opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, pMax = 4.1, maxSteps = 250, newtonOptions = opt_newton, detectFold = true, detectBifurcation = 1, nev = 15)

br, u1 = @time PALC.continuation(
			(x,p) -> F_chan(x,p, 0.01),
			(x,p) -> (Jac_mat(x,p, 0.01)),
			printSolution = (x,p)->norm(x,Inf64),
			out,a,opts_br0,plot = false, verbosity = 0)

opts_br0 = ContinuationPar(dsmin = 0.01, dsmax = 0.15, ds= 0.01, pMax = 4.1, maxSteps = 250, newtonOptions = NewtonPar(tol =1e-8), detectFold = true, detectBifurcation = 1, nev = 15)

br, u1 = @time PALC.continuation(
			(x,p) -> F_chan(x,p, 0.01),
			(x,p) -> (Jac_mat(x,p, 0.01)),
			printSolution = (x,p)->norm(x,Inf64),
			out,a,opts_br0,plot = false, verbosity = 0)


text:
import Base: +, -, *, /
import Base.LinAlg.norm

typealias Point{N,T} NTuple{N,T}

typealias Point1{T} Tuple{T}
typealias Point2{T} Tuple{T,T}
typealias Point3{T} Tuple{T,T,T}
typealias Point4{T} Tuple{T,T,T,T}

typealias Point1f Point1{Float64}
typealias Point2f Point2{Float64}
typealias Point3f Point3{Float64}
typealias Point4f Point4{Float64}
typealias Point1i Point1{Int}
typealias Point2i Point2{Int}
typealias Point3i Point3{Int}
typealias Point4i Point4{Int}

# Point{T}(x::T, y::T, z::T) = (x,y,z)

precompile(Point, (Float64,Float64,Float64))

for op in (:+, :-, :*, :/, :.+, :.-, :.*, :./)
    @eval begin
        Base.$op{T}(x::Point2{T}, y::Point2{T}) = $op(x[1], y[1]), $op(x[2], y[2])
        Base.$op{T}(x::Point3{T}, y::Point3{T}) = $op(x[1], y[1]), $op(x[2], y[2]), $op(x[3], y[3])
        Base.$op{T}(x::Point2{T}, a::Number) = $op(x[1], a), $op(x[2], a)
        Base.$op{T}(x::Point3{T}, a::Number) = $op(x[1], a), $op(x[2], a), $op(x[3], a)
        Base.$op{T}(a::Number, x::Point2{T}) = $op(x,a)
        Base.$op{T}(a::Number, x::Point3{T}) = $op(x,a)
        precompile($op, (Point2f, Point2f))
        precompile($op, (Point3f, Point3f))
        precompile($op, (Point2f, Float64))
        precompile($op, (Point3f, Float64))
    end
end

norm(x::Point2f) = sqrt(x[1]*x[1] + x[2]*x[2])
norm(x::Point3f) = sqrt(x[1]*x[1] + x[2]*x[2] + x[3]*x[3])

precompile(norm, (Point2f,))
precompile(norm, (Point3f,))

function Base.(:*){T}(A::Array{T,2}, x::Point3{T})
    A[1,1]*x[1] + A[1,2]*x[2] + A[1,3]*x[3], A[2,1]*x[1] + A[2,2]*x[2] + A[2,3]*x[3], A[3,1]*x[1] + A[3,2]*x[2] + A[3,3]*x[3]
end

precompile(*, (Array{Float64,2}, Point3f))

@doc """
`dist(u, v)` computes the Euclidean distance betweens points `u` and `v`.
""" ->
function dist{T}(u::Point3{T}, v::Point3{T})
    dx = u[1]-v[1]
    dy = u[2]-v[2]
    dz = u[3]-v[3]
    sqrt(dx*dx + dy*dy + dz*dz)
end

dist{T}(u::AbstractVector{Point3{T}}, v::Point3{T}) = T[dist(u[i],v) for i in 1:length(u)]
dist{T}(u::Point3{T}, v::AbstractVector{Point3{T}}) = T[dist(u,v[i]) for i in 1:length(v)]

function dist!{T}(dst::AbstractMatrix{T}, u::AbstractVector{Point3{T}}, v::AbstractVector{Point3{T}})
    @assert size(dst,1) == length(u) && size(dst,2) == length(v)
    for i in 1:length(v)
        for j in 1:length(u)
            dst[j,i] = dist(u[j],v[i])
        end
    end
end

function dist{T}(u::AbstractVector{Point3{T}}, v::AbstractVector{Point3{T}})
    dst = Array{T}(length(u),length(v))
    dist!(dst,u,v)
    dst
end

function dist!{T}(dst::AbstractMatrix{T}, u::AbstractVector{Point3{T}})
    n = length(u)
    @assert size(dst,1) == size(dst,2) == n
    for i in 1:n
        for j in 1:i
            dst[j,i] = dist(u[j],u[i])
        end
    end
end

function dist{T}(u::AbstractVector{Point3{T}})
    dst = Array{T}(length(u),length(u))
    dist!(dst,u)
    Symmetric(dst)
end

function randp(d::Int, dims...)
    dst = Array{NTuple{d,Float64}}(dims...)
    for i in eachindex(dst)
        dst[i] = ntuple(i->rand(),d)
    end
    dst
end

precompile(dist , (Point3f, Point3f))
precompile(dist , (Vector{Point3f}, Point3f))
precompile(dist , (Point3f, Vector{Point3f}))
precompile(dist!, (Matrix{Point3f}, Vector{Point3f}, Vector{Point3f}))
precompile(dist , (Vector{Point3f}, Vector{Point3f}))
precompile(dist!, (Matrix{Point3f}, Vector{Point3f}))
precompile(dist , (Vector{Point3f},))

max_stars_repo_name: johncsnyder/Element


text:
"""
This module provides the [`newconnection`](@ref) function with support for:
- Opening TCP and SSL connections.
- Reusing connections for multiple Request/Response Messages

This module defines a [`Connection`](@ref)
struct to manage the lifetime of a connection and its reuse.
Methods are provided for `eof`, `readavailable`,
`unsafe_write` and `close`.
This allows the `Connection` object to act as a proxy for the
`TCPSocket` or `SSLContext` that it wraps.

The [`POOL`](@ref) is used to manage connection pooling. Connections
are identified by their host, port, whether they require
ssl verification, and whether they are a client or server connection.
If a subsequent request matches these properties of a previous connection
and limits are respected (reuse limit, idle timeout), and it wasn't otherwise
remotely closed, a connection will be reused.
"""
module ConnectionPool

export Connection, newconnection, getrawstream, inactiveseconds

using ..IOExtras, ..Sockets

import ..@debug, ..DEBUG_LEVEL, ..taskid
import ..@require, ..precondition_error, ..@ensure, ..postcondition_error
using MbedTLS: SSLConfig, SSLContext, setup!, associate!, hostname!, handshake!
import NetworkOptions

const default_connection_limit = 8
const default_pipeline_limit = 16
const nolimit = typemax(Int)

include("connectionpools.jl")
using .ConnectionPools

"""
    Connection{T <: IO}

A `TCPSocket` or `SSLContext` connection to a HTTP `host` and `port`.

Fields:
- `host::String`
- `port::String`, exactly as specified in the URI (i.e. may be empty).
- `idle_timeout`, No. of seconds to maintain connection after last request/response.
- `peerip`, remote IP adress (used for debug/log messages).
- `peerport`, remote TCP port number (used for debug/log messages).
- `localport`, local TCP port number (used for debug messages).
- `io::T`, the `TCPSocket` or `SSLContext.
- `clientconnection::Bool`, whether the Connection was created from client code (as opposed to server code)
- `buffer::IOBuffer`, left over bytes read from the connection after
   the end of a response header (or chunksize). These bytes are usually
   part of the response body.
- `timestamp`, time data was last received.
"""
mutable struct Connection <: IO
    host::String
    port::String
    idle_timeout::Int
    require_ssl_verification::Bool
    peerip::IPAddr # for debugging/logging
    peerport::UInt16 # for debugging/logging
    localport::UInt16 # debug only
    io::IO
    clientconnection::Bool
    buffer::IOBuffer
    timestamp::Float64
    readable::Bool
    writable::Bool
end

"""
    connectionkey

Used for "hashing" a Connection object on just the key properties necessary for determining
connection re-useability. That is, when a new request calls `newconnection`, we take the
request parameters of what socket type, the host and port, and if ssl
verification is required, and if an existing Connection was already created with the exact
same parameters, we can re-use it (as long as it's not already being used, obviously).
"""
connectionkey(x::Connection) = (typeof(x.io), x.host, x.port, x.require_ssl_verification, x.clientconnection)

Connection(host::AbstractString, port::AbstractString,
           idle_timeout::Int,
           require_ssl_verification::Bool, io::IO, client=true) =
    Connection(host, port, idle_timeout,
                require_ssl_verification,
                safe_getpeername(io)..., localport(io),
                io, client, PipeBuffer(), time(), false, false)

Connection(io; require_ssl_verification::Bool=true) =
    Connection("", "", 0, require_ssl_verification, io, false)

getrawstream(c::Connection) = c.io

inactiveseconds(c::Connection)::Float64 = time() - c.timestamp

Base.unsafe_write(c::Connection, p::Ptr{UInt8}, n::UInt) =
    unsafe_write(c.io, p, n)

Base.isopen(c::Connection) = isopen(c.io)

"""
    flush(c::Connection)

Flush a TCP buffer by toggling the Nagle algorithm off and on again for a socket.
This forces the socket to send whatever data is within its buffer immediately,
rather than waiting 10's of milliseconds for the buffer to fill more.
"""
function Base.flush(c::Connection)
    # Flushing the TCP buffer requires support for `Sockets.nagle()`
    # which was only added in Julia v1.3
    @static if VERSION >= v"1.3"
        sock = tcpsocket(c.io)
        # I don't understand why uninitializd sockets can get here, but they can
        if sock.status ∉ (Base.StatusInit, Base.StatusUninit) && isopen(sock)
            Sockets.nagle(sock, false)
            Sockets.nagle(sock, true)
        end
    end
end

Base.isreadable(c::Connection) = c.readable
Base.iswritable(c::Connection) = c.writable

function Base.eof(c::Connection)
    @require isreadable(c) || !isopen(c)
    if bytesavailable(c) > 0
        return false
    end
    try
        return eof(c.io)
    catch IOError
        # suppose a reset connection, assume eof is true
        return true
    end
end

Base.bytesavailable(c::Connection) = bytesavailable(c.buffer) +
                                     bytesavailable(c.io)

function Base.read(c::Connection, nb::Int)
    nb = min(nb, bytesavailable(c))
    bytes = Base.StringVector(nb)
    GC.@preserve bytes unsafe_read(c, pointer(bytes), nb)
    return bytes
end

function Base.read(c::Connection, ::Type{UInt8})
    if bytesavailable(c.buffer) == 0
        read_to_buffer(c)
    end
    return read(c.buffer, UInt8)
end

function Base.unsafe_read(c::Connection, p::Ptr{UInt8}, n::UInt)
    l = bytesavailable(c.buffer)
    if l > 0
        nb = min(l, n)
        unsafe_read(c.buffer, p, nb)
        p += nb
        n -= nb
        c.timestamp = time()
    end
    if n > 0
        # try-catch underlying errors here
        # as the Connection object, we don't really care
        # if the underlying socket was closed/terminated
        # or just plain reached EOF, so we wrap any
        # Base.IOErrors and just throw as EOFError
        # that way we get more consistent errors thrown
        # at the headers/body parsing level
        try
            unsafe_read(c.io, p, n)
            c.timestamp = time()
        catch e
            e isa Base.IOError && throw(EOFError())
            rethrow(e)
        end
    end
    return nothing
end

function read_to_buffer(c::Connection, sizehint=4096)
    buf = c.buffer

    # Reset the buffer if it is empty.
    if bytesavailable(buf) == 0
        buf.size = 0
        buf.ptr = 1
    end

    # Wait for data.
    if eof(c.io)
        throw(EOFError())
    end

    # Read from stream into buffer.
    n = min(sizehint, bytesavailable(c.io))
    buf = c.buffer
    Base.ensureroom(buf, n)
    GC.@preserve buf unsafe_read(c.io, pointer(buf.data, buf.size + 1), n)
    buf.size += n
end

"""
Read until `find_delimiter(bytes)` returns non-zero.
Return view of bytes up to the delimiter.
"""
function Base.readuntil(c::Connection, f::Function #=Vector{UInt8} -> Int=#,
                                        sizehint=4096)::ByteView
    buf = c.buffer
    if bytesavailable(buf) == 0
        read_to_buffer(c, sizehint)
    end
    while (bytes = readuntil(buf, f)) === nobytes
        read_to_buffer(c, sizehint)
    end
    return bytes
end

"""
    startwrite(::Connection)
"""
function IOExtras.startwrite(c::Connection)
    @require !iswritable(c)
    c.writable = true
    @debug 2 "👁  Start write:$c"
    return
end

"""
    closewrite(::Connection)

Signal that an entire Request Message has been written to the `Connection`.
"""
function IOExtras.closewrite(c::Connection)
    @require iswritable(c)
    c.writable = false
    @debug 2 "🗣  Write done: $c"
    flush(c)
    return
end

"""
    startread(::Connection)
"""
function IOExtras.startread(c::Connection)
    @require !isreadable(c)
    c.timestamp = time()
    c.readable = true
    @debug 2 "👁  Start read: $c"
    return
end

"""
Wait for `c` to receive data or reach EOF.
Close `c` on EOF or if response data arrives when no request was sent.
"""
function monitor_idle_connection(c::Connection)
    if eof(c.io)                                  ;@debug 2 "💀  Closed:     $c"
        close(c.io)
    end
end

"""
    closeread(::Connection)

Signal that an entire Response Message has been read from the `Connection`.
"""
function IOExtras.closeread(c::Connection)
    @require isreadable(c)
    c.readable = false
    @debug 2 "✉  Read done:  $c"
    if c.clientconnection
        release(POOL, connectionkey(c), c)
        # Ignore SSLContext as it already monitors idle connections for TLS close_notify messages
        !(c.io isa SSLContext) && @async monitor_idle_connection(c)
    end
    return
end

Base.wait_close(c::Connection) = Base.wait_close(tcpsocket(c.io))

function Base.close(c::Connection)
    if iswritable(c)
        closewrite(c)
    end
    if isreadable(c)
        closeread(c)
    end
    try
        close(c.io)
        if bytesavailable(c) > 0
            purge(c)
        end
    catch
        # ignore errors closing underlying socket
    end
    return
end

"""
    purge(::Connection)

Remove unread data from a `Connection`.
"""
function purge(c::Connection)
    @require !isopen(c.io)
    while !eof(c.io)
        readavailable(c.io)
    end
    c.buffer.size = 0
    c.buffer.ptr = 1
    @ensure bytesavailable(c) == 0
end

"""
    closeall()

Close all connections in `pool`.
"""
function closeall()
    ConnectionPools.reset!(POOL)
    return
end

"""
    POOL

Global connection pool keeping track of active connections.
"""
const POOL = Pool(Connection)

"""
    getconnection(type, host, port) -> Connection

Find a reusable `Connection` in the `pool`,
or create a new `Connection` if required.
"""
function newconnection(::Type{T},
                       host::AbstractString,
                       port::AbstractString;
                       connection_limit=default_connection_limit,
                       pipeline_limit=default_pipeline_limit,
                       idle_timeout=typemax(Int),
                       reuse_limit=nolimit,
                       require_ssl_verification::Bool=NetworkOptions.verify_host(host, "SSL"),
                       kw...)::Connection where {T <: IO}
    return acquire(
            POOL,
            (T, host, port, require_ssl_verification, true);
            max_concurrent_connections=Int(connection_limit),
            idle_timeout=Int(idle_timeout)) do
        Connection(host, port,
            idle_timeout, require_ssl_verification,
            getconnection(T, host, port;
                require_ssl_verification=require_ssl_verification, kw...)
        )
    end
end

function keepalive!(tcp)
    @debug 2 "setting keepalive on tcp socket"
    err = ccall(:uv_tcp_keepalive, Cint, (Ptr{Nothing}, Cint, Cuint),
                                          tcp.handle, 1, 1)
    err != 0 && error("error setting keepalive on socket")
    return
end

struct ConnectTimeout <: Exception
    host
    port
end

function getconnection(::Type{TCPSocket},
                       host::AbstractString,
                       port::AbstractString;
                       keepalive::Bool=false,
                       connect_timeout::Int=0,
                       readtimeout::Int=0,
                       kw...)::TCPSocket

    p::UInt = isempty(port) ? UInt(80) : parse(UInt, port)

    @debug 2 "TCP connect: $host:$p..."

    addrs = Sockets.getalladdrinfo(host)

    connect_timeout = connect_timeout == 0 && readtimeout > 0 ? readtimeout : connect_timeout

    lasterr = ErrorException("unknown connection error")

    for addr in addrs
        if connect_timeout == 0
            try
                tcp = Sockets.connect(addr, p)
                keepalive && keepalive!(tcp)
                return tcp
            catch err
                lasterr = err
                continue # to next ip addr
            end
        else
            tcp = Sockets.TCPSocket()
            Sockets.connect!(tcp, addr, p)

            timeout = Ref{Bool}(false)
            @async begin
                sleep(connect_timeout)
                if tcp.status == Base.StatusConnecting
                    timeout[] = true
                    tcp.status = Base.StatusClosing
                    ccall(:jl_forceclose_uv, Nothing, (Ptr{Nothing},), tcp.handle)
                    #close(tcp)
                end
            end
            try
                Sockets.wait_connected(tcp)
                keepalive && keepalive!(tcp)
                return tcp
            catch err
                if timeout[]
                    lasterr = ConnectTimeout(host, port)
                else
                    lasterr = err
                end
                continue # to next ip addr
            end
        end
    end
    # If no connetion could be set up, to any address, throw last error
    throw(lasterr)
end

const nosslconfig = SSLConfig()
default_sslconfig = nothing
noverify_sslconfig = nothing

function global_sslconfig(require_ssl_verification::Bool)::SSLConfig
    global default_sslconfig
    global noverify_sslconfig
    if default_sslconfig === nothing
        default_sslconfig = SSLConfig(true)
        noverify_sslconfig = SSLConfig(false)
    end
    return require_ssl_verification ? default_sslconfig : noverify_sslconfig
end

function getconnection(::Type{SSLContext},
                       host::AbstractString,
                       port::AbstractString;
                       kw...)::SSLContext

    port = isempty(port) ? "443" : port
    @debug 2 "SSL connect: $host:$port..."
    tcp = getconnection(TCPSocket, host, port; kw...)
    return sslconnection(tcp, host; kw...)
end

function sslconnection(tcp::TCPSocket, host::AbstractString;
                       require_ssl_verification::Bool=NetworkOptions.verify_host(host, "SSL"),
                       sslconfig::SSLConfig=nosslconfig,
                       kw...)::SSLContext

    if sslconfig === nosslconfig
        sslconfig = global_sslconfig(require_ssl_verification)
    end

    io = SSLContext()
    setup!(io, sslconfig)
    associate!(io, tcp)
    hostname!(io, host)
    handshake!(io)
    return io
end

function sslupgrade(c::Connection,
                    host::AbstractString;
                    require_ssl_verification::Bool=NetworkOptions.verify_host(host, "SSL"),
                    kw...)::Connection
    # first we release the original connection, but we don't want it to be
    # reused in the pool, because we're hijacking the TCPSocket
    release(POOL, connectionkey(c), c; return_for_reuse=false)
    # now we hijack the TCPSocket and upgrade to SSLContext
    tls = sslconnection(c.io, host;
                        require_ssl_verification=require_ssl_verification,
                        kw...)
    conn = Connection(host, "", 0, require_ssl_verification, tls)
    return acquire(POOL, connectionkey(conn), conn)
end

function Base.show(io::IO, c::Connection)
    nwaiting = applicable(tcpsocket, c.io) ? bytesavailable(tcpsocket(c.io)) : 0
    print(
        io,
        tcpstatus(c), " ",
        "$(lpad(round(Int, time() - c.timestamp), 3))s ",
        c.host, ":",
        c.port != "" ? c.port : Int(c.peerport), ":", Int(c.localport),
        bytesavailable(c.buffer) > 0 ?
            " $(bytesavailable(c.buffer))-byte excess" : "",
        nwaiting > 0 ? " $nwaiting bytes waiting" : "",
        DEBUG_LEVEL[] > 1 && applicable(tcpsocket, c.io) ?
            " $(Base._fd(tcpsocket(c.io)))" : "")
end

function tcpstatus(c::Connection)
    if !applicable(tcpsocket, c.io)
        return ""
    end
    s = Base.uv_status_string(tcpsocket(c.io))
        if s == "connecting" return "🔜🔗"
    elseif s == "open"       return "🔗 "
    elseif s == "active"     return "🔁 "
    elseif s == "paused"     return "⏸ "
    elseif s == "closing"    return "🔜💀"
    elseif s == "closed"     return "💀 "
    else
        return s
    end
end

end # module ConnectionPool


text:
using Base.Test
using Base.Test: AbstractTestSet, DefaultTestSet
using LightGraphs

import Base.Test: record, finish

export @testprop, QuickTestSet

const ntests_default = 100::Int
const maxntries_default_factor = 4::Int

abstract type Element{P} end
abstract type ElementVector{P} end
abstract type Expression{P} end

abstract type Elem{P} end
abstract type ElemV{P} end
abstract type Exp{P} end

type QuickTestSet <: AbstractTestSet
  testset::DefaultTestSet
  function QuickTestSet(desc; var_vals=Tuple{String,Symbol}[])
    for (s,v) in var_vals
      desc = string(desc, "  $s = $v")
    end
    new(DefaultTestSet(desc))
  end
end

record(ts::QuickTestSet,r) = record(ts.testset,r)
finish(ts::QuickTestSet) = finish(ts.testset)

function testprop_(ntests_::Union{Int,Void}, prop::Expr, type_asserts::Vector{Expr})
  if isa(ntests_, Void)
    ntests = ntests_default
  else
    ntests = ntests_
  end

  if isempty(type_asserts)
    description = string(prop)
  else
    description = string(prop, " | ", [string(ta, " ") for ta in type_asserts]...)
  end

  types = Dict{Symbol,Expr}()
  extract_types_from_asserts!(types, type_asserts)
  if isempty(types)
    extract_types_from_expression!(types, prop)
  else
    check_types_from_property(types, prop)
  end
  types = topological_sort_types(types)

  Expr(:for, Expr(:(=), :testsize, Expr(:(:), 1, ntests)),
       Expr(:block,
            init_test(:testsize, types)...,
            Expr(:macrocall, Symbol("@testset"),
                 esc(:QuickTestSet),
                 Expr(:(=), esc(:var_vals), esc(var_vals(types))),
                 string(description, ": case"),
                 Expr(:block, Expr(:macrocall, Symbol("@test"), esc(prop)))
                )
           )
      )
end

function testcondprop_(ntests::Union{Int,Void}, maxntries::Union{Int,Void}, prop::Expr, cond_and_type_asserts)
  cond_and_type_asserts = Vector{Expr}(collect(cond_and_type_asserts))

  if !isempty(cond_and_type_asserts) && cond_and_type_asserts[1].head != :(::)
    cond = cond_and_type_asserts[1]
    type_asserts = cond_and_type_asserts[2:length(cond_and_type_asserts)]
    return testcondprop_(ntests, nothing, prop, cond, type_asserts)
  elseif isa(maxntries, Void)
    type_asserts = cond_and_type_asserts
    return testprop_(ntests, prop, type_asserts)
  else
    error("maximum number of tries given, but no test condition present")
  end
end

function testcondprop_(ntests_::Union{Int,Void}, maxntries_::Union{Int,Void}, prop::Expr, cond::Expr, type_asserts::Vector{Expr})
  if isa(ntests_, Void)
    ntests = ntests_default
  else
    ntests = ntests_
  end

  if isa(maxntries_, Void)
    maxntries = maxntries_default_factor * ntests
  else
    maxntries = maxntries_
  end

  if isempty(type_asserts)
    description = string(prop, " | ", cond)
  else
    description = string(prop, " | ", cond, " ; ",
                         [string(ta, " ") for ta in type_asserts]...)
  end

  types = Dict{Symbol,Expr}()
  extract_types_from_expression!(types, cond)
  extract_types_from_asserts!(types, type_asserts)
  if isempty(types)
    extract_types_from_expression!(types, prop)
  else
    check_types_from_property(types, prop)
  end
  types = topological_sort_types(types)

  Expr(:block,
       Expr(:(=), :testsize, 1),
       Expr(:(=), :curtest, 1),
       Expr(:for, Expr(:(=), :xtry, Expr(:(:), 1, maxntries)),
            Expr(:block,
                 Expr(:if, Expr(:call, :(>), :xtry, Expr(:call, :(*), :maxntries_default_factor, :testsize)),
                      Expr(:(+=), :testsize, 1),
                     ),
                 init_test(:testsize, types)...,
                 Expr(:(||), esc(cond), :(continue)),
                 Expr(:macrocall, Symbol("@testset"),
                      esc(:QuickTestSet),
                      Expr(:(=), esc(:var_vals), esc(var_vals(types))),
                      string(description, ": case"),
                      Expr(:block, Expr(:macrocall, Symbol("@test"), esc(prop)))
                     ),
                 Expr(:(+=), :testsize, 1),
                 Expr(:(+=), :curtest, 1),
                 Expr(:if, Expr(:call, :(>), :curtest, ntests),
                      Expr(:break)
                )
           )
       ),
       Expr(:if, Expr(:call, :(<=), :curtest, ntests),
            Expr(:call, :warn, string("Did not exhaust number of test cases for ", description))
           )
      )
end

function topological_sort_types(types)
  g = DiGraph()
  symbol_labels = Dict{Symbol,Int}()
  symbol_list = Symbol[]
  for (s,t) in types
    add_vertex!(g)
    symbol_labels[s] = nv(g)
    push!(symbol_list,s)
  end
  for (s,t) in types
    if t.head == :curly && t.args[1] in [:Element, :Elem, :ElementVector, :ElemV]
      if haskey(symbol_labels, t.args[2])
        add_edge!(g, symbol_labels[t.args[2]], symbol_labels[s])
      end
    elseif t.head == :call && t.args[1].head == :curly && t.args[1].args[1] in [:Expression, :Exp]
      for ts in t.args[1].args[2:end]
        if haskey(symbol_labels, ts)
          add_edge!(g, symbol_labels[ts], symbol_labels[s])
        end
      end
    end
  end

  sorted_types = Tuple{Symbol,Expr}[]
  for sx in topological_sort_by_dfs(g)
    s = symbol_list[sx]
    push!(sorted_types, (s, types[s]))
  end
  return sorted_types
end

function init_test(size_symbol, types)
  exprs = Expr[]
  for (s,t) in types
    if t.head == :curly && t.args[1] in [:Element, :Elem]
      push!(exprs, Expr(:(=), esc(s), Expr(:call, generate_test_value, esc(t.args[2]), size_symbol)))
    elseif t.head == :curly && t.args[1] in [:ElementVector, :ElemV]
      push!(exprs, Expr(:(=), esc(s), Expr(:call, generate_test_vector, esc(t.args[2]), size_symbol)))
    elseif t.head == :call && t.args[1].head == :curly && t.args[1].args[1] in [:Expression, :Exp]
      push!(exprs, Expr(:(=), esc(s), esc(t.args[2])))
    else
      push!(exprs, Expr(:(=), esc(s), Expr(:call, generate_test_value, esc(t), size_symbol)))
    end
  end
  return exprs
end

function var_vals(types)
  return Expr(:vect, [Expr(:tuple, string(s), esc(s)) for (s,t) in types]...)
end

macro testprop(prop::Expr, cond_and_type_asserts...)
  return testcondprop_(nothing, nothing, prop, cond_and_type_asserts)
end

macro testprop(ntests::Int, prop::Expr, cond_and_type_asserts...)
  return testcondprop_(ntests, nothing, prop, cond_and_type_asserts)
end

macro testprop(ntests::Int, maxntries::Int, prop::Expr, cond_and_type_asserts...)
  return testcondprop_(ntests, maxntries, prop, cond_and_type_asserts)
end

function extract_types_from_expression!(types::Dict{Symbol,Expr}, e::Expr)
  stack = [e]
  while !isempty(stack)
    e = pop!(stack)
    if isa(e.args, Array)
      append!(stack, filter((f) -> isa(f,Expr), e.args))
    end

    if e.head == :(::)
      add_type_to_dict!(types, e)
    end
  end
end

function extract_types_from_asserts!(types::Dict{Symbol,Expr}, es::Array{Expr,1})
  for e in es
    if e.head != :(::)
      throw( ArgumentError(string("Type assertions must be of the form ",
                                  "symbol::type; ",
                                  "found instread: ", string(e))) )
    end

    add_type_to_dict!(types, e)
  end
end

function check_types_from_property(types::Dict{Symbol,Expr}, prop::Expr)
  prop_types = Dict{Symbol,Expr}()
  extract_types_from_expression!(prop_types, prop)
  for (s,t) in prop_types
    if haskey(types, s) && types[s] != t
      throw( ArgumentError(string("Test asserts that $s is of type $t, ",
                                  "but property requires $(types[s])")) )
    end
  end
end

function add_type_to_dict!(types::Dict{Symbol,Expr}, e::Expr)
  if isa(e.args[2], Symbol)
    t = Expr(:block, e.args[2])
  else
    t = e.args[2]
  end
  add_type_to_dict!(types, e.args[1], t)
end

function add_type_to_dict!(types::Dict{Symbol,Expr}, s::Symbol, t::Expr)
  if haskey(types, s)
    if types[s] != t
      throw( ArgumentError(string("Incompatible type symbols ", types[s], " and ", t,
                                  " specified for symbol ", s)) )
    end
  else
    types[s] = t
  end
end


text:
# From equations in "THE LOXODROME ON AN ELLIPSOID"
# http://www.mygeodesy.id.au/documents/Loxodrome%20on%20Ellipsoid.pdf

"""
    function loxodrome_inverse(lon1, lat1, lon2, lat2, a=6378137.0, f=0.0033528106647474805)

Compute the inverse problem of a loxodrome on the ellipsoid.

Given latitudes and longitudes of P1 and P2 on the ellipsoid, compute the azimuth a12 of
the loxodrome P1P2, the arc length s along the loxodrome curve.

Args:

- `lon1, lat1, lon2, lat2`: - longitude and latitude of starting and end points (degrees).
- `a` - major axis of the ellipsoid (meters). Default values for WGS84
- `f` - flattening od the ellipsoid (default = 1 / 298.257223563)

### Returns
- Distance (meters) and azimuth from P1 to P2

## Example: Compute the distance and azimuth beyween points (0,0) and (5,5)

    dist, azim = loxodrome_inverse(0,0,5,5)
"""
function loxodrome_inverse(lon1, lat1, lon2, lat2, a=6378137.0, f=0.0033528106647474805)
	D2R = pi / 180
	lon1 *= D2R;	lat1 *= D2R
	lon2 *= D2R;	lat2 *= D2R
	e2 = f * (2 - f)

	isolat1 = isometric_lat(lat1, e2)
	isolat2 = isometric_lat(lat2, e2)

	# Compute changes in isometric latitude and longitude between P1 and P2
	Az12 = atan((lon2 - lon1), (isolat2 - isolat1))		# The azimuth
	# Compute distance along loxodromic curve
	m1 = meridian_dist(lat1, a, e2)
	m2 = meridian_dist(lat2, a, e2)
	lox_s = (m2 - m1) / cos(Az12);
	return lox_s, Az12 / D2R
end

# -------------------------------------------------------------------------------------------------------
"""
    loxodrome_direct(lon, lat, azimuth, distance, a=6378137.0, f=0.0033528106647474805)

Compute the direct problem of a loxodrome on the ellipsoid.

Given latitude and longitude of P1, azimuth a12 of the loxodrome P1P2 and the arc length s along the
loxodrome curve, compute the latitude and longitude of P2.

Args:

- `lon, lat`: - longitude, latitude (degrees) of starting point.
- `azimuth`:  - azimuth (degrees)
- `distance`: - distance to move from (lat,lon) in meters
- `a` - major axis of the ellipsoid (meters). Default values for WGS84
- `f` - flattening od the ellipsoid (default = 1 / 298.257223563)

### Returns
- [lon lat] of destination after moving for [distance] metres in [azimuth] direction.

## Example: Compute the end point at a bearing of 45 degrees 10000 meters from point 0,0

    loxo = loxodrome_direct(0,0,45, 10000)
"""
function loxodrome_direct(lon, lat, azim, dist, a=6378137.0, f=0.0033528106647474805)
	(rem(azim, 90) == 0) && (azim -= 0.0001)	# Equator is a degenerated case. Comparing with GeographicLib sows 4 dec is ~good
	D2R = pi / 180
	lon *= D2R;		lat *= D2R;		azim *= D2R
	e2 = f * (2 - f)
	m1 = meridian_dist(lat, a, e2)
	m2 = dist * cos(azim) + m1
	dm = m2 - m1

	A0, A2, A4, A6, A8, A10 = f_phi_coeff(e2)
	ae2 = a * (1 - e2)
	phi_n, dphi = lat, 1e3
	n = 0
	while (dphi > 1e-10 && n < 6)
		phi_n1 = phi_n - (ae2 * f_phi(phi_n, A0, A2, A4, A6, A8, A10) - dm) / (ae2 * f_prime_phi(phi_n, A0, A2, A4, A6, A8, A10))
		dphi = abs(phi_n1 - phi_n)
		phi_n = phi_n1
		n += 1
	end

	pi2 = pi/2
	if (-pi2 <= phi_n <= pi2)		# For long dist we may have out off bounds lats
		isolat1 = isometric_lat(lat, e2)
		isolat2 = isometric_lat(phi_n, e2)
		dlon = (isolat2 - isolat1) * tan(azim)
		lon2 = lon + dlon
		return [lon2/D2R phi_n/D2R]
	else
		return [NaN NaN]
	end
end

function f_phi(lat, A0, A2, A4, A6, A8, A10)
	A0*lat - (A2/2)*sin(2*lat) + (A4/4)*sin(4*lat) - (A6/6)*sin(6*lat) + (A8/8)*sin(8*lat) - (A10/10)*sin(10*lat)
end

function f_prime_phi(lat, A0, A2, A4, A6, A8, A10)
	A0 - A2 * cos(2*lat) + A4 * cos(4*lat) - A6 * cos(6*lat) + A8 * cos(8*lat) - A10 * cos(10*lat)
end

function meridian_dist(lat, a, e2)
	A0, A2, A4, A6, A8, A10 = f_phi_coeff(e2)
	a*(1-e2) * f_phi(lat, A0, A2, A4, A6, A8, A10)
end

function f_phi_coeff(e2)
	e4 = e2*e2;			# powers of eccentricity
	e6 = e4*e2;
	e8 = e6*e2;
	e10 = e8*e2;
	# coefficients of series expansion for meridian distance
	A0 = 1 + (3/4)*e2 + (45/64)*e4 + (175/256)*e6 + (11025/16384)*e8 + (43659/65536)*e10
	A2 =     (3/4)*e2 + (15/16)*e4 + (525/512)*e6 + (2205/2048)*e8   + (72765/65536)*e10
	A4 =                (15/64)*e4 + (105/256)*e6 + (2205/4096)*e8   + (10395/16384)*e10
	A6 =                             (35/512)*e6  + (315/2048)*e8    + (31185/131072)*e10
	A8 =                                            (315/16384)*e8   + (3465/65536)*e10
	A10 =                                                              (693/131072)*e10
	return A0, A2, A4, A6, A8, A10
end

function isometric_lat(lat, ecc2)
	# Compute isometric latitude. ECC2 is the squared eccentricity = f*(2-f), where f is the flattening
	ecc = sqrt(ecc2)
	x = ecc * sin(lat)
	y = (1-x) / (1+x)
	z = pi/4 + lat/2
	log(tan(z) * (y^(ecc/2)))
end


text:
using Documenter
using CorrelationFunctions
using CorrelationTrackers
using MaterialReconstruction

makedocs(sitename = "MaterialReconstruction.jl documentation",
         format = Documenter.HTML(prettyurls = false))

deploydocs(repo = "github.com/shamazmazum/MaterialReconstruction.jl.git")


text:
module Visualization
    using Plots
    using Measures

    function plot_synth_QA(ImTrue, ImInit, ImSynth; fname=nothing)

        # -------- define plot1 to append plot to a list
        function plot1(ps, image; clim=nothing, bin=1.0, fsz=16, label=nothing)
            # ps    - list of plots to append to
            # image - image to heatmap
            # clim  - tuple of color limits (min, max)
            # bin   - bin down by this factor, centered on nx/2+1
            # fsz   - font size
            marg   = 1mm
            nx, ny = size(image)
            nxb    = nx/round(Integer, 2*bin)

            # -------- center on nx/2+1
            i0 = max(1,round(Integer, (nx/2+2)-nxb-1))
            i1 = min(nx,round(Integer, (nx/2)+nxb+1))
            lims = [i0,i1]
            subim = image[i0:i1,i0:i1]
            push!(ps, heatmap(image, aspect_ratio=:equal, clim=clim,
                xlims=lims, ylims=lims, size=(400,400),
                legend=false, xtickfontsize=fsz, ytickfontsize=fsz,#tick_direction=:out,
                rightmargin=marg, leftmargin=marg, topmargin=marg, bottommargin=marg))
            if label != nothing
                annotate!(lims'*[.96,.04],lims'*[.09,.91],text(label,:left,:white,32))
            end
            return
        end

        # -------- initialize array of plots
        ps   = []
        clim  = (minimum(ImTrue), maximum(ImTrue))
        resit = ImInit-ImTrue
        resir = ImInit-ImSynth
        resrt = ImSynth-ImTrue
        clim2 = (minimum([minimum(resit), minimum(resir), minimum(resrt)]), maximum([maximum(resit), maximum(resir), maximum(resrt)]))
        println(clim, clim2)

        # -------- 6 panel QA plot

        plot1(ps, ImTrue, clim=clim, label="True")
        plot1(ps, ImSynth, clim=clim, label="Synth")
        plot1(ps, ImInit, clim=clim, label="Init")
        plot1(ps, resit, clim=clim2, label="Init-True")
        plot1(ps, resir, clim=clim2, label="Init-Synth")
        plot1(ps, resrt, clim=clim2, label="Synth-True")
        #=
        heatmap(ImTrue, clim=clim, label="True")
        heatmap(ImSynth, clim=clim, label="Synth")
        heatmap(ImInit, clim=clim, label="Init")
        heatmap(resit, clim=clim2, label="Init-True")
        heatmap(resir, clim=clim2, label="Init-Synth")
        heatmap(resrt, clim=clim2, label="Synth-True")
        =#
        myplot = plot(ps..., layout=(3,2), size=(1400,2000))
        if fname!=nothing
            savefig(myplot, fname)
        end
    end

    function plot_diffscales(images, titles; fname=nothing)
        pl12 = plot(
            heatmap(images[1], title=titles[1]),
            heatmap(images[2], title=titles[2]),
            heatmap(images[3],title= titles[3]),
            heatmap(images[4], title=titles[4]);
            layout=4,
        )
        plot!(size=(1000, 800))
        if fname!=nothing
            savefig(pl12, fname)
        end
    end
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
test = """
    sesenwnenenewseeswwswswwnenewsewsw
    neeenesenwnwwswnenewnwwsewnenwseswesw
    seswneswswsenwwnwse
    nwnwneseeswswnenewneswwnewseswneseene
    swweswneswnenwsewnwneneseenw
    eesenwseswswnenwswnwnwsewwnwsene
    sewnenenenesenwsewnenwwwse
    wenwwweseeeweswwwnwwe
    wsweesenenewnwwnwsenewsenwwsesesenwne
    neeswseenwwswnwswswnw
    nenwswwsewswnenenewsenwsenwnesesenew
    enewnwewneswsewnwswenweswnenwsenwsw
    sweneswneswneneenwnewenewwneswswnese
    swwesenesewenwneswnwwneseswwne
    enesenwswwswneneswsenwnewswseenwsese
    wnwnesenesenenwwnenwsewesewsesesew
    nenewswnwewswnenesenwnesewesw
    eneswnwswnwsenenwnwnwwseeswneewsenese
    neswnwewnwnwseenwseesewsenwsweewe
    wseweeenwnesenwwwswnew
    """

function parse_input(input::AbstractString)
    lines = split(input, "\n"; keepempty = false)
    return lines
end

function parse_path(input::AbstractString)::Vector{Symbol}
    isempty(input) && return []
    l = input[1] ∈ ['n', 's'] ? 2 : 1
    return [Symbol(input[1:l]); parse_path(input[l + 1:end])]
end

function move(path::AbstractString;
              init::Tuple = (0, 0))::Tuple
    return move(parse_path(path); init)
end

function move(path::AbstractVector{<: Symbol};
               init::Tuple = (0, 0))::Tuple
    return reduce(path; init) do pos, dir
        return move(dir, pos)
    end
end

move(dir::Symbol, i::Tuple) = move(Val(dir), i)
move(::Val{:w},  i::Tuple) = i .- (0, 1)
move(::Val{:e},  i::Tuple) = i .+ (0, 1)
move(::Val{:nw}, i::Tuple) = i .+ (-1, mod(i[1], 2) - 1)
move(::Val{:sw}, i::Tuple) = i .+ ( 1, mod(i[1], 2) - 1)
move(::Val{:ne}, i::Tuple) = i .+ (-1, mod(i[1], 2))
move(::Val{:se}, i::Tuple) = i .+ ( 1, mod(i[1], 2))

q1(paths) = length(flip_all(paths))

function flip_all(paths::AbstractVector)::Set
    blacks = Set{Tuple}()
    for tile in move.(paths)
        (tile ∈ blacks ? delete! : push!)(blacks, tile)
    end
    return blacks
end

function q2(paths::AbstractVector; steps::Integer = 100)::Int
    blacks = flip_all(paths)
    for _ in 1:steps
        step!(blacks)
    end
    return length(blacks)
end

function neighbours(pos::Tuple)::Vector{Tuple}
    return map(dir -> move(dir, pos), [:w, :nw, :ne, :e, :se, :sw])
end

function step!(blacks::Set)
    (minrow, maxrow) = blacks .|> first |> extrema
    (mincol, maxcol) = blacks .|> last  |> extrema
    tl = CartesianIndex(minrow - 1, mincol - 1)
    br = CartesianIndex(maxrow + 1, maxcol + 1)
    whites = setdiff(Tuple.(tl:br), blacks)

    to_white = [tile for tile in blacks
                     if count(∈(blacks), neighbours(tile)) ∉ [1, 2]]
    to_black = [tile for tile in whites
                     if count(∈(blacks), neighbours(tile)) == 2]

    union!(setdiff!(blacks, to_white), to_black)
end

@assert move("nwwswee") == (0, 0)
@assert move("esew") == (1, 0)

let v = parse_input(test)
    @assert q1(v) == 10
    @assert q2(v) == 2208
end

v = parse_input(read("day24.in", String))
println("Q1: ", q1(v))
println("Q2: ", q2(v))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
# Unit test functions
function unit_test_selectionsort_empty()
    a = Vector{Float64}(undef, 0)
    @test selectionsort!(a) == a

    return nothing
end

function unit_test_selectionsort_one()
    a = [10.0]
    @test selectionsort!(a) == a

    a = [10]
    @test selectionsort!(a) == a

    return nothing
end

function unit_test_selectionsort_many()
    # known example list
    a = [10.0, 1.0, 10.0, 3.0, 5.0, 6.0]
    as = [1.0, 3.0, 5.0, 6.0, 10.0, 10.0]
    @test selectionsort!(a) == as

    return nothing
end

function unit_test_selectionsort_real()
    n = 10^3
    # random example list
    a = rand((0:n), n)
    as = sort(a)
    @test selectionsort!(a) == as

    return nothing
end

function unit_test_selectionsort_string()
    # known example list
    a = "A search has no name."
    as = sort(a) # use built in for answer
    @test selectionsort!(a, "A") == as

    return nothing
end

function unit_test_insertionsort_empty()
    a = Vector{Float64}(undef, 0)
    @test insertionsort!(a) == a

    return nothing
end

function unit_test_insertionsort_one()
    a = [10.0]
    @test insertionsort!(a) == a

    a = [10]
    @test insertionsort!(a) == a

    return nothing
end

function unit_test_insertionsort_many()
    # known example list
    a = [10.0, 1.0, 10.0, 3.0, 5.0, 6.0]
    as = [1.0, 3.0, 5.0, 6.0, 10.0, 10.0]
    @test insertionsort!(a) == as

    return nothing
end

function unit_test_insertionsort_real()
    n = 10^3
    # random example list
    a = rand((0:n), n)
    as = sort(a)
    @test insertionsort!(a) == as

    return nothing
end

function unit_test_insertionsort_string()
    # known example list
    a = "A search has no name."
    as = sort(a) # use built in for answer
    @test insertionsort!(a, "A") == as

    return nothing
end

function unit_test_mergesort_empty()
    a = Vector{Float64}(undef, 0)
    @test mergesort(a) == a

    return nothing
end

function unit_test_mergesort_one()
    a = [10.0]
    @test mergesort(a) == a

    a = [10]
    @test mergesort(a) == a

    return nothing
end

function unit_test_mergesort_many()
    # known example list
    a = [10.0, 1.0, 10.0, 3.0, 5.0, 6.0]
    as = [1.0, 3.0, 5.0, 6.0, 10.0, 10.0]
    @test mergesort(a) == as

    return nothing
end

function unit_test_mergesort_real()
    n = 10^4
    # random example list
    a = rand((0:n), n)
    as = sort(a)
    @test mergesort(a) == as

    return nothing
end

function unit_test_mergesort_string()
    # known example list
    a = "A search has no name."
    as = sort(a) # use built in for answer
    @test mergesort(a, "A") == as

    return nothing
end

function unit_test_quicksort_empty()
    a = Vector{Float64}(undef, 0)
    @test quicksort!(a) == a

    return nothing
end

function unit_test_quicksort_one()
    a = [10.0]
    @test quicksort!(a) == a

    a = [10]
    @test quicksort!(a) == a

    return nothing
end

function unit_test_quicksort_many()
    # known example list
    a = [10.0, 1.0, 10.0, 3.0, 5.0, 6.0]
    as = [1.0, 3.0, 5.0, 6.0, 10.0, 10.0]
    @test quicksort!(a) == as

    return nothing
end

function unit_test_quicksort_real()
    n = 10^4
    # random example list
    a = rand((0:n), n)
    as = sort(a)
    @test quicksort!(a) == as

    return nothing
end

function unit_test_quicksort_string()
    # known example list
    a = "A search has no name."
    as = sort(a) # use built in for answer
    @test quicksort!(a, "A") == as

    return nothing
end

# Run the unit tests
@testset "Selection sort" begin
    unit_test_selectionsort_empty()
    unit_test_selectionsort_one()
    unit_test_selectionsort_many()
    unit_test_selectionsort_real()
    # unit_test_selectionsort_string()
end

@testset "Insertion sort" begin
    unit_test_insertionsort_empty()
    unit_test_insertionsort_one()
    unit_test_insertionsort_many()
    unit_test_insertionsort_real()
    # unit_test_insertionsort_string()
end

@testset "Merge sort " begin
    unit_test_mergesort_empty()
    unit_test_mergesort_one()
    unit_test_mergesort_many()
    unit_test_mergesort_real()
    # unit_test_mergesort_string()
end

@testset "Quicksort " begin
    unit_test_quicksort_empty()
    unit_test_quicksort_one()
    unit_test_quicksort_many()
    unit_test_quicksort_real()
    # unit_test_quicksort_string()
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# Run a local search algorithm to optimise ambulance deployment, given some objective function (objFn).
# The local search is simple hill climbing, and the neighbourhood of a solution
# (deployment) involves moving one ambulance from its station to another station.
# The local search stops when a local optimum (as evaluated by simulating) is found.
# Ambulances are assumed to be identical, and station capacities are ignored.

using JEMSS
using Random
using Statistics
using StatsFuns
using Base.Iterators: Stateful

const StationsNumAmbs = Vector{Int} # type alias

function deploymentLocalSearch!(sim::Simulation, deployments::Vector{Deployment}; outputFolder::String = "")
	isdir(outputFolder) || mkdir(outputFolder)
	@assert(isdir(outputFolder))
	@assert(!isempty(deployments))

	# parameters
	global solFilename = "$outputFolder/solutions.csv" # final solutions from each local search
	global deploymentsOutputFilename = "$outputFolder/deployments.csv" # solutions, stored as deployments
	global logFilename = "$outputFolder/log.csv" # log search progress to file
	global sense = :max # :min or :max; direction of optimisation for objective function
	global conf = 0.95 # statistical confidence level
	global doPrint = true
	# parallel (multi-threading):
	global parallel = false
	global numThreads = 1

	# some parameter checks
	@assert(sense == :min || sense == :max)
	@assert(0 <= conf < 1) # should be between 0 and 1, but not equal to 1 otherwise confidence interval is infinite
	@assert(isa(parallel, Bool))
	@assert(numThreads >= 1)

	global nullObjVal = MeanAndHalfWidth(NaN,NaN)

	# keep track of station ambulance counts tried, and their objective values
	global stationsNumAmbsObjVal = Dict{StationsNumAmbs, MeanAndHalfWidth}()
	global stationsNumAmbsStats = Dict{StationsNumAmbs, Dict{String, Any}}()

	# save the locally optimal station ambulance counts found for each search
	global stationsNumAmbsSols = Vector{StationsNumAmbs}() # stationsNumAmbsSols[i] is solution for deployments[i]

	repeatedLocalSearch!(sim, deployments)
end

function objValLookup(stationsNumAmbs::StationsNumAmbs)::MeanAndHalfWidth
	return get(stationsNumAmbsObjVal, stationsNumAmbs, nullObjVal) # return objective value if found, otherwise nullObjVal
end

function meanAndHalfWidth(x::Vector{T}; conf::Float = conf)::MeanAndHalfWidth where T <: Real
	x = convert(Vector{Float}, x)
	return MeanAndHalfWidth(mean(x), tDistrHalfWidth(x; conf = conf))
end

# Objective value for a single period
function objFn(period::SimPeriodStats)::Float
	return period.call.numResponsesInTime / period.call.numCalls # fraction of calls reached in time
end

# For completed simulation replications, calculate and return the objective function value mean and half-width of the mean.
function objFn(sim::Simulation)::MeanAndHalfWidth
	@assert(all(rep -> rep.complete, sim.reps))
	periods = getRepsPeriodStatsList(sim.reps)
	return meanAndHalfWidth([objFn(p) for p in periods])
end

# apply stationsNumAmbs to sim replications, and return objective value (and half-width) and whether lookup was used
# uses look-up if possible
# mutates: sim, stationsNumAmbsObjVal
function simObjVal!(sim::Simulation, stationsNumAmbs::StationsNumAmbs)::Tuple{MeanAndHalfWidth,Bool}
	global stationsNumAmbsObjVal, parallel, numThreads
	objVal = objValLookup(stationsNumAmbs)
	if objVal != nullObjVal return (objVal, true) end
	function runRep!(rep::Simulation)
		resetRep!(rep)
		applyStationsNumAmbs!(rep, stationsNumAmbs)
		simulateRep!(rep)
		@assert(getStationsNumAmbs(rep) == stationsNumAmbs) # check that rep had the current deployment applied
	end
	if parallel == true
		runParallel!(runRep!, sim.reps...; numThreads = numThreads)
	else
		for rep in sim.reps runRep!(rep) end
	end
	objVal = objFn(sim)
	stationsNumAmbsObjVal[stationsNumAmbs] = objVal
	return (objVal, false)
end

# print number of ambulances at each station
function printStationsNumAmbs(stationsNumAmbs::StationsNumAmbs)
	for i = 1:length(stationsNumAmbs)
		println("station ", i, ": ", stationsNumAmbs[i], " ambulance(s)")
	end
end

function getSimStats(sim::Simulation)
	# some sim statistics to write to log file
	global periodDuration, simStatsKeys, simStatsEmpty
	if all(r -> r.complete, sim.reps)
		periods = getRepsPeriodStatsList(sim.reps)
		@assert(all(p -> isapprox(p.duration, periodDuration), periods))
		stats = flatten(statsDictFromPeriodStatsList(periods))
		return merge(Dict(["$(key)_mean" => stats[key].mean for key in simStatsKeys]), Dict(["$(key)_halfWidth" => stats[key].halfWidth for key in simStatsKeys]))
	end
	return deepcopy(simStatsEmpty)
end

global simStatsEmpty = flatten(statsDictFromPeriodStatsList(SimPeriodStats[]))
global simStatsKeys = collect(keys(simStatsEmpty))
global simStatsEmpty = merge(Dict(["$(key)_mean" => NaN for key in simStatsKeys]), Dict(["$(key)_halfWidth" => NaN for key in simStatsKeys]))

global solFileHeader = vcat(["search", "objVal", "objValHalfWidth"], sort(collect(keys(simStatsEmpty)))) # ... and stationsNumAmbs
global logFileHeader = vcat(["search", "iter", "i", "j", "usedLookup", "usedMove", "searchDurationSeconds", "objVal", "objValHalfWidth", "bestObjVal"], sort(collect(keys(simStatsEmpty)))) # ... and stationsNumAmbs
global logFileDict = Dict{String,Any}([s => "" for s in logFileHeader])

# mutates: file
function fileWriteDlmLine!(file::IOStream, fileHeader::Vector{String}, data::Dict{String,T}, stationsNumAmbs::StationsNumAmbs) where T <: Any
	line = []
	for header in fileHeader
		push!(line, data[header])
	end
	line = vcat(line, stationsNumAmbs)
	writeDlmLine!(file, line...)
	flush(file)
end

# perform local search, starting at each of the deployments provided
# mutates: sim
function repeatedLocalSearch!(sim::Simulation, deployments::Vector{Deployment})
	global doPrint

	# open files for writing solution
	global solFilename, logFilename
	solFile = open(solFilename, "w")
	logFile = open(logFilename, "w")

	# write misc data to files
	global warmUpDuration, periodDuration, conf
	numSearches = length(deployments)
	miscTable = Table("miscData", ["numAmbs", "numStations", "numSearches", "conf"];
		rows = [[sim.numAmbs, sim.numStations, numSearches, conf]]) # note that numCalls may vary between sim replications
	repsTable = Table("repsData", ["numReps", "warmUpDuration", "periodDuration"];
		rows = [[sim.numReps, warmUpDuration, periodDuration]])
	writeTablesToFile!(logFile, [miscTable, repsTable])
	writeTablesToFile!(solFile, [miscTable, repsTable])

	# write file headers
	writeDlmLine!(solFile, "solutions")
	writeDlmLine!(logFile, "log")
	writeDlmLine!(solFile, solFileHeader..., ["station_$i numAmbs" for i = 1:sim.numStations]...)
	writeDlmLine!(logFile, logFileHeader..., ["station_$i numAmbs" for i = 1:sim.numStations]...)
	flush(solFile)
	flush(logFile)

	for i = 1:numSearches
		logFileDict["search"] = i

		# initial deployment
		stationsNumAmbs = deploymentToStationsNumAmbs(deployments[i], sim.numStations)

		# perform local search until no improvements can be made
		doPrint && println()
		doPrint && println("Search ", i, " (of ", numSearches, ")")
		doPrint && printStationsNumAmbs(stationsNumAmbs)
		stationsNumAmbs = localSearch!(sim, stationsNumAmbs, logFile)
		doPrint && println("Finished local search $i (of $numSearches); solution:")
		doPrint && printStationsNumAmbs(stationsNumAmbs)

		# save best station ambulance count found for this search
		global stationsNumAmbsSols
		push!(stationsNumAmbsSols, stationsNumAmbs)

		# write solution to file
		objValMeanAndHalfWidth = objValLookup(stationsNumAmbs)
		solFileDict = Dict("search" => i, "objVal" => objValMeanAndHalfWidth.mean, "objValHalfWidth" => objValMeanAndHalfWidth.halfWidth)
		merge!(solFileDict, stationsNumAmbsStats[stationsNumAmbs])
		fileWriteDlmLine!(solFile, solFileHeader, solFileDict, stationsNumAmbs)

		global deploymentsOutputFilename
		writeDeploymentsFile(deploymentsOutputFilename, map(stationsNumAmbsToDeployment, stationsNumAmbsSols), sim.numStations)
	end

	# print out all results from each finished local search
	doPrint && println()
	doPrint && println("Station ambulance counts from completed local searches:")
	for i = 1:numSearches
		doPrint && println()
		doPrint && println("Search: ", i)
		stationsNumAmbs = stationsNumAmbsSols[i]
		doPrint && printStationsNumAmbs(stationsNumAmbs)
		doPrint && println("objective value = ", objValLookup(stationsNumAmbs).mean)
	end

	close(solFile)
	close(logFile)
end

# Local search of ambulance deployment, starting at stationsNumAmbs.
# Move ambulances between stations, making changes that optimise
# the objective function (objFn) for the given sense (:min or :max).
# Continue search until no improvement can be made.
# mutates: sim, logFile
function localSearch!(sim::Simulation, stationsNumAmbs::StationsNumAmbs, logFile::IOStream)::StationsNumAmbs

	global logFileHeader, logFileDict, stationsNumAmbsStats, sense, doPrint

	# track iterations
	iter = 1 # count iterations performed
	startTime = time()
	getSearchDuration() = round(time() - startTime, digits = 2)

	# calculate objective value for starting point
	objValMeanAndHalfWidth, usedLookup = simObjVal!(sim, stationsNumAmbs)
	(objVal, objValHalfWidth) = (objValMeanAndHalfWidth.mean, objValMeanAndHalfWidth.halfWidth)
	doPrint && println("starting objective value: ", objVal)
	bestObjVal = objVal # current best objective value
	bestStationsNumAmbs = copy(stationsNumAmbs)

	# write starting point
	for (key, value) in logFileDict
		if key != "search"
			logFileDict[key] = "" # reset value
		end
	end
	stats = getSimStats(sim)
	if !usedLookup stationsNumAmbsStats[stationsNumAmbs] = stats end
	merge!(logFileDict, stats)
	merge!(logFileDict, Dict("iter" => iter, "usedLookup" => Int(usedLookup), "searchDurationSeconds" => getSearchDuration(),
		"objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal))
	fileWriteDlmLine!(logFile, logFileHeader, logFileDict, stationsNumAmbs)

	# local search
	doPrint && println("starting local search")
	endPoint = nothing
	iter += 1
	numStations = sim.numStations # shorthand
	while true
		for i = 1:numStations, j = [i+1:numStations; 1:i-1;]
			newBestFound = false
			if bestStationsNumAmbs[i] > 0
				# test moving an ambulance from station i to j

				doPrint && print("move: i = ", i, ", j = ", j)

				# change station ambulance counts, later keep if improvement made
				stationsNumAmbs = copy(bestStationsNumAmbs)
				stationsNumAmbs[i] -= 1
				stationsNumAmbs[j] += 1

				objValMeanAndHalfWidth, usedLookup = simObjVal!(sim, stationsNumAmbs)
				(objVal, objValHalfWidth) = (objValMeanAndHalfWidth.mean, objValMeanAndHalfWidth.halfWidth)
				if !usedLookup
					stats = stationsNumAmbsStats[stationsNumAmbs] = getSimStats(sim)
					merge!(logFileDict, stats)
				else
					merge!(logFileDict, stationsNumAmbsStats[stationsNumAmbs])
					doPrint && print("; done")
				end

				doPrint && println("; objective value: ", objVal, " (best: ", bestObjVal, ")")

				usedMove = false
				if (sense == :max && bestObjVal < objVal) || (sense == :min && bestObjVal > objVal)
					bestObjVal = objVal
					# keep change (move ambulance from station i to j)
					bestStationsNumAmbs = stationsNumAmbs
					usedMove = true
					newBestFound = true

					doPrint && println("moved: i = ", i, ", j = ", j, "; new objective value: ", bestObjVal)
				end

				merge!(logFileDict, Dict("iter" => iter, "i" => i, "j" => j, "usedLookup" => Int(usedLookup), "usedMove" => Int(usedMove),
					"searchDurationSeconds" => getSearchDuration(), "objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal))
				fileWriteDlmLine!(logFile, logFileHeader, logFileDict, stationsNumAmbs)

				iter += 1
			end

			if newBestFound || endPoint == nothing
				endPoint = (i,j)
			elseif endPoint == (i,j)
				doPrint && println("Found local optimum.")
				return bestStationsNumAmbs
			end
		end
	end
end

# mutates: sim
function setSimStatsCapture!(sim::Simulation, periodDurationsIter::Stateful, warmUpDuration::Float)
	stats = sim.stats = deepcopy(sim.backup.stats)
	stats.doCapture = true
	(stats.periodDurationsIter, stats.warmUpDuration) = (periodDurationsIter, warmUpDuration)
	stats.nextCaptureTime = sim.startTime + (stats.warmUpDuration > 0 ? stats.warmUpDuration : first(stats.periodDurationsIter))
	sim.backup.stats = deepcopy(sim.stats)
end

# parameters
configFilename = "sim_config.xml"
outputFolder = "output"
numSearches = 1
deploymentRng = Random.MersenneTwister(0) # useful for reproducing results, if using random restarts
# replication parameters
warmUpDuration = 0.1
periodDuration = 1.0

# parameter checks
@assert(numSearches >= 1)
@assert(warmUpDuration >= 0)
@assert(periodDuration > 0)

# initialise sim
t = time()
@assert(isfile(configFilename))
println("Initialising simulation from config: ", configFilename)
sim = initSim(configFilename, createBackup = true)
@assert(sim.numReps >= 2, "Need at least two simulation replications (sim.reps) for statistics.")

# set sim stats capturing
setSimStatsCapture!(sim, Stateful(periodDuration), warmUpDuration)
periodEndTime = sim.startTime + warmUpDuration + periodDuration
for rep in sim.reps
	# statistics should stop collecting before last call arrives (if not earlier), to avoid cool-down period
	@assert(periodEndTime <= rep.calls[end].arrivalTime)
end

makeRepsRunnable!(sim) # copy changes to sim (stats) to sim.reps

# generate deployments
deployments = makeRandDeployments(sim, numSearches; rng = deploymentRng)

# run
deploymentLocalSearch!(sim, deployments, outputFolder = outputFolder)
println("total runtime: ", round(time()-t, digits = 2), " seconds")


text:
function ordinal_string( n )  #=

converts a given integer into its ordinal-number string. It is assumed that the argument is an integer, and no error handling is implemented.

Author: Taku Noda
Started on Oct. 31, 2019
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

r = n % 10
if r == 1
  str = string( n )*"st"
elseif r == 2
  str = string( n )*"nd"
elseif r == 3
  str = string( n )*"rd"
else
  str = string( n )*"th"
end

return str

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
module RoboChariot

export
  version

function version()
  v"0.0.0+"
end

end


text:
#   Unit "crossspectra" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements Welch cross-spectral estimates using FFTW

#############################################################################
# Generic constructors of CrossSpectra objects:
#
# Enable construction giving only `y`, `sr`, `wl`, `DC`, `taper`
# and `non-linear` argument.
# `flabels` is generated automatically, `smoothing` is set to `noSmoother`
# and `tril` is set to false
CrossSpectra(y, sr, wl, DC, taper, nonlinear) =
        CrossSpectra(y, sr, wl, DC, taper, fdf(sr, wl; DC=DC), nonlinear, noSmoother, false)
#
# As above, but setting by default also `DC` and `nonlinear` to false
CrossSpectra(y, sr, wl, taper) =
    CrossSpectra(y, sr, wl, false, taper, fdf(sr, wl; DC=false), false, noSmoother, false)
#############################################################################

"""
```
(1)
function crossSpectra( X    :: Matrix{T},
                       sr   :: Int,
                       wl   :: Int;
                  tapering  :: Union{Taper, TaperKind} = harris4,
                  planner   :: Planner                 = getplanner,
                  slide     :: Int                     = 0,
                  DC        :: Bool                    = false,
                  nonlinear :: Bool                    = false,
                  smoothing :: Smoother                = noSmoother,
                  tril      :: Bool                    = false,
                  ⏩       :: Bool                    = true) where T<:Real

(2)
function crossSpectra( 𝐗    :: Vector{Matrix{T}},
              < same argument sr, ..., ⏩ of method (1) > where T<:Real

```

(1)

Construct a [CrossSpectra](@ref) objects from real multivariate data
using the Welch method.

Given sampling rate `sr` and epoch length `wl`, compute the cross-spectral
matrices of dimension ``n``x``n`` for a multivariate data matrix
`X` of dimension ``t``x``n``, where ``t`` is the number of samples (rows) and
`n>1` is the number of time-series (columns).

The cross-spectral matrices are hold in the `.y` vector field of the created
object. The length of `.y` depends upon
the `wl` argument and `DC` optional keyword argument (see below).

**Optional Keyword Arguments**:

`sr`, `wl`, `tapering`, `planner` and `slide` have the same meaning as for the
[`spectra`](@ref) function.

`DC`: if true the cross-spectral matrix of the DC level is returned
in the first position of `y` (see the fields of the [CrossSpectra](@ref)
object), otherwise (default) the matrices in `y` start with the first positive
discrete frequency, that is, ``sr/wl`` Hz.

`nonlinear`: if true, the amplitude information is eliminated from the
DFT coefficients, that is, they are normalized by their modulus before being
averaged.
This leads to non-linear estimates (Congedo, 2018;
Pascual-Marqui 2007) where the diagonal elements of the cross-spectral matrices
(the spectra) are 1.0 for all frequencies. By default, it is false.

`smoothing`: apply a smoothing function of type [Smoother](@ref)
to the cross-spectral matrices across frequencies.
By default no smoothing is applied.

`tril`: if false (default), the whole cross-spectra matrices will be computed,
otherwise only their lower triangular part (see below).

`⏩`: if true (default), the method is run in multi-threaded mode
across the series in `X` if the number of series is at least twice
the number of threads Julia is instructed to use. See [Threads](@ref).

**Return**

If `tril` is false (default), the output is of type `Array{Hermitian,1}`,
which is the `ℍVector` type used in package [PosDefManifold](https://github.com/Marco-Congedo/PosDefManifold.jl).
Since cross-spectral estimates are Hermitian
positive definite, they can be straightaway used as
argument to PosDefManifold's functions, e.g., for computing matrix moves on
[geodesics](https://marco-congedo.github.io/PosDefManifold.jl/dev/riemannianGeometry/#Geodesic-equations-1),
matrix [distances](https://marco-congedo.github.io/PosDefManifold.jl/dev/riemannianGeometry/#PosDefManifold.distance),
etc. and the the whole vector output to compute
matrix [means](https://marco-congedo.github.io/PosDefManifold.jl/dev/riemannianGeometry/#Means-1),
[spectral embedding](https://marco-congedo.github.io/PosDefManifold.jl/dev/riemannianGeometry/#PosDefManifold.spectralEmbedding)
and more.

If `tril` is true, the output is of type `Array{LowerTriangular,1}`,
which is the `𝕃Vector` type used in PosDefManifold, that is, only the
lower triangle of the cross-spectra is computed in order to save time
and memory.

(2)

Construct a [CrossSpectraVector](@ref) object from a vector of
real multivariate data matrices. Compute the cross-spectral matrices
using the Welch method as per method (1) for all ``k`` data
matrices in `𝐗`.

The ``k`` matrices in `𝐗` must have the same number of columns
(i.e., the same number of time-series), but may have any number of (at least
`wl`) rows (samples).
All other arguments have the same meaning as in method (1),
with the following difference:

`⏩`: if true (default), the method is run in multi-threaded mode across the
``k`` data matrices if ``k`` is at least twice the number of threads Julia
is instructed to use, otherwise this method attempts to run each cross-spectral
estimation in multi-threaded mode across series as per method (1).
See [Threads](@ref).

If a `planner` is not explicitly passed as an argument,
the FFTW plan is computed once and applied for all cross-spectral
estimations.

**Return**

If `tril` is false, the output is of type `Array{Array{Hermitian,1},1}`,
which is the `ℍVector₂` type used in
[PosDefManifold](https://github.com/Marco-Congedo/PosDefManifold.jl).

If `tril` is true, the output is of type `Array{Array{LowerTriangular,1},1}`,
which is the `𝕃Vector₂` type used in PosDefManifold.

**See**: [CrossSpectra](@ref).

**See also**: [`spectra`](@ref), [`coherence`](@ref).

**Examples**:
```
using FourierAnalysis, Plots, LinearAlgebra

function generateSomeData(sr::Int, t::Int; noise::Real=1.)
    # four sinusoids of length t samples and sr sampling rate
    # peak amplitude: 0.7, 0.6, 0.5, 0.4
    # frequency:        5,   7,  13,  27
    # phase:            0, π/4, π/2,   π
    v1=sinusoidal(0.7, 5,  sr, t, 0)
    v2=sinusoidal(0.6, 7,  sr, t, π/4)
    v3=sinusoidal(0.5, 13, sr, t, π/2)
    v4=sinusoidal(0.4, 27, sr, t, π)
    return hcat(v1, v2, v3, v4) + (randn(t, 4)*noise)
end

sr, wl, t = 128, 512, 8192

# (1)

X=generateSomeData(sr, t) # multivariate data matrix 8192x4

# cross-spectra using default harris4 tapering window
S=crossSpectra(X, sr, wl)

# check the cross-spectral matrix at frequency 5Hz
S.y[f2b(5, sr, wl)]

# check only the diagonal part of this matrix as a vector
diag(Diagonal(S.y[f2b(5, sr, wl)]))

# cross-spectra using hann tapering window
S=crossSpectra(X, sr, wl; tapering=hann)

# using Slepian's multi-tapering
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl))

# compute non-linear cross-spectra
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl), nonlinear=true)

# compute only the lower triangle of cross-spectral matrices
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl), tril=true)

# smooth a-posteriori the cross-spectra
S2=smooth(blackmanSmoother, S)

# or compute cross-spectra already smoothed
S=crossSpectra(X, sr, wl;
               tapering=slepians(sr, wl), tril=true, smoothing=blackmanSmoother)

# mean cross-spectral matrix in 8Hz-12Hz range
M=mean(S, (8, 12)) # or also M=mean(S, (8.0, 12.0))

# extract all cross-spectral matrices in 8Hz-12Hz range
E=extract(S, (8, 12))

# cross-spectral matrices averaged in 2Hz band-pass regions
B=bands(S, 2)

# Get the spectra from a CrossSpectra object
PowerSpectra=Spectra(S)

# Get the amplitude spectra from a CrossSpectra object
AmpSpectra=Spectra(S, func=√)

# Get the log10-spectra from a CrossSpectra object
log10Spectra=Spectra(S, func=log10)

# plot the spectra (see recipes.jl)
plot(AmpSpectra; fmax=32, xspace=4, ytitle="Amplitude")

# (2)
# generate 3 multivariate data matrices 8192x4
X=[generateSomeData(sr, t) for i=1:3]

# The examples here below use exactly the same syntax as the previous method.
# However, since the input here is a vector of data matrices
# and not a single data matrix, the examples here below create a vector
# of the object created by the examples above.
# For example:

# cross-spectra using the default harris4 tapering window
# this creates a CrossSpectraVector object
S=crossSpectra(X, sr, wl)

# check the cross-spectral matrix at fr. 5Hz for the first CrossSpectra object
S[1].y[f2b(5, sr, wl)]

# check only the diagonal part of this matrix as a vector
diag(Diagonal(S[1].y[f2b(5, sr, wl)]))

# cross-spectra using Hamming's tapering window
S=crossSpectra(X, sr, wl; tapering=hamming)

# using Slepian's multi-tapering
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl))

# compute non-linear cross-spectra
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl), nonlinear=true)

# compute only the lower triangle of cross-spectral matrices
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl), tril=true)

# smooth a-posteriori all CrossSpectra objects in S
S2=smooth(blackmanSmoother, S)

# or compute them all already smoothed
S=crossSpectra(X, sr, wl; tapering=parzen, smoothing=blackmanSmoother)

# mean cross-spectral matrix in 8Hz-12Hz range for all CrossSpectra (CS) objects
M=mean(S, (8, 12)) # or also M=mean(S, (8.0, 12.0))

# grand-average mean of the above across all CS objects
meanM=mean(mean(S, (8, 12)))

# extract all cross-spectral matrices in 8Hz-12Hz range for all CS objects
E=extract(S, (8, 12))

# grand average of cross-spectral matrices in 8Hz-12Hz range for all CS objects
meanE=mean(extract(S, (8, 12)))

# cross-spectral matrices averaged in 2Hz band-pass regions for all CS objects
B=bands(S, 2)

# Get and plot the spectra from a CrossSpectra object
plot(Spectra(S[1]); fmax=32, xspace=4)

# Pre-compute a FFT planner and pass it as argument
# (this interesting if the function is to be called repeatedly).
plan=Planner(plan_exhaustive, 10.0, wl, eltype(X[1])) # wait 10s
S=crossSpectra(X, sr, wl; planner=plan)

# how faster is this?
using BenchmarkTools
@benchmark(crossSpectra(X, sr, wl))
@benchmark(crossSpectra(X, sr, wl; planner=plan))
...
...

```
"""
function crossSpectra( X    :: Matrix{T},
                       sr   :: Int,
                       wl   :: Int;
                  tapering  :: Union{Taper, TaperKind} = harris4,
                  planner   :: Planner                 = getplanner,
                  slide     :: Int                     = 0,
                  DC        :: Bool                    = false,
                  nonlinear :: Bool                    = false,
                  smoothing :: Smoother                = noSmoother,
                  tril      :: Bool                    = false,
                  ⏩       :: Bool                    = true) where T<:Real

    # function to set all parameters
    n, 𝚙, 𝜏, s, e, 𝜎, η, wl½, wl½_, cT = _paramSP!(X, tapering, wl, planner, slide, DC)

    # allocate memory for DFTs (𝙵)
    𝙵 = [zeros(cT, wl½_, n) for i=1:𝜏.n]

    # allocate memory for cross-spectra and initialize them with zeros (𝙎)
    𝙎 = Vector{LowerTriangular}([LowerTriangular(zeros(cT, n, n)) for i=1:wl½_])

    # compute the descrete Fourier transform (DFT) of a Vector or of the column vectors
    # of matrix `X` and put it into the columns of `𝙵`. The DFT is computed on the epoch
    # of duration `t` samples starting at sample `at`. It does not return anything.
    # `𝙵` is a `n𝜏`-Vector of DFT matrices, holding a DFT matrix (# discr. freq. x time-series),
    # where `n𝜏` is the number of tapering window (>1 only for Slepian tapers).
    # ARGUMENTS:
    # `η` is a calibration Vector to be applied to the FFT coefficients, computed
    #    by function `_fftNormalization`. It is applied element-wise to the FFT
    #    coefficient vector(s) to obtain spectral estimates in correct physical units.
    #    It is applied only if `nonlinear` is false.
    # `n` (Int) is the number of columns in `X`; it must be 1 if `X` is a vector.
    # `𝜏` is a tapering window, ``rectangular` if no tapering is applied
    # `𝚙` is the FFTW plan performing the FFT.
    # `DC` (bool): if true the DC level(s) is(are) computed and stored in the first matrix of 𝚂
    _dft!(at, 𝚗) =
        if nonlinear
            𝜏.kind == rectangular ? ( 𝙵[1][:, 𝚗] = ((𝚙*X[at:at+wl-1, 𝚗])[1+!DC:end]) ) : # no tapering
                                    ( for i=1:𝜏.n 𝙵[i][:, 𝚗] = (𝚙*(𝜏.y[:, i].*X[at:at+wl-1, 𝚗]))[1+!DC:end] end ) # with tapering
        else
            𝜏.kind == rectangular ? ( 𝙵[1][:, 𝚗] = η.*((𝚙*X[at:at+wl-1, 𝚗])[1+!DC:end]) ) : # no tapering
                                    ( for i=1:𝜏.n 𝙵[i][:, 𝚗] = η.*(𝚙*(𝜏.y[:, i].*X[at:at+wl-1, 𝚗]))[1+!DC:end] end ) # with tapering
        end

    # function computing DFTs (discrete Fourier transform)
    dft!(𝚎, 𝚗) = _dft!((𝚎*s)+1, 𝚗)

    # function cumulating the avarage cross-spectra. The mean is over tapering windows (>1 only for Slepians)
    cs!(𝚏) = for j=1:n, i=j:n 𝙎[𝚏][i, j] += mean(𝙵[l][𝚏, i] * conj(𝙵[l][𝚏, j]) for l=1:𝜏.n) end

    # average cross-spectra across sliding windows
    for 𝚎=0:e-1
        _thread(⏩, n) ? (@threads for 𝚗=1:n dft!(𝚎, 𝚗) end) : (for 𝚗=1:n dft!(𝚎, 𝚗) end)
        nonlinear ? (for i=1:length(𝙵) 𝙵[i]./=(abs.(𝙵[i])*𝜎) end) : nothing
        _thread(⏩, wl½_) ? (@threads for 𝚏=1:wl½_ cs!(𝚏) end) : (for 𝚏=1:wl½_ cs!(𝚏) end)
    end

    # if `tril`=false return the full matrices, otherwise the lower-triangular matrices
    # Typecasting output: 𝔸=Array; 𝕄=Matrix; ℍ=Hermitian; 𝕃=LowerTriangular
    H=Hermitian; V=Vector; M=Matrix; CS=CrossSpectra; sm=smoothing
    args=(sr, wl, DC, taperinfo(𝜏), fdf(sr, wl; DC=DC), nonlinear, sm, tril)
    return tril ? smooth(sm, CS(𝙎, args...)) :
                  smooth(sm, CS(V{H}([H(M(𝙎[𝚏]), :L) for 𝚏=1:wl½_]), args...))
end

function crossSpectra( 𝐗    :: Vector{Matrix{T}},
                       sr   :: Int,
                       wl   :: Int;
                  tapering  :: Union{Taper, TaperKind} = harris4,
                  planner   :: Planner                 = getplanner,
                  slide     :: Int                     = 0,
                  DC        :: Bool                    = false,
                  nonlinear :: Bool                    = false,
                  smoothing :: Smoother                = noSmoother,
                  tril      :: Bool                    = false,
                  ⏩       :: Bool                    = true) where T<:Real

    # set some parameters. k=# input data matrices, n=# time series
    k, n, wl½, thr, cT = length(𝐗), size(𝐗[1], 2), (wl>>1), nthreads(), typeof(Complex(T(0)))

    # get the FFTW plan (for one time-series)
    planner ≠ getplanner ? plan=planner : plan = _Planner(plan_rfft(zeros(T, wl)))

    # get the tapering window if not provided
    tapering isa Taper ? 𝜏=tapering : 𝜏=taper(tapering, wl)

    # allocate memory for cross-spectra
    𝓢 = CrossSpectraVector(undef, k)

    # define internal function to compute cross-spectra for the ith input data matrix `𝐗[i]`
    cs!(i, ⏩) = 𝓢[i]=crossSpectra(𝐗[i],
                                   sr,
                                   wl;
                               tapering  = 𝜏,
                               slide     = slide,
                               DC        = DC,
                               planner   = plan,
                               nonlinear = nonlinear,
                               smoothing = smoothing,
                               tril      = tril,
                               ⏩       = ⏩)

    # threads over k if k is large, otherwise let crossSpectra function decide if to thread over n
    _thread(⏩, k) ? (@threads for i=1:k cs!(i, false) end) : (for i=1:k cs!(i, ⏩) end)

    return 𝓢
end

# ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, S::CrossSpectra)
    r=size(S.y[1], 1)
    l=length(S.flabels)
    ll=length(S.y)
    println(io, titleFont, "▦ ⋯ ▦ CrossSpectra type; $l freq. x $(r)²-matrices")
    #println(io, "□  □    □      □        □           □", defaultFont)
    println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
    println(io, "sampling rate    (.sr): $(S.sr)")
    println(io, "epoch length     (.wl): $(S.wl)")
    println(io, "DC level         (.DC): $(S.DC)")
    println(io, "taper kind (.tapering): $(S.taper)")
    println(io, "freq. lab.  (.flabels): $(l)-", typeof(S.flabels))
    println(io, "data              (.y): $(l)-", typeof(S.y))
    # println(io, S.y)
    println(io, "non-linear(.nonlinear): ", string(S.nonlinear))
    println(io, "smoother  (.smoothing): ", string(S.smoothing))
    println(io, "triang. form   (.tril): ", string(S.tril))
    ll≠l && @warn "number of frequency labels does not match the CrossSpectra data size" l ll
end

function Base.show(io::IO, ::MIME{Symbol("text/plain")}, 𝐒::CrossSpectraVector)
    println(io, titleFont, "▦ ⋯ ▦ CrossSpectraVector Type")
    println(io, titleFont, "  ⋯")
    println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
    println(io, "$(length(𝐒))-element Vector{CrossSpectra}")
end

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #


text:
mutable struct Friction{T} <: Contact{T}
    ainv3::Adjoint{T,SVector{3,T}}
    p::SVector{3,T}
    D::SMatrix{2,3,T,6}
    cf::T
    offset::SVector{3,T}
    b::SVector{2,T}

    function Friction(body::Body{T}, normal::AbstractVector, cf::Real; p::AbstractVector = zeros(3), offset::AbstractVector = zeros(3)) where T
        @assert cf>0

        # Derived from plane equation a*v1 + b*v2 + distance*v3 = p - offset
        V1, V2, V3 = orthogonalcols(normal) # gives two plane vectors and the original normal axis
        A = [V1 V2 V3]
        Ainv = inv(A)
        ainv3 = Ainv[3,SA[1; 2; 3]]'
        D = [V1 V2]'

        new{T}(ainv3, p, D, cf, offset, zeros(2)), body.id
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::Friction{T}) where {T}
    summary(io, constraint)
    println(io,"")
    println(io, " ainv3:  "*string(constraint.ainv3))
    println(io, " D:      "*string(constraint.D))
    println(io, " cf:     "*string(constraint.cf))
    println(io, " offset: "*string(constraint.offset))
end

@inline function calcb0(friction::Friction, m, Δt, x::AbstractVector, q::UnitQuaternion, vc, ωc, dyn)
    p = friction.p
    np2 = norm(p)^2
    vp = vc + vrotate(cross(ωc,p),q)
    if np2 > 0
        return friction.D * (I*m*(-vp)/Δt + dyn[SA[1;2;3]] + vrotate(cross(dyn[SA[4;5;6]]/2,p),q)/np2)
    else
        return friction.D * (I*m*(-vp)/Δt + dyn[SA[1;2;3]])
    end
end

@inline function frictionmap(friction::Friction, x::AbstractVector, q::UnitQuaternion)
    Fp = friction.D' * friction.b
    Fc = Fp
    τc = torqueFromForce(vrotate(Fp,inv(q)), friction.p) # in local coordinates
    return [Fc;2*τc]
end

@inline additionalforce(friction::Friction, x::AbstractVector, q::UnitQuaternion) = frictionmap(friction, x, q)

@inline function calcFrictionForce!(mechanism, ineqc, friction::Friction, i, body::Body)
    calcFrictionForce!(mechanism, friction, body, ineqc.γsol[2][i])
    return
end

##Friction
# TODO this might be wrong, will be done differently in future anyways
@inline function calcFrictionForce!(mechanism::Mechanism{T}, friction::Contact, body::Body, γ) where T
    cf = friction.cf
    D = friction.D
    state = body.state

    d = state.d
    v = state.vsol[2]
    ω = state.ωsol[2]
    state.vsol[2] = state.vc
    state.ωsol[2] = state.ωc
    dyn = dynamics(mechanism, body) + frictionmap(friction, posargsk(body.state)...)
    state.vsol[2] = v
    state.ωsol[2] = ω
    state.d = d

    b0 = calcb0(friction, body.m, mechanism.Δt, posargsk(body.state)..., state.vc, state.ωc, dyn)

    if norm(b0) > cf*γ
        friction.b = b0/norm(b0)*cf*γ
    else
        friction.b = b0
    end
    return
end


text:
import Random

import ALifeBenchmark
using Test

Random.seed!(1234)

include("analysis/neutral_networks_test.jl")
include("analysis/estimator_test.jl")

include("geb/inputs_tests.jl")
include("geb/actions_tests.jl")
include("geb/network_tests.jl")
include("geb/nodes_tests.jl")
include("geb/rules_tests.jl")
include("geb/utils_tests.jl")

include("tierra/scheduler_tests.jl")
include("tierra/instruction_tests.jl")
include("tierra/ancestor_tests.jl")
include("tierra/memory_utils_tests.jl")
include("tierra/memory_tests.jl")

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using StructuralCausalModels
#using StatsPlots
#gr(size=(700,600))

ProjDir = @__DIR__
cd(ProjDir) #do

df = DataFrame!(CSV.File(scm_path("..", "data", "WaffleDivorce.csv"), delim=';'));
df = DataFrame(
  :s => df[:, :South],
  :a => df[:, :MedianAgeMarriage],
  :m => df[:, :Marriage],
  :w => df[:, :WaffleHouses],
  :d => df[:, :Divorce]
)
first(df, 5) |> display

if isdefined(Main, :StatsPlots)
  StatsPlots.cornerplot(Array(df), label=names(df))
  savefig("$(ProjDir)/SR6.4.3.png")
end

d = OrderedDict(
  [:w, :m, :a] => :s,
  :d => [:a, :w, :m],
  :m => [:a]
);

dag = DAG("sr6_4_3", d, df=df);
show(dag)

fname = ProjDir * "/sr6.4.3.dot"
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

display(dag.s); println()

println("Basis set:")
bs = basis_set(dag)
display(bs); println()

t = shipley_test(dag)
display(t); println()

f = :w; s = :d;

e = d_separation(dag, [f], [s], c=[:m, :a])
println("d_separation($(dag.name), $f, $s, c=[:m, :a]) = $e\n")

#ap = all_paths(dag, f, s)
ap = all_paths(dag, f, s)
ap |> display
println()

bp = backdoor_paths(dag, ap, f)
bp |> display
println()

adjustmentsets = adjustment_sets(dag, :w, :d)
println("Adjustment sets:")
adjustmentsets |> display

#end


text:
@testset "Vulkan" begin
    include("vulkan/gen.jl")
end


text:
using PyPlot, PyCall

inset_locator = pyimport("mpl_toolkits.axes_grid1.inset_locator")
inset_axes = inset_locator.inset_axes

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[6, 3])

axins1 = inset_axes(ax1,
                    width="50%",  # width = 50% of parent_bbox width
                    height="5%",  # height : 5%
                    loc="upper right")

im1 = ax1.imshow([[1, 2], [2, 3]])
fig.colorbar(im1, cax=axins1, orientation="horizontal", ticks=[1, 2, 3])
axins1.xaxis.set_ticks_position("bottom")

axins2 = inset_axes(ax2,
                   width="5%",  # width = 5% of parent_bbox width
                   height="100%",  # height : 50%
                   loc="lower left",
                   bbox_to_anchor=(1.05, 0., 1, 1),
                   bbox_transform=ax2.transAxes,
                   borderpad=0,
                   )

# Controlling the placement of the inset axes is basically same as that
# of the legend. You may want to play with the borderpad value and
# the bbox_to_anchor coordinate.

im = ax2.imshow([[1, 2], [2, 3]])
fig.colorbar(im, cax=axins2, ticks=[1, 2, 3])

plt.show()


text:
using SQLite

if (!("test-faulty.sqlite" in readdir()))

    db = SQLite.DB("test-faulty.sqlite")

    SQLite.execute(
        db,
        """
        CREATE TABLE partey
        (
            party_id INTEGER NOT NULL PRIMARY KEY,
            party_shorthand TEXT
        );
        """
    )

    SQLite.execute(
        db,
        """
        CREATE TABLE statament
        (
            statement_id INTEGER NOT NULL PRIMARY KEY,
            statement TEXT
        );
        """
    )

    SQLite.execute(
        db,
        """
        CREATE TABLE opinioon
        (
            party_id INTEGER NOT NULL,
            statement_id INTEGER NOT NULL,
            position INTEGER,
            FOREIGN KEY(party_id) REFERENCES party(party_id),
            FOREIGN KEY(statement_id) REFERENCES statement(statement_id),
            PRIMARY KEY(party_id, statement_id)
        );
        """
    )

    SQLite.execute(
        db,
        """
        INSERT INTO partey
        VALUES
            (1, \"TP_1\"),
            (2, \"TP_2\"),
            (3, \"TP_3\");
        """
    )

    SQLite.execute(
        db,
        """
        INSERT INTO statament
        VALUES
            (1, \"This is the first test statement.\"),
            (2, \"This is the second test statement.\"),
            (3, \"This is the third test statement.\");
        """
    )

    SQLite.execute(
        db,
        """
        INSERT INTO opinioon
        VALUES
            (1, 1, 1),
            (1, 2, -1),
            (1, 3, 0),
            (2, 1, -1),
            (2, 2, 0),
            (2, 3, -1),
            (3, 1, 1),
            (3, 2, 1),
            (3, 3, 0);
        """
    )

end