text:
module mymod
    const g = 9.81
end

module mymod2
    const g = 9.81/6
end

struct Constants1 end
(c::Constants1)(::Val{s}) where s = getproperty(mymod, s)
cst1 = Constants1()
cst1(Val(:g))
@code_llvm (() -> 2*Constants1()(Val(:g)))()

f1(::Val{cst}) where cst = 2*cst(Val(:g))
f1(Val(cst1))
@code_llvm f1(Val(cst1))

struct Constants2 end
(c::Constants2)(::Val{s}) where s = getproperty(mymod2, s)
cst3 = Constants2()
f1(Val(cst3))
@code_llvm f1(Val(cst3))

struct ConstantsB end
(c::ConstantsB)(::Val{s}) where s = getproperty(Bigleaf.BigleafConstantsDef, s)
cstB = ConstantsB()
cstB(Val(:g))
f1(Val(cstB))
@code_llvm f1(Val(cstB))

using Bigleaf
f3(;cst=BigleafConstants()) = 2 * cst.g
f3()
@code_llvm f3()
@code_llvm f3(cst=BigleafConstants())
@code_warntype f3()
@code_warntype f3(cst=BigleafConstants())
@code_warntype bigleaf_constants.g


text:
include("../src/hiscatter.jl")

using Distributions
using Plots, RecipesBase

data = rand(Pareto(2),100)

hiscatter(data,10)
hiscatter(data,20)

ext: jl


text: largez03() = largez(3)


text:
######################################################################
# Types UTs
# -----
# Licensed under MIT License
module TestTypes
using Test
using KirUtil
import ..Helpers: Immutable

@testset "KirUtil Types" begin
    @testset "Ident" begin
        identtest(::Ident{:foo}) = 42
        identtest(::Ident{:bar}) = 69.69
        identtest(ident::Symbol) = identtest(Ident{ident}())

        @test identtest(:foo) == 42
        @test identtest(:bar) == 69.69
    end

    @testset "Mutable" begin
        let immutable = Immutable(42, false)
            immutable.mutable[] = 42
            @test immutable.mutable[] == 42
        end
    end

    @testset "Optional" begin
        @testset "explicit" begin
            struct WithExplicitOptional{T}
                opt::Optional{:explicit, T}
            end
            WithExplicitOptional{T}() where T = WithExplicitOptional{T}(unknown)

            @test WithExplicitOptional{Real}().opt[] === unknown
            @test WithExplicitOptional{Real}(Optional()).opt[] === unknown
            @test WithExplicitOptional{Integer}(24).opt[] == 24
            @test WithExplicitOptional{Integer}(42.).opt[] == 42
            @test WithExplicitOptional{Integer}(Optional(42.)).opt[] == 42
            let obj = WithExplicitOptional{Real}(24)
                obj.opt[] = 69.69
                @test obj.opt[] == 69.69
            end
            let obj = WithExplicitOptional{Integer}()
                obj.opt[] = 420
                @test obj.opt[] == 420
            end
        end

        @testset "generic" begin
            struct WithGenericOptional
                opt::Optional{:generic}
            end
            WithGenericOptional() = WithGenericOptional(unknown)

            @test WithGenericOptional().opt[] === unknown
            @test WithGenericOptional(24).opt[] == 24
            @test WithGenericOptional(Optional(24)).opt[] == 24
            let obj = WithGenericOptional()
                obj.opt[] = 24
                @test obj.opt[] == 24
            end
        end

        @testset "symbolless" begin
            struct WithSymbollessOptional{T}
                opt::Optional{S, T} where S
            end
            WithSymbollessOptional{T}() where T = WithSymbollessOptional{T}(unknown)

            @test WithSymbollessOptional{Real}().opt[] === unknown
            @test WithSymbollessOptional{Real}(Optional(42)).opt[] == 42
            @test WithSymbollessOptional{Integer}(Optional(42.)).opt[] == 42
        end

        @testset "implicit load" begin
            KirUtil.load(opt::Optional{:implicit_load}) = opt.value = 42
            let opt = Optional{:implicit_load}()
                @assert isunknown(opt)
                @test opt[] == 42
                @test !isunknown(opt)
            end
        end

        @testset "explicit load" begin
            KirUtil.load(io::IO, opt::Optional{:explicit_load, T}) where T = opt.value = read(io, T)
            let buff = IOBuffer(), opt = Optional{:explicit_load, Int}()
                write(buff, 42)
                seek(buff, 0)
                @test load(buff, opt) == 42
                @test opt[] == 42
            end
            let buff = IOBuffer(), opt = Optional{:explicit_load, Float32}()
                write(buff, 24.f0)
                seek(buff, 0)
                @test load(buff, opt) ≈ 24
                @test opt[] ≈ 24
            end
        end
    end

    @testset "Dirty" begin
        let dirty = Dirty(42)
            @test !isdirty(dirty)
            dirty[] += 1
            @test isdirty(dirty)
            @test dirty[] == 43
        end
    end
end
end # module TestTypes


text:
Aug16 = Dates.days(Date(2016, 7, 27)-today())/365
Sep16 = Dates.days(Date(2016, 8, 26)-today())/365
Oct16 = Dates.days(Date(2016, 9, 26)-today())/365
Nov16 = Dates.days(Date(2016, 10, 27)-today())/365
Dec16 = Dates.days(Date(2016, 11, 25)-today())/365
Jan17 = Dates.days(Date(2016, 12, 23)-today())/365
Feb17 = Dates.days(Date(2017, 1, 27)-today())/365
Mar17 = Dates.days(Date(2017, 2, 24)-today())/365
Apr17 = Dates.days(Date(2017, 3, 27)-today())/365
May17 = Dates.days(Date(2017, 4, 26)-today())/365
Jun17 = Dates.days(Date(2017, 5, 26)-today())/365
Jul17 = Dates.days(Date(2017, 6, 26)-today())/365
Aug17 = Dates.days(Date(2017, 7, 27)-today())/365
Sep17 = Dates.days(Date(2017, 8, 25)-today())/365
Oct17 = Dates.days(Date(2017, 9, 26)-today())/365
Nov17 = Dates.days(Date(2017, 10, 27)-today())/365
Dec17 = Dates.days(Date(2017, 11, 24)-today())/365
Jan18 = Dates.days(Date(2017, 12, 27)-today())/365
Feb18 = Dates.days(Date(2018, 1, 26)-today())/365
Mar18 = Dates.days(Date(2018, 2, 23)-today())/365
Apr18 = Dates.days(Date(2018, 3, 27)-today())/365
May18 = Dates.days(Date(2018, 4, 26)-today())/365
Jun18 = Dates.days(Date(2018, 5, 25)-today())/365
Jul18 = Dates.days(Date(2018, 6, 26)-today())/365
Aug18 = Dates.days(Date(2018, 7, 27)-today())/365
Sep18 = Dates.days(Date(2018, 8, 27)-today())/365
Oct18 = Dates.days(Date(2018, 9, 26)-today())/365
Nov18 = Dates.days(Date(2018, 10, 26)-today())/365
Dec18 = Dates.days(Date(2018, 11, 26)-today())/365

Q416 = [Oct16, Nov16, Dec16]
Q117 = [Jan17, Feb17, Mar17]
Q217 = [Apr17, May17, Jun17]
Q317 = [Jul17, Aug17, Sep17]
Q417 = [Oct17, Nov17, Dec17]
Q118 = [Jan18, Feb18, Mar18]
Q218 = [Apr18, May18, Jun18]
Q318 = [Jul18, Aug18, Sep18]
Q418 = [Oct18, Nov18, Dec18]

W16 = [Q416, Q117]
S17 = [Q217, Q317]
W17 = [Q417, Q118]
S18 = [Q218, Q318]

# Date(2017,12) - Dates.Day(7)
# Dates.dayofweek(Date(2017,10) - Dates.Day(5))
# Date(2017,7)
#
# Dates.days(Date(2016,7) - Dates.Day(5)
#
# if (Dates.dayofweek(Date(2016,7) - Dates.Day(5)) == 7)


text:
using FinancialToolbox, HyperDualNumbers, DualNumbers, FinancialFFT, FinancialMonteCarlo

A = 600.0;
N = 20000;

method = LewisMethod(A, N);
# S0 = 100.0;
S0 = hyper(100.0, 1.0, 1.0, 0.0);
K = 100.0;
# r = hyper(0.02, 1.0, 1.0, 0.0);
r = 0.02;
T = 1.0;
d = 0.0;
#sigma=dual(0.2,1.0);
sigma = 0.2;
# sigma = hyper(0.2, 1.0, 1.0, 0.0);
mu1 = 0.03;
sigma1 = 0.02;
p = 0.3;
lam = 5.0;
lamp = 30.0;
lamm = 20.0;
zero_rate = ZeroRate(r);

Model = MertonProcess(sigma, lam, mu1, sigma1, Underlying(S0, d));
# Model = KouProcess(sigma, lam, p, lamp, lamm, Underlying(S0, d));
#Model = BlackScholesProcess(sigma, Underlying(S0, d));
# Model = VarianceGammaProcess(sigma, 0.01, 0.01, Underlying(S0, d));
#Model = NormalInverseGaussianProcess(sigma, 0.01, 0.01, Underlying(S0, d));

EUData = EuropeanOption(T, K)
Nsim = 10_000;
Nstep = 30;
mc = MonteCarloConfiguration(Nsim, Nstep);

@show pricer(Model, zero_rate, method, EUData);
@show pricer(Model, zero_rate, mc, EUData);
#typeof(Model) <: BlackScholesProcess ? @time(@show(blsprice(S0, K, r, T, sigma, d))) : nothing;


text:
#####
##### Benchmark from Pascal Pons
#####

# A benchmark to evaluate connect-four agents is available at:
#     http://blog.gamesolver.org/solving-connect-four/02-test-protocol/
# It features 6000 positions along with their expected scores.

# Let `ne` the number of elapsed moves. A game stage is either
# `begin` (ne <= 14), `middle` (14 < ne <= 28) or `end` (ne > 28)
const STAGES = [:begin, :middle, :end]

# Let `nr` the number of remaining moves. A difficulty level is either
# `easy` (nr < 14), `medium` (14 <= nr < 28) or `hard` (nr > 28).
const DIFFICULTIES = [:easy, :medium, :hard]

# Return either `nothing` or a `(stage, level)` tuple
function parse_test_filename(name)
  try
    m = match(r"^Test_L(\d)_R(\d)$", name)
    L = parse(Int, m.captures[1])
    R = parse(Int, m.captures[2])
    return (STAGES[L], DIFFICULTIES[R])
  catch e
    rethrow(e)
    return nothing
  end
end

struct Benchmark
  stage :: Symbol
  difficulty :: Symbol
  entries :: Vector{Tuple{String, Int}}
end

function load_benchmarks(dir)
  @assert isdir(dir)
  benchmarks = Benchmark[]
  files = readdir(dir)
  for bf in files
    meta = parse_test_filename(bf)
    f = joinpath(dir, bf)
    if !isnothing(meta) && isfile(f)
      stage, difficulty = meta
      entries = []
      for L in readlines(f)
        try
          L = split(L)
          push!(entries, (L[1], parse(Int, L[2])))
        catch e
          rethrow(e)
        end
      end
    push!(benchmarks, Benchmark(stage, difficulty, entries))
    end
  end
  return benchmarks
end

const BENCHMARKS_DIR = joinpath(@__DIR__, "..", "benchmark")
const BENCHMARKS = load_benchmarks(BENCHMARKS_DIR)

####
#### Load a session and run the test
####

using AlphaZero
include("../main.jl")
using .ConnectFour
import .ConnectFour: Solver

using ProgressMeter

function state_of_string(str)
  g = Game()
  for c in str
    a = parse(Int, c)
    GI.play!(g, a)
  end
  return g
end

function optimal_on(solver, player, e)
  s = state_of_string(e[1])
  as, π = AlphaZero.think(player, s)
  api = argmax(π)
  qs = [Solver.qvalue(solver, s, a) for a in as]
  return sign(qs[api]) == maximum(sign, qs)
end

function test_player_on(player, benchmark, progress)
  solver = Solver.Player()
  won = 0
  for e in benchmark.entries
    if optimal_on(solver, player, e)
      won += 1
    end
    next!(progress)
  end
  return won / length(benchmark.entries)
end

function test_player(player)
  for bench in BENCHMARKS
    p = Progress(length(bench.entries))
    score = test_player_on(player, bench, p)
    println("($(bench.stage), $(bench.difficulty)): $score")
  end
end

player = MinMax.Player{Game}(depth=5, τ=0)
#player = AlphaZero.RandomPlayer{Game}()
test_player(player)

#= Scores for minmax (depth 5)
(begin, easy):     0.983
(begin, medium):   0.927
(begin, hard):     0.766
(middle, easy):    0.993
(middle, medium):  0.908
(end, easy):       0.986
=#


text:
@nullable mutable struct BattleContext
    battle::Battle = Battle()
    parsed_battle::ParsedBattle
    data::VsContextData = VsContextData()
end

@nullable mutable struct ParsingContext
    battles::Vector{Battle} = Battle[]
    parsed_battles::Vector{ParsedBattle} = ParsedBattle[]
    current::BattleContext
end


text:
#!/usr/bin/julia

using DBI
using MySQL
using Dates

# Needs to getenv

host = get(ENV, "DB_HOST", "localhost")
user = get(ENV, "DB_USER", "")
passwd = get(ENV, "DB_PASSWD", "")
dbname = get(ENV, "DB_NAME", "")
port = parse(get(ENV, "DB_PORT", "3306"))
unix_socket = get(ENV, "DB_UNIX_SOCKET", C_NULL)
client_flag = 0

db = connect(MySQL5, host, user, passwd, dbname, port, unix_socket, client_flag)

# The analysis is very limited for the time being:
# The distribution of times required for departing from one station and
# arriving at the next station.
#
# The process is performed by looking at the schedule for all train numbers
#

type Stop
  station::String
  arrival::DateTime
  departure::DateTime
end

type TrainRun
  TrainNum::Int
  schedule::Tuple{Stop}
  actual::Tuple{Stop}
end

# Generates a schedule object
function read_schedules()
end

# Start by getting the list of trains
get_trains_query = "SELECT TrainNum, OrigTime FROM trains;"

get_trains = prepare(db, get_trains_query)

execute(get_trains)

disconnect(db)


text:
# Commutative semirings have commutative multiplication (addition is commutative by definition).
@traitdef IsCommutative{S<:AbstractSemiringElement}
@traitimpl IsCommutative{S} <- iscommutative(S)
is_commutative(x) = false

# Idempotent semirings are semirings with idempotent addition: `a+a = a`
@traitdef IsIdempotent{S<:AbstractSemiringElement}
@traitimpl IsIdempotent{S} <- isidempotent(S)
isidempotent(::Type{<:AbstractSemiringElement}) = false

# Star semirings have an additional unary operator `*`: `a^* = \sum_{i = 0}^{\infty} a^i`.
@traitdef IsStar{S<:AbstractSemiringElement}
@traitimpl IsStar{S} <- isstar(S)
isstar(::Type{<:AbstractSemiringElement}) = false


text:
using Aqua
using Ripserer

Aqua.test_all(Ripserer; ambiguities=false)


text:
"""

Set the color for drawing operations on `window`.

"""
setcolor!(window::Window, r::Int, g::Int, b::Int, a::Int) = SDL.SetRenderDrawColor(window.render_ptr, r, g, b, a)
setcolor!(window::Window, color::Colors.Color) = setcolor!(window, rgba_fromcolor(color)...)
setcolor!(window::Window, color::Colors.Color, a::Int) = setcolor!(window, rgba_fromcolor(color, a)...)

"""
    clear!(window)

Fill `window` with the currently selected color.

"""
clear!(window::Window) = SDL.RenderClear(window.render_ptr)

"""
    present!(window)

Update `window` to display all render operations since the last call
to `present!`.

"""
present!(window::Window) = SDL.RenderPresent(window.render_ptr)

##################################################
# Generic rendering
##################################################

function render!(layer::AbstractLayer, source::Union{AbstractGraphics,AbstractShape}, x, y, θ = 0.0; scale = 1.0, offset_x = 0, offset_y = 0, flip::Symbol = :none, color = colorant"black")
    push!(layer.render_tasks, RenderTask(source, x, y, θ, scale, offset_x, offset_y, flip, color))
    return layer
end

function render!(layer::AbstractLayer, text::Text, x, y, θ = 0.0; scale = 1.0, offset_x = 0, offset_y = 0, flip::Symbol = :none, color = colorant"black")
    push!(layer.render_tasks, RenderTask(text.graphics, x, y, θ, scale, offset_x, offset_y, flip, color))
    return layer
end

##################################################
# AbstractGraphics
##################################################

function render!(window::Window, layer::AbstractLayer, r::RenderTask{<:AbstractGraphics})
    x = r.x*layer.scale + layer.x
    y = r.y*layer.scale + layer.y
    scale = layer.scale*r.scale
    render!(window, r.source, round(Int, x), round(Int, y), r.θ, scale, scale, r.offset_x, r.offset_y, r.flip)
end

function render!(window::Window, source::Texture{SDL.Texture}, x::Int, y::Int, θ::Float64 = 0.0, scale_x::Float64 = 1.0, scale_y::Float64 = 1.0, offset_x::Int = 0, offset_y::Int = 0, flip::Symbol = :none)
    center_x, center_y = (source.center_x + offset_x)*scale_x, (source.center_y + offset_y)*scale_y
    rect = SDL.Rect(round(Int, x - center_x), round(Int, y - center_y),
                    round(Int, source.width*scale_x), round(Int, source.height*scale_y))
    point = SDL.Point(round(Int, center_x), round(Int, center_y))
    flip = UInt32(flip === :horizontal ? SDL.FLIP_HORIZONTAL :
                  flip === :vertical ? SDL.FLIP_VERTICAL :
                  SDL.FLIP_NONE)
    GC.@preserve rect point SDL.RenderCopyEx(window.render_ptr, source.ptr,
                                             C_NULL, pointer_from_objref(rect),
                                             θ, pointer_from_objref(point),
                                             flip)
    return window
end

##################################################
# Shapes
##################################################

"""
"""
function drawrect!(window::Window, x::Int, y::Int, w::Int, h::Int)
    rect = SDL.Rect(x, y, w, h)
    GC.@preserve rect SDL.RenderDrawRect(window.render_ptr, pointer_from_objref(rect))
    return window
end

"""
"""
function fillrect!(window::Window, x::Int, y::Int, w::Int, h::Int)
    rect = SDL.Rect(x, y, w, h)
    GC.@preserve rect SDL.RenderFillRect(window.render_ptr, pointer_from_objref(rect))
    return window
end

"""
"""
function drawpoint!(window::Window, x::Int, y::Int)
    SDL.RenderDrawPoint(window.render_ptr, Int32(x), Int32(y))
    return window
end

"""
"""
function drawline!(window::Window, x1::Int, y1::Int, x2::Int, y2::Int)
    SDL.RenderDrawLine(window.render_ptr, Int32(x1), Int32(y1), Int32(x2), Int32(y2))
    return window
end

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

function render!(window::Window, p::Point{2}, x, y, θ = 0.0; scale = 1.0, color)
    setcolor!(window, color)
    drawpoint!(window, round(Int, x + first(p).x), round(Int, y + first(p).y))
end

function render!(window::Window, l::AbstractLineShape{N,2}, x, y, θ = 0.0; scale = 1.0, color) where N
    setcolor!(window, color)
    for e in edges(l |> rotate(θ) |> translate(x, y))
        @inbounds x1 = e[1].x*scale
        @inbounds x2 = e[2].x*scale
        @inbounds y1 = e[1].y*scale
        @inbounds y2 = e[2].y*scale
        drawline!(window, round(Int, x1), round(Int, y1), round(Int, x2), round(Int, y2))
    end
end

function render!(window::Window, l::HalfPlane, x, y, θ = 0.0; scale = 1.0, color) where N
    setcolor!(window, color)
    for e in edges(l |> rotate(θ) |> translate(x, y))
        @inbounds x1 = e[1].x*scale
        @inbounds x2 = e[2].x*scale
        @inbounds y1 = e[1].y*scale
        @inbounds y2 = e[2].y*scale
        drawline!(window, round(Int, x1), round(Int, y1), round(Int, x2), round(Int, y2))
    end
end

function render!(window::Window, layer::AbstractLayer, r::RenderTask{<:AbstractShape})
    x = r.x*layer.scale + layer.x
    y = r.y*layer.scale + layer.y
    scale = layer.scale*r.scale
    render!(window, r.source, x, y, r.θ, scale=scale, color=r.color)
end


text:
using Optim

mutable struct HLT
	α::Float64
	β::Float64
	l₀::Float64
	b₀::Float64

	HLT() = new(0, 0, 0, 0)
	HLT(α::Number, β::Number, l₀::Number, b₀::Number) = new(Float64(α), Float64(β),
		Float64(l₀), Float64(b₀))
end

function loss(model::HLT, time_series)
	α, β, l₀, b₀ = model.α, model.β, model.l₀, model.b₀

	N = length(time_series)

	l_t, l_t_, b_t = 0, 0, 0 # l_t_ is the variable to save l(t-1)
	loss = 0

	for t in 1:N
		if t == 1
			l_t = l₀
			b_t = b₀
		else
			l_t = time_series[t - 1] * α + (l_t + b_t) * (1 - α) #b_t is taking b(t-1) value
		end

		l_t_ = l_t
		y_pred = l_t + b_t
		loss += (time_series[t] - y_pred)^2
	end

	return loss
end

function fit(model::HLT, y)
	lower = [-Inf, -0.001, -Inf, -Inf]
	upper = [1., 1., Inf, Inf]
	initial = [model.α, model.β, model.l₀, model.b₀]

	function loss_(parameters::Array{Float64, 1})
		α, β, l₀, b₀ = parameters
		return loss(HLT(α, β, l₀, b₀), y)
	end
	res = Optim.optimize(loss_, lower, upper, initial)
	optimal = Optim.minimizer(res)
	return HLT(optimal[1], optimal[2], optimal[3], optimal[4])
end

function forecast(model::HLT, time_series, forecast_length)
	N = length(time_series)
	α, β, l₀, b₀ = model.α, model.β, model.l₀, model.b₀
	l_t, l_t_, b_t = 0, 0, 0
	pred = Array{Float64, 1}(undef, forecast_length)

    #go through the whole time series making the point by point estimate
	for t in 1:N
		if t == 1
			l_t = l₀
			b_t = b₀
		else
			l_t = time_series[t - 1] * α + (l_t + b_t) * (1 - α) #b_t "is" b(t-1)
			b_t = β * (l_t - l_t_) + (1 - β) * b_t
		end
		l_t_ = l_t
	end

    #The parameter´s values to make the forecast are those estimated in the last step of the time series
	l_t = time_series[end] * α + (l_t + b_t) * (1 - α)
	b_t = β * (l_t - l_t_) + (1 - β) * b_t

	for i in 1:forecast_length
		#y_pred = l_t + b_t * i
		#push!(pred, y_pred)
		pred[i] = l_t + b_t * i
	end

	return pred
end


text:
# Optimized math and matrix types
include("linalg/main.jl")

# Statistics for the DQMC simulation runtime. This has nothin to do with
# measurements/observables.
include("statistics.jl")

# Some tapes, checkerboard stuff (which has been neglected a lot)
include("abstract.jl")

# contains `DQMCParameters` which holds parameters relevant to DQMC and the
# general simulation.
include("parameters.jl")

# There are a bunch of functions that require the DQMC type, but logically fit
# with one of its constituents. For example `initialize_stack` requires
# a ::DQMC to figure out how large a bunch of matrices need to be, but
# logically fits in `stack.jl`.
mutable struct DQMC{
        M <: Model, CB <: Checkerboard, FT <: AbstractField, RT <: AbstractRecorder,
        Stack <: AbstractDQMCStack, UTStack <: AbstractDQMCStack,
        US <: AbstractUpdateScheduler
    } <: MonteCarloFlavor

    model::M
    field::FT
    last_sweep::Int

    stack::Stack
    ut_stack::UTStack
    scheduler::US
    parameters::DQMCParameters
    analysis::DQMCAnalysis

    recorder::RT
    thermalization_measurements::Dict{Symbol, AbstractMeasurement}
    measurements::Dict{Symbol, AbstractMeasurement}
    lattice_iterator_cache::LatticeIteratorCache

    function DQMC{M, CB, FT, RT, Stack, UTStack, US}(args...) where {
            M <: Model, CB <: Checkerboard, FT <: AbstractField,
            RT <: AbstractRecorder,
            Stack <: AbstractDQMCStack, UTStack <: AbstractDQMCStack,
            US <: AbstractUpdateScheduler
        }

        @assert isconcretetype(M)
        @assert isconcretetype(FT)
        @assert isconcretetype(Stack)
        @assert isconcretetype(UTStack)
        @assert isconcretetype(RT)
        @assert isconcretetype(US)

        new{M, CB, FT, RT, Stack, UTStack, US}(args..., LatticeIteratorCache())
    end
end

# Simplified constructor
function DQMC(
        CB, model::M, field::FT, last_sweep,
        stack::Stack, ut_stack::UTStack, scheduler::US,
        parameters, analysis,
        recorder::RT,
        thermalization_measurements, measurements
    ) where {M, FT, RT, Stack, UTStack, US}

    DQMC{M, CB, FT, RT, Stack, UTStack, US}(
        model, field, last_sweep, stack, ut_stack,
        scheduler, parameters, analysis, recorder,
        thermalization_measurements, measurements
    )
end

# copy constructor
function DQMC(
        mc::DQMC{x, CBT};
        CB = CBT, model::M = mc.model, field::FT = mc.field,
        last_sweep = mc.last_sweep,
        stack::Stack = mc.stack, ut_stack::UTStack = mc.ut_stack,
        scheduler::US = mc.scheduler, parameters = mc.parameters,
        analysis = mc.analysis, recorder::RT = mc.recorder,
        thermalization_measurements = mc.thermalization_measurements,
        measurements = mc.measurements
    ) where {x, CBT, M, FT, RT, Stack, UTStack, US}

    DQMC{M, CB, FT, RT, Stack, UTStack, US}(
        model, field, last_sweep, stack, ut_stack,
        scheduler, parameters, analysis, recorder,
        thermalization_measurements, measurements
    )
end

# Contains different auxiliary field types
include("fields.jl")

# Contains `DQMCStack`, `propagate` and related functions. The stack keeps track
# of all the matrices that go into calculating the path integral. `propagate`
# moves imaginary time along, between `l = 1` and `l = slices`.
# This also contains some functions for calculating greens functions.
include("stack.jl")

# Contains the `UnequalTimeStack`, related greens calculation methods and
# iterators (i.e. CombinedGreens). This is used for unequal time measurements.
include("unequal_time_stack.jl")

# Contains functions for computations using the matrices representing a time
# slice. Also contains special methods for Checkerboard decomposed matrices.
include("slice_matrices.jl")

# Contains global update schedulers, the updates themselves as well as some
# functions necessary to calculate and perform those updates.
# include("global.jl")
include("updates/scheduler.jl")
include("updates/local_updates.jl")
include("updates/global_updates.jl")
include("updates/parallel_updates.jl")

# Contains the "main/top-level" simulation code. Specifically construction and
# initialization of the DQMC type, convenience functions, `run!` and `replay!`,
# and local updates.
include("DQMC.jl")

# Contains mandatory and optional method definitions expected to exist for a
# given model.
include("DQMC_interface.jl")

# saving and loading of DQMC
include("FileIO.jl")

# Contains some greens methods
# some others are in stack.jl and unequal_time_stack.jl
include("greens.jl")

# Contains code related to make measurements. Specifically:
# The overall structure
include("measurements/generic.jl")
# Quick constructers and measurement kernels (i.e. applied Wicks theorem)
include("measurements/measurements.jl")
# Superfluid stiffness stuff
include("measurements/superfluid_stiffness.jl")
# Contains some post processing tools
include("measurements/extensions.jl")
# structs and conversions from the old system
include("measurements/deprecated.jl")


text:
using Random
using Gadfly
using LinearAlgebra
using DataFrames

hinge_loss(x :: Float64 ) = max(0,1-x)

mutable struct HyperPlan{d}

    w :: Vector{Float64}
    b :: Float64

    function HyperPlan{d}() where {d}
        new(zeros(d), 0.0)
    end

    function HyperPlan{d}(w :: Vector{Float64}) where {d}
        n = length(w)
        if d >= n
            w = [w;zeros(d-n)]
        else
            w = w[1:d]
        end
        new(w, 0.0)
    end

    function HyperPlan{d}(w :: Vector{Float64},
                          b :: Float64) where {d}
        n = length(w)
        d >= n && (w = [w;zeros(d-n)])
        d >= n || (w = w[1:d])
        new(w, b)
    end
end

Base.ndims(::HyperPlan{d}) where {d} = d

function (h :: HyperPlan)( x :: Vector{Float64})
    dot(h.w, x) + h.b
end

import Base:+

function +(h1 :: HyperPlan{d}, h2 :: HyperPlan{d}) where {d}

    HyperPlan{d}(h1.w + h2.w, h1.b + h2.b)

end

import Base:-

function -(h1 :: HyperPlan{d}, h2 :: HyperPlan{d}) where {d}

    HyperPlan{d}(h1.w - h2.w, h1.b - h2.b)

end

function external_prod(h :: HyperPlan{d},
                       a :: Float64) where {d}

    HyperPlan{d}(h.w .* a, h.b * a)

end

square(h :: HyperPlan) = dot(h.w, h.w) + h.b^2

function gradient_loss(h     :: HyperPlan{d},
                       x     :: Vector{Float64},
                       y     :: Float64,
                       lambd :: Float64) where {d}
    if y * h(x) < 1
        v1 = lambd .* h.w .- y .* x
        v2 = lambd * h.b - y
    else
        v1 = lambd .* h.w
        v2 = lambd * h.b
    end
    HyperPlan{d}(v1, v2)
end

function SGD(h     :: HyperPlan{d},
             x     :: Array{Float64,2},
             y     :: Vector{Float64},
             lambd :: Float64,
             niter     = 1000,
             learning  = 0.1,
             averaging = 1) where {d}

    i0       = abs(niter-averaging)
    liste    = HyperPlan{d}[]
    push!(liste, h)
    res0, res1 = h, h
    n = length(y)
    for i in 1:niter
        it = rand(1:n)
        vt = gradient_loss(res0, x[it,:], y[it], lambd)
        res0 = res0 - external_prod(vt,learning/(1+learning*lambd*i))
        push!(liste, res0)
        mut = 1/max(1,i-i0)
        res1 = (external_prod(res1, 1-mut)
              + external_prod(res0, mut))
        res0 = res1
    end
    res1, liste
end

using LinearAlgebra

mutable struct Polynomial

    weights :: Array{Float64, 2}
    void    :: Array{Float64, 2}

    function Polynomial(w)
        d = size(w)[1]
        void = zeros(Float64, (d,d))
        for i in 1:d
            for j in 1:d-i+1
                void[i,j] = 1
            end
        end
        weights = void .* w
        new( weights, void)
    end
end

degree(p :: Polynomial) = size(p.weights)[1]-1

function (p :: Polynomial)( x :: Vector{Float64} )
    @assert length(x) > 1
    X1 = [x[1]^exp for exp in 0:degree(p)]
    X2 = [x[2]^exp for exp in 0:degree(p)]
    ((p.weights * X1)' * X2)[1,1]
end

# Les paramètres :
# dimension de l'espace des x
h = HyperPlan{2}()
# degré du polynôme de départ
deg = 1
# on introduit ce petit pourcentage d'erreurs dans les données
seuil_erreur = 0.
# taille de l'échantillon de données d'entraînement
m=100
# taille de l'échantillon de données de test
mtest=1000
# bornes sur l'espace des données x
xmax =  1.25
xmin = -1.25
# valeur de lambda pour le soft-SVM
lambd=1.
# valeur du taux d'apprentissage initialement pour le soft-SVM
learning  = 0.1
# taille du stencil pour l'averaging dans la SGD
averaging = 10

# fabrication des données
Random.seed!(1234)
w   = reshape(collect(0:(deg+1)*(deg+1)-1),deg+1,deg+1)
w[1,1] = 0.0

pf = Polynomial(w)
#x  = xmin .+ (xmax .- xmin) .* rand(Float64,(m, ndims(h)))
x1 = collect(range(xmin, stop=xmax, length=10))  .* ones(10)'
x2 = collect(range(xmin, stop=xmax, length=10))' .* ones(10)
x  = hcat(vec(x1),vec(x2))

u  = rand(Float64,m)
y  = zeros(m)
y[ u .>  seuil_erreur ] .=   sign.([pf(x[i,:]) for i in 1:m if u[i] >   seuil_erreur])
y[ u .<= seuil_erreur ] .= - sign.([pf(x[i,:]) for i in 1:m if u[i] <=  seuil_erreur])

plot(x = x[:,1], y = x[:,2], color = y)

Random.seed!(5678)
xtest = xmin .+ (xmax - xmin) .* rand(Float64,( mtest, ndims(h)))
ytest = sign.([pf(xtest[i,:]) for i in 1:mtest])
plot(x = xtest[:,1], y = xtest[:,2], color = ytest)

# Résolution du problème soft-SVM

h, liste = SGD(h, x, y, lambd, 200, learning, averaging)

h

function affiche(x, y, h :: HyperPlan)
    m = length(y)
    yresult = [sign(h(x[i,:])) for i in 1:m]

    p1 = plot( x = x[:,1], y = x[:,2], color=y)
    p2 = plot( x = x[:,1], y = x[:,2], color=yresult)

    hstack(p1, p2)
end

affiche(x, y, h)

affiche(xtest,ytest,h)

function loss(h :: HyperPlan{d},
     x :: Vector{Float64},
     y :: Float64) where {d}

    hinge_loss(y * h(x))
end

function risk(h     :: HyperPlan{d},
              x     :: Array{Float64, d},
              y     :: Array{Float64, 1},
              lambd :: Float64) where {d}

    (lambd/2 * square(h)
        +1/length(y)*sum([loss(h, x[i,:], y[i])
                for i in 1:size(x)[1]]))
end

function risk_vrai(h :: HyperPlan{d},
                   x :: Array{Float64,d},
                   y :: Vector{Float64}) where {d}
    n = length(y)
    @assert n == size(x)[1]
    sum([sign.(h(x[i,:])) != y[i] for i in 1:n])/n
end

function affiche_risques(liste, x, y, lambd)
    y1 = [risk(h, x, y, lambd) for h in liste]
    y2 = [risk_vrai(h, x, y) for h in liste]
    df1 = DataFrame( x = collect(1:length(liste)),
        y = y1, label = :risque)
    df2 = DataFrame( x = collect(1:length(liste)),
        y = y2, label = :risque_vrai)
    df = vcat(df1,df2)
    plot(df, x=:x, y=:y, Geom.line,
        color= :label,
        Scale.color_discrete_manual("blue","red"))
end

affiche_risques(liste,xtest,ytest,lambd)


text:
## A small driver file to run the saturation tests and combine
#run like this
#julia -p 4 -L testSaturation_Harness.jl test_stauration.jl
#ARGS[1] is numRun per batch

using Distributed

const spath = "./Results/satExp"
const param_path = "./Results/param_portExp_mtn2.csv"

n_grid = [2^i for i = 5:17]
numRuns = parse(Int, ARGS[1])
alpha = parse(Float64, ARGS[2])

start_time = time_ns()
a = @spawn test_saturation(spath, numRuns, n_grid, 8675309, param_path, alpha=alpha)
b = @spawn test_saturation(spath, numRuns, n_grid, 5164290, param_path, alpha=alpha)
c = @spawn test_saturation(spath, numRuns, n_grid, 1236, param_path, alpha=alpha)
d = @spawn test_saturation(spath, numRuns, n_grid, 51666, param_path, alpha=alpha)

######
file_a = fetch(a)
file_b = fetch(b)
file_c = fetch(c)
file_d = fetch(d)

time_stamp = (time_ns() - start_time) * 1e-9

##read everyone in, throw away a line
data, header = readdlm(file_a, ',', header=true)

data_t = readdlm(file_b, ',', skipstart=1)

data_t[:, 1] .+= numRuns
data = vcat(data, data_t)

data_t = readdlm(file_c, ',', skipstart=1)
data_t[:, 1] .+= 2numRuns
data = vcat(data, data_t)

data_t = readdlm(file_d, ',', skipstart=1)
data_t[:, 1] .+= 3numRuns
data = vcat(data, data_t)

#strip the name of file_a to make the numbers better
println("Filea \t", file_a)
println("spath \t", spath)

indx =  findfirst(spath, file_a)[end] + 1

f = open("$(spath)_$(file_a[indx:end])_full_$(round(alpha, digits=6))_$(4numRuns).csv", "w")
writedlm(f,  header, ',')
writedlm(f,  data, ',')
close(f)

println("Num Paths: \t $(numRuns) \t Time:", time_stamp )

max_stars_repo_path: Normal Model/test_saturation.jl

max_issues_repo_path: Normal Model/test_saturation.jl

max_forks_repo_path: Normal Model/test_saturation.jl


text:
import MechanicalSketch
import MechanicalSketch: color_with_lumin, empty_figure, background, sethue, O, EM, FS, finish, PALETTE
import MechanicalSketch: dimension_aligned, settext, arrow, setfont
import MechanicalSketch: ComplexQuantity, generate_complex_potential_vortex, @import_expand, string_polar_form
import MechanicalSketch: place_image, x_y_iterators_at_pixels, ColSchemeNoMiddle
import MechanicalSketch: BinLegend, draw_legend, latexify, @layer, fontsize

let
if !@isdefined m²
    @import_expand(m, s)
end

BACKCOLOR = color_with_lumin(PALETTE[8], 80)
function restart()
    empty_figure(filename = joinpath(@__DIR__, "test_21.png"))
    background(BACKCOLOR)
    sethue(PALETTE[5])
end
restart()

"Vortex position"
p_vortex = (0.0, 2.5)m |> ComplexQuantity

"Vorticity, 2d flow"
K = 1.0m²/s

"2d velocity potential function. Complex quantity domain, real quantity range. "
ϕ_vortex = generate_complex_potential_vortex(; pos = p_vortex, vorticity = K)

physwidth = 20m
physheight = physwidth / 3
xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
A = [ϕ_vortex(complex(x, y)) for y in ys, x in xs]
legend = BinLegend(;maxlegend = maximum(A), minlegend = minimum(A), noofbins = 20,
    name = :Potential)

colormat = legend.(A);
ulp, lrp = place_image(O, colormat);
legendpos = O + ((lrp[1] - ulp[1]) / 2 + EM, (ulp[2] - lrp[2]) / 2)
@layer begin
    fontsize(0.8FS)
    draw_legend(legendpos, legend)
end

# Add some decoration to the plot
sethue(color_with_lumin(PALETTE[5], 30))
dimension_aligned(O + (-physwidth / 2, physheight / 2), O + (physwidth / 2, physheight / 2))
dimension_aligned(O + p_vortex, O)
dimension_aligned(O + (-physwidth / 2, - physheight / 2 ),  O +  (-physwidth / 2, physheight / 2 ))

@layer begin # Leader for a value
    sethue(color_with_lumin(PALETTE[5], 90))
    p0 = (0.0 + 1.0im)m
    arrow(O + p0 - (3EM, -2EM), O + p0 )
    fvalue = ϕ_vortex(p0)
    strvalue = string(round(typeof(fvalue ), fvalue , digits = 3))
    strargument = string(p0)
    strargumentpol = string_polar_form(round(typeof(p0 ), p0 , digits = 3))
    txt =  "ϕvortex( $strargument ) = ϕvortex( $strargumentpol ) =$strvalue"
    settext(txt, O + p0 + (-12EM, 3EM), markup=true)
end

# Convert the function to latex and display
la1 = latexify(quote
    ϕ_vortex(θ:: Angle)::Potential = $K * θ
end);
pt = O + (-14.0, -5.0)m
ptul, ptlr, scalefactor = place_image(pt, la1, height = 1EM)

# Following text
#setfont("DejaVu Sans", FS)
str = "is, surprisingly, irrotational outside of the centre"
settext(str, ptlr + (EM, 0), markup = true)

# Convert the outer function to latex
la2 = latexify(quote
    ϕ_vortex(p::ComplexQuantity)::Potential = ϕ_vortex(-angle(p - $p_vortex))
end)
place_image(pt + (0, 3EM), la2; scalefactor)

finish()
nothing
end # let


text:
_isequal(p, q) = all(JuMP.isequal_canonical.(coefficients(p), coefficients(q)))
function _isequal(x::AbstractArray, y::AbstractArray)
    size(x) == size(y) && all(_isequal.(x, y))
end

@testset "@constraint macro with polynomials" begin
    m = Model()
    setpolymodule!(m, TestPolyModule)
    @variable m α
    @variable m β
    @polyvar x y
    p = α * x*y + β * x^2
    #q = MatPolynomial([α β; β α], [x])
    q = α*x^2 + β*x*y + α*y^2
    @test_macro_throws ErrorException @constraint(m, p)
    @test_macro_throws ErrorException @constraint(m, begin p >= 0 end)
    @test_macro_throws ErrorException @constraint(m, +(p, p, p))
    @test_macro_throws ErrorException @constraint(m, p >= 0, 1)
    #@test_macro_throws ErrorException @constraint(m, p >= 0, domain = (@set x >= -1 && x <= 1, domain = y >= -1 && y <= 1))
    @test_macro_throws ErrorException @constraint(m, p + 0, domain = (@set x >= -1 && x <= 1))

    function testcon(m, cref, S::Type, jump_set::PolyJuMP.PolynomialSet,
                     p, ineqs, eqs, basis=MB.MonomialBasis, kwargs=[])
        @test cref isa JuMP.ConstraintRef{
            Model, <:MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},
                                         <:S}}
        c = JuMP.constraint_object(cref)
        set = JuMP.moi_set(c)
        @test set isa S
        if set isa PolyJuMP.PlusMinusSet
            set = set.set
        end
        expected_str = string(JuMP.function_string(REPLMode, p), ' ',
                              JuMP._math_symbol(REPLMode, :in), ' ', jump_set)
        @test sprint(show, MIME"text/plain"(), cref) == expected_str
        expected_str = string("\$ ", JuMP.function_string(IJuliaMode, p), ' ',
                              JuMP._math_symbol(IJuliaMode, :in), ' ', jump_set,
                              " \$")
        @test sprint(show, MIME"text/latex"(), cref) == expected_str
        # @test set.basis isa basis #breaks for TestPosDefMatrix
        if !isempty(kwargs)
            @test length(set.kwargs) == length(kwargs)
            for (i, kw) in enumerate(set.kwargs)
                @test kw.first == kwargs[i][1]
                @test kw.second == kwargs[i][2]
            end
        end
        # == between JuMP affine expression is not accurate, e.g. β + α != α + β
        # == 0 is not defined either
        # c.p and p can be matrices
        @test _isequal(JuMP.reshape_vector(JuMP.jump_function(c),
                                           JuMP.shape(c)), p)
        if isempty(ineqs)
            if isempty(eqs)
                @test isa(set.domain, FullSpace)
            else
                @test isa(set.domain, AlgebraicSet)
                @test equalities(set.domain) == eqs
            end
        else
            @test isa(set.domain, BasicSemialgebraicSet)
            @test inequalities(set.domain) == ineqs
            @test equalities(set.domain) == eqs
        end
    end

    f(x, y) = @set x + y == 2
    dom = @set x^2 + y^2 == 1 && x^3 + x*y^2 + y >= 1

    @testset "Printing" begin
        in_sym = JuMP._math_symbol(REPLMode, :in)
        eqref = @constraint(m, p == q)
        @test sprint(show, MIME"text/plain"(), eqref) == "(β - α)x² + (α - β)xy + (-α)y² $in_sym PolyJuMP.ZeroPoly()"
        @test sprint(show, MIME"text/latex"(), eqref) == "\$ (β - α)x^{2} + (α - β)xy + (-α)y^{2} \\in PolyJuMP.ZeroPoly() \$"
        sdref = @constraint(m, [p q; q p] in PSDCone())
        if VERSION < v"1.4-"
            @test sprint(show, MIME"text/plain"(), sdref) == "[(β)x² + (α)xy          (α)x² + (β)xy + (α)y²;\n (α)x² + (β)xy + (α)y²  (β)x² + (α)xy        ] $in_sym Main.TestPolyModule.TestPosDefMatrix()"
        else
            @test sprint(show, MIME"text/plain"(), sdref) == "[(β)x² + (α)xy          (α)x² + (β)xy + (α)y²;\n (α)x² + (β)xy + (α)y²  (β)x² + (α)xy] $in_sym Main.TestPolyModule.TestPosDefMatrix()"
        end
        @test sprint(show, MIME"text/latex"(), sdref) == "\$ \\begin{bmatrix}\n(β)x^{2} + (α)xy & (α)x^{2} + (β)xy + (α)y^{2}\\\\\n(α)x^{2} + (β)xy + (α)y^{2} & (β)x^{2} + (α)xy\\\\\n\\end{bmatrix} \\in Main.TestPolyModule.TestPosDefMatrix() \$"
    end

    @testset "NonNeg" begin
        S = TestPolyModule.NonNeg
        jump_set = TestPolyModule.TestNonNeg()
        testcon(m, @constraint(m, p >= q + 1, domain = @set y >= 1 && dom),
                S, jump_set, p - q - 1, [y - 1, x^3 + x*y^2 + y - 1],
                [x^2 + y^2 - 1])
        testcon(m, @constraint(m, p <= q),
                S, jump_set, -p + q, [], [])
        testcon(m, @constraint(m, q - p in PolyJuMP.NonNegPoly()),
                S, jump_set, q - p, [], [])
        testcon(m, @constraint(m, p + q >= 0, domain = @set x == y^3),
                S, jump_set, p + q, [], [x - y^3])
        @testset "Custom keyword" begin
            testcon(m, @constraint(m, p <= q, maxdegree=1),
                    S, jump_set, -p + q, [], [], MB.MonomialBasis,
                    [(:maxdegree, 1)])
        end
    end
    @testset "ZeroPolynomialSet" begin
        jump_set = PolyJuMP.ZeroPoly()
        @testset "ZeroPolynomialSet{FullSpace}" begin
            S = PolyJuMP.ZeroPolynomialSet{FullSpace}
            testcon(m, @constraint(m, p == q),
                    S, jump_set, p - q, [], [])
            testcon(m, @constraint(m, p - q in PolyJuMP.ZeroPoly()),
                    S, jump_set, p - q, [], [])
        end
        S = PolyJuMP.ZeroPolynomialSet
        testcon(m, @constraint(m, p == q, domain = @set x == 1 && f(x, y)),
                S, jump_set, p - q, [], [x + y - 2, x - 1])
        testcon(m, @constraint(m, p - q in PolyJuMP.ZeroPoly(), domain = @set x == 1 && f(x, y)),
                S, jump_set, p - q, [], [x + y - 2, x - 1])
        S = PolyJuMP.PlusMinusSet
        testcon(m, @constraint(m, p == q, domain = dom),
                S, jump_set, p - q, [x^3 + x*y^2 + y - 1],
                [x^2 + y^2 - 1])
    end
    @testset "PosDefMatrix" begin
        jump_set = TestPolyModule.TestPosDefMatrix()
        testcon(m, @SDconstraint(m, [p q; q 0] ⪰ [0 0; 0 p]),
                TestPolyModule.PosDefMatrix, jump_set, [p q; q -p], [], [])
    end
end


text:
module AllnixThinkJulia
include("ch15module.jl")
include("ch16module.jl")
include("ch17module.jl")
include("ch18module.jl")
end


text:
α = 0.025

@testset "helpers" begin
  Σ = [1 0.5 0.5 0.6; 0.5 1 0.5 0.6; 0.5 0.5 1. 0.6; 0.6 0.6 0.6 1.]
  Random.seed!(43)
  x = Array(transpose(rand(MvNormal(Σ),500000)))
  y = norm2unifind(x, [1,2], "frechet")
  @test pvalue(ExactOneSampleKSTest(y[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(y[:,2], Uniform(0,1))) > α
  @test corspearman(y) ≈ [1. 0.; 0. 1.] atol=1.0e-3
  y = norm2unifind(x, [1,2])
  @test pvalue(ExactOneSampleKSTest(y[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(y[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(y[:,3], Uniform(0,1))) > α
  @test corspearman(y) ≈ [1. 0. 0.; 0. 1. 0.; 0. 0. 1.] atol=1.0e-3
  s= [1. 0.4 0.3; 0.4 1. 0.2; 0.3 0.2 1.]
  @test meanΣ(s) == 0.3
  @test mean_outer(s, [[1,2]]) == 0.25
  p = parameters(s, [[1,2]])
  @test p == ([0.4], 0.25)
  @test are_parameters_good(p...) == true
  @test Σ_theor([0.5], 0.3, [[1,2], [4]]) == [1.0  0.5  0.3; 0.5  1.0  0.3; 0.3  0.3  1.0]
  Random.seed!(42)
  @test frechet(0.5, [1. 0.2; 0.4 .6]) == [1.0  0.2; 0.6 0.6]
  Random.seed!(42)
  x = rand(1000, 5)
  c = getcors_advanced(x)
  @test c[1] == [[1, 4], [2, 3, 5]]
  @test c[2] ≈ [0.04729, 0.0195683] atol=1.0e-4
  @test c[3] ≈ -0.021774209774209772
  x = frechet(0.6, rand(100000, 4))
  Σ = cor(x)
  @test Σ[1,2] ≈ 0.6 atol=1.0e-2
  @test Σ[3,2] ≈ 0.6 atol=1.0e-2
  @test tail(x[:,1], x[:,2], "r") ≈ 0.6 atol=1.0e-1
  @test tail(x[:,1], x[:,2], "l") ≈ 0.6 atol=1.0e-1
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
end

@testset "convert sub-copula to archimedean" begin
  Random.seed!(43)
  Σ = cormatgen(25)
  Σ1 =0.8*ones(25,25) + 0.2*Matrix(1.0I, 25, 25)
  Σ = 0.3*Σ + 0.7*Σ1
  S = rand([0.8, 0.9, 1, 1.1, 1.2], 25)
  y = rand(MvNormal(Σ), 100000)'
  y = y.*S'
  d=["clayton" => [1,2,3,4,9]]
  x = gcop2arch(y, d)
  @test pvalue(ExactOneSampleKSTest(x[:,1], Normal(0,S[1]))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Normal(0,S[3]))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,4], Normal(0,S[4]))) > α
  @test norm(cor(y)-cor(x))/norm(cor(y)) < 0.045
  @test norm(cov(y)-cov(x))/norm(cov(y)) < 0.045
  @test maximum(abs.(cor(y)-cor(x))) < 0.11
  x1 = gcop2arch(y, d; notnested = true)
  @test pvalue(ExactOneSampleKSTest(x1[:,1], Normal(0,S[1]))) > α
  @test pvalue(ExactOneSampleKSTest(x1[:,3], Normal(0,S[3]))) > α
  @test pvalue(ExactOneSampleKSTest(x1[:,4], Normal(0,S[4]))) > α
  @test norm(cor(y)-cor(x1))/norm(cor(y)) < 0.045
  @test norm(cov(y)-cov(x1))/norm(cov(y)) < 0.045
  @test maximum(abs.(cor(y)-cor(x1))) < 0.11
  x2 = gcop2arch(y, d; naive = true)
  @test pvalue(ExactOneSampleKSTest(x2[:,1], Normal(0,S[1]))) > α
  @test pvalue(ExactOneSampleKSTest(x2[:,3], Normal(0,S[3]))) > α
  @test pvalue(ExactOneSampleKSTest(x2[:,4], Normal(0,S[4]))) > α
  @test maximum(abs.(cov(y[:,1:4])-cov(x2[:,1:4]))) < 0.05
  @test_throws AssertionError gcop2arch(y, ["clayton" => [1,1,3,4]])
end

@testset "convert sub-copula to t-Student" begin
  Random.seed!(42)
  Σ = cormatgen(25)
  S = rand([0.8, 0.9, 1, 1.1, 1.2], 25)
  mu = rand([0.8, 0.9, 1, 1.1, 1.2], 25)
  y = rand(MvNormal(Σ), 100000)'
  y = y.*S'.+mu'
  x = gcop2tstudent(y, [1,2,3,4], 10)
  @test pvalue(ExactOneSampleKSTest(x[:,1], Normal(mu[1],S[1]))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Normal(mu[3],S[3]))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,4], Normal(mu[4],S[4]))) > α
  @test norm(cor(y)-cor(x))/norm(cor(y)) < 0.015
  @test norm(cov(y)-cov(x))/norm(cov(y)) < 0.015
  @test maximum(abs.(cor(y)-cor(x))) < 0.02
  @test_throws AssertionError gcop2tstudent(y, [1,1,3,4], 10)
  x2 = gcop2tstudent(y, [1,2,3,4], 10; naive = true)
  @test pvalue(ExactOneSampleKSTest(x2[:,3], Normal(mu[3],S[3]))) > α
  @test pvalue(ExactOneSampleKSTest(x2[:,4], Normal(mu[4],S[4]))) > α
  @test maximum(abs.(cov(y[:,1:4])-cov(x2[:,1:4]))) < 0.0055
end

@testset "convert sub-copula to Frechet" begin
  Random.seed!(42)
  Σ = cormatgen(25)
  S = rand([0.8, 0.9, 1, 1.1, 1.2], 25)
  mu = rand([0.8, 0.9, 1, 1.1, 1.2], 25)
  y = rand(MvNormal(Σ), 100000)'
  y = y.*S'.+mu'
  x = gcop2frechet(y, [1,2,3])
  @test pvalue(ExactOneSampleKSTest(x[:,1], Normal(mu[1],S[1]))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Normal(mu[2],S[2]))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Normal(mu[3],S[3]))) > α
  @test norm(cor(y)-cor(x))/norm(cor(y)) < 0.15
  @test norm(cov(y)-cov(x))/norm(cov(y)) < 0.15
  @test maximum(abs.(cor(y)-cor(x))) < 0.25
  @test_throws AssertionError gcop2frechet(y, [1,1,3,4])
  Random.seed!(42)
  x2 = gcop2frechet(y, [1,2,3,4]; naive = true)
  @test pvalue(ExactOneSampleKSTest(x2[:,1], Normal(mu[1],S[1]))) > α
  @test pvalue(ExactOneSampleKSTest(x2[:,2], Normal(mu[2],S[2]))) > α
  @test pvalue(ExactOneSampleKSTest(x2[:,3], Normal(mu[3],S[3]))) > α
  @test maximum(abs.(cov(y[:,1:4])-cov(x2[:,1:4]))) < 0.075
end

@testset "convert sub-copula to Marshall-Olkin" begin
  Random.seed!(42)
  Σ = cormatgen(25)
  S = rand([0.8, 0.9, 1, 1.1, 1.2], 25)
  mu = rand([0.8, 0.9, 1, 1.1, 1.2], 25)
  y = rand(MvNormal(Σ), 100000)'
  y = y.*S'.+mu'
  x = gcop2marshallolkin(y, [1,2])
  @test pvalue(ExactOneSampleKSTest(x[:,1], Normal(mu[1],S[1]))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Normal(mu[2],S[2]))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Normal(mu[3],S[3]))) > α
  @test maximum(abs.(cov(y[:,1:2])-cov(x[:,1:2]))) < 0.03
  @test maximum(abs.(cor(y[:,1:2])-cor(x[:,1:2]))) < 0.03
  @test_throws AssertionError gcop2marshallolkin(y, [1,1])
  @test_throws DomainError gcop2marshallolkin(y, [1,2], 1., -1.)
  @test_throws DomainError gcop2marshallolkin(y, [1,2], -1., 1.)
  @test_throws AssertionError gcop2marshallolkin(y, [1,3,4])
  Random.seed!(42)
  x2 = gcop2marshallolkin(y, [1,2]; naive = true)
  @test pvalue(ExactOneSampleKSTest(x2[:,1], Normal(mu[1],S[1]))) > α
  @test pvalue(ExactOneSampleKSTest(x2[:,2], Normal(mu[2],S[2]))) > α
  @test pvalue(ExactOneSampleKSTest(x2[:,3], Normal(mu[3],S[3]))) > α
  @test maximum(abs.(cov(y[:,1:2])-cov(x2[:,1:2]))) < 0.03
end


text:
###################################################
## projectors/avgradius.jl
## implementation of network projector with radius projection
###################################################

"""
    `AvgRadius` : Projects trips to nodes within given radius, give the associated prediction
    If no available nodes within radius, projects to (at most) 3 nearest neighbors
"""
mutable struct AvgRadius <: NetworkProjector
    # compulsory attributes
    network::Network
    trips::GeoData

    "radius, in meters"
    radius::Float64
    "projection KD-tree (contains rides)"
    tree::KDTree
    "precomputation (node->trip list)"
    nodeList::Vector{Vector{Tuple{Int,Int}}}

    function AvgRadius(n::Network, r::Real)
        obj = new()
        obj.radius = r
        obj.network = n
        obj.trips = GeoData[]
        obj.nodeList = Vector{Tuple{Int,Int}}[]
        nNodePairs = length(n.nodes) * length(n.nodes)
        println("Creating KDTree...")
        dataPos = Array{Float64}(undef,(4, nNodePairs))
        for (i, startNode) in enumerate(n.nodes), (j, endNode) in enumerate(n.nodes)
            dataPos[1, (i-1) * length(n.nodes) + j] = startNode.x
            dataPos[2, (i-1) * length(n.nodes) + j] = startNode.y
            dataPos[3, (i-1) * length(n.nodes) + j] = endNode.x
            dataPos[4, (i-1) * length(n.nodes) + j] = endNode.y
        end
        obj.tree = KDTree(dataPos)
        return obj
    end
end

"""
    `AvgRadius` extended constructor: also preloads the data
"""
function AvgRadius(n::Network, r::Real, trips::GeoData)
    ar = AvgRadius(n, r)
    preloadData!(ar, trips)
    return ar
end

"""
    `preloadData!`: project all trips onto the node pairs within the given radius
"""
function preloadData!(ar::AvgRadius, trips::GeoData)
    nTrips = length(trips)
    println("Projecting $nTrips geo-trips onto network-nodes...")
    #initializing containers
    ar.trips = trips
    ar.nodeList = Array{Vector{Tuple{Int,Int}}}(undef, nTrips)

    for (i,t) in enumerate(trips)
        if i%100_000 == 0
            @printf("\r%.2f%% trips projected     ",100*i/nTrips)
        end
        pX, pY = toENU(t.pLon, t.pLat, ar.network)
        dX, dY = toENU(t.dLon, t.dLat, ar.network)
        tripLocation = [pX, pY, dX, dY]
        nodePairs = inrange(ar.tree, tripLocation, ar.radius)
        if length(nodePairs) == 0
            nodePairs, dists = knn(ar.tree, tripLocation, 3)
        end
        nodePairs = [decipherNodePairIndex(node, length(ar.network.nodes)) for node in nodePairs]
        ar.nodeList[i] = nodePairs[map(isValidNodePair, nodePairs)]
    end
    println("\r100.00% trips projected     ")
    return ar
end

"""
    `getNetworkTrips`: constructs list of NetworkTrip objects
"""
function getNetworkTrips(ar::AvgRadius, tIds::AbstractArray{Int64,1})
    if length(ar.trips) < length(tIds)
        error("Trips not loaded yet")
    end
    # Store result in dictionary : (orig,dest)->NetworkTrip
    netTrips = Dict{Tuple{Int,Int},NetworkTrip}()
    for id in tIds
        t = ar.trips[id]
        nodePairs = ar.nodeList[id]
        # compute weight of trip = inverse of number of ODs (network trips) it will be added to
        tripWeight = 1/length(nodePairs)
        for od in nodePairs
            #if same origin/destination
            if od[1] == od[2]
                continue
            end
            # check if network trip for this OD pair already exists
            if haskey(netTrips, od)
                nt = netTrips[od]
                # compute new time = weighted geometric average of times
                time = nt.time ^ (nt.weight/(nt.weight + tripWeight)) * t.time ^ (tripWeight/(nt.weight + tripWeight))
                # update time and wait
                netTrips[od] = NetworkTrip(nt.orig, nt.dest, time, nt.weight + tripWeight, false)
            else
                # create new network trip
                netTrips[od] = NetworkTrip((od[1], od[1], 0.), (od[2], od[2], 0.), t.time, tripWeight, false)
            end
        end
    end
    return collect(values(netTrips))
end

"""
    `getTripTiming`: returns a trip timing estimation from a NetworkTimings result
"""
function getTripTiming(ar::AvgRadius, timings::NetworkTimings, tId::Int)
    nodePairs = ar.nodeList[tId]
    # geometric average over neighboring nodes
    time = 1.0
    for od in nodePairs
        time *= timings.pathTimes[od[1], od[2]]
    end
    return time^(1.0/length(nodePairs))
end
function getTripTiming(ar::AvgRadius, timings::NetworkTimings,
        pLon::Float32, pLat::Float32, dLon::Float32, dLat::Float32)
    # find nodes within radius
    pX, pY = toENU(pLon, pLat, ar.network)
    dX, dY = toENU(dLon, dLat, ar.network)
    tripLocation = [pX, pY, dX, dY]
    nodePairs = inrange(ar.tree, tripLocation, ar.radius)
    nodePairs = [decipherNodePairIndex(node, length(ar.network.nodes)) for node in nodePairs]
    nodePairs = nodePairs[map(isValidNodePair, nodePairs)]
    # geometric average of travel times
    time = 1.0
    for od in nodePairs
        time *= timings.pathTimes[od[1], od[2]]
    end
    return time^(1.0/length(nodePairs))
end

"""
    `isValidNodePair`: helper function to check different start and end
"""
isValidNodePair(nodePair::Tuple{Int,Int}) = (nodePair[1] != nodePair[2])

"""
    `decipherNodePairIndex`: helper function to go from single index in KDTree to node pair
"""
function decipherNodePairIndex(idx::Int, nNodes::Int)
    startNode = Int(div(idx - 1, nNodes)) + 1
    endNode = idx % nNodes
    if endNode == 0
        endNode = nNodes
    end
    return (startNode, endNode)
end


text:
texstring{T<:Number}(a::T) = string(a)
texstring(a::Float64) = string(@sprintf("%1.2f",a))
texstring(a::Rational) = ( a.den == one(a) ? string(a.num) : string("\\frac{", a.num, "}{", a.den, "}") )

# convert a polynomial to tex format, e.g.,  2 x_1 x_2 + x_3
function latex(p::Poly)
    if p.m == 0 return("0") end

    s = ""
    for i=1:p.m
        if i==1
            s *= (p.c[i]>=0 ? "" : "- ")
        else
            s *= (p.c[i]>=0 ? "+ " : "- ")
        end

        ci = abs(p.c[i])
        if sum(p.alpha[i,:]) == 0
            s *= texstring(ci)
        else
            if ci != one(ci) s *= texstring(ci) * " " end
            for k=1:p.n
                if p.alpha[i,k] > 0
                    s *= p.syms.names[k]
                    if p.alpha[i,k] > 1
                        s *= "^{$(p.alpha[i,k])} "
                    else
                        s *= " "
                    end
                end
            end
        end
    end
    s
end

# convert a polynomial to vectorized tex format, e.g.,  2 y_{110} + y_{001}
function latex(p::Poly, imap::Dict{Array{Int,2},Int}, texsymbols::Array{String,1})
    if p.m == 0 return("0") end

    s = ""
    for i=1:p.m
        if i==1
            s *= (p.c[i]>=0 ? "" : "- ")
        else
            s *= (p.c[i]>=0 ? "+ " : "- ")
        end

        ci = abs(p.c[i])
        if sum(p.alpha[i,:]) == 0
            s *= texstring(ci) * " "
        else
            if ci != one(ci) s *= texstring(ci) * " " end
            s *= texsymbols[imap[p.alpha[i,:]]] * " "
        end
    end
    s
end

function latex{T}(M::Array{Poly{T}}, imap::Dict{Array{Int,2},Int}, texsymbols::Array{String,1}, matrixlimit = 10)
    n = size(M,1)
    s = "\\left[\n\\begin{array}{*{$(min(matrixlimit,n))}c}\n"
    if matrixlimit < n
        for i=1:matrixlimit-2
            s *= join([ [ latex(m, imap, texsymbols) for m=M[i,1:matrixlimit-2] ],
                            ["\\dots", latex(M[i,n], imap, texsymbols)] ], " & ") * "\\\\\n"
        end
        s *= join([ ["\\vdots" for k=1:matrixlimit-2], "\\ddots", "\\vdots\\\\\n"], " & ")

        s *= join([ [ latex(m, imap, texsymbols) for m=M[n,1:matrixlimit-2] ],
                        ["\\dots", latex(M[n,n], imap, texsymbols)] ], " & ") * "\\\\\n"
    else
        for i=1:n
            s *= join([ latex(m, imap, texsymbols) for m=M[i,:]], " & ") * "\\\\\n"
        end
    end
    string(s, "\\end{array}\n\\right]\n")
end

function latex_dual(prob::MomentProb, imap::Dict{Array{Int,2},Int}, texsymbols::Array{String,1}, matrixlimit = 10)

    s = string("Dual problem:\n\\[\n\\begin{array}{ll}\n\\text{minimize} &",
               latex(dot(prob.obj, prob.basis), imap, texsymbols), "\\\\\n")
    s *= "\\text{subject to} "

    for k=1:length(prob.mom)
        nk = int(sqrt(size(prob.mom[k],1)))
        Mk = reshape(prob.basis'*prob.mom[k]', (nk, nk))
        Mk = Mk + Mk' - diagm(diag(Mk))

        if nk > 1
            s = string(s, " & ", latex(Mk, imap, texsymbols, matrixlimit), "\\succeq 0\\\\\n\n")
        else
            s = string(s, " & ", latex(Mk[1], imap, texsymbols), " \\geq 0\\\\\n")
        end
    end

    for k=1:length(prob.eq)
        nk = int(sqrt(size(prob.eq[k],1)))
        Mk = reshape(prob.basis'*prob.eq[k]', (nk, nk))
        Mk = Mk + Mk' - diagm(diag(Mk))

        if nk > 1
            s = string(s, " & ", latex(Mk, imap, texsymbols, matrixlimit), " = 0\\\\\n\n")
        else
            s = string(s, " & ", latex(Mk[1], imap, texsymbols), " = 0\\\\\n")
        end
    end

    string(s, "\\end{array}\n\\]\n")
end

function latex_dual(prob::MomentProb, vectorizednames::Bool, matrixlimit::Int)

    imap = Polyopt.indexmap(prob.basis)
    syms = Array(String, length(imap))

    if vectorizednames == true
        for (key,k) = imap
            syms[k] = "y_{" * join([string(i) for i=key]) * "}"
        end
    else
        for k=1:length(imap)
            syms[k] = latex(prob.basis[k])
        end
    end

    latex_dual(prob, imap, syms, matrixlimit)
end

function latex_probstats(prob::MomentProb)
    sdpblksize  = [ int(sqrt(size(m,1))) for m=prob.mom ]
    numsdpblk   = length(sdpblksize)
    maxblksize  = maximum(sdpblksize)
    avgblksize  = int(mean(sdpblksize))

    eqsize      = [ int(sqrt(size(m,1))) for m=prob.eq ]
    numfree     = sum(Int[ n*(n+1) >> 1 for n=eqsize ])
    numvar      = sum([ n*(n+1) >> 1 for n=sdpblksize ]) + numfree
    numcon      = length(prob.basis)-1

    maxrhs  = maximum(prob.obj)
    maxcoef = zeros(length(prob.mom) + length(prob.eq))
    for k=1:length(prob.mom)
        maxcoef[k] = norm(prob.mom[k].nzval, Inf)
    end

    for k=1:length(prob.eq)
        maxcoef[k+numsdpblk] = norm(prob.eq[k].nzval, Inf)
    end

    s = "\\begin{quote}\n\\begin{tabular}{ll}\n"
    s = string(s,"Problem statistics\\\\\n\\hline\n")
    s = string(s,"constraints & $(numcon) \\\\\n")
    s = string(s,"sdp blocks & $(numsdpblk) \\\\\n")
    s = string(s,"largest block size & $(maxblksize) \\\\\n")
    s = string(s,"avg. block size & $(avgblksize) \\\\\n")
    s = string(s,"free variables & $(numfree) \\\\\n")
    s = string(s,"total number of scalar vars& $(numvar) \\\\\n")
    s = string(s, @sprintf("largest matrix coeff & %1.1e \\\\", maximum(maxcoef)))
    s = string(s, @sprintf("largest rhs & %1.1e\\\\\n", maxrhs))
    s = string(s, "\\hline\n")
    string(s,"\\end{tabular}\n\\end{quote}\n\n")
end

function latex_header()
    s = "\\documentclass{standalone}\n"
    s = string(s,"\\usepackage{varwidth}\n")
    s = string(s,"\\usepackage[paperwidth=\\maxdimen,paperheight=\\maxdimen]{geometry}\n")
    s = string(s,"\\usepackage{amsmath}\n")
    s = string(s,"\\begin{document}\n")
    string(s,"\\begin{varwidth}{\\linewidth}\n")
end

function latex_footer()
    s = "\\end{varwidth}\n"
    string(s, "\\end{document}")
end

function latex(texfile::String, prob::MomentProb, header="", body="", vectorizednames=false, matrixlimit=10)

    f = open(texfile,"w")

    tex = string(latex_header(),
                 header,
                 latex_probstats(prob),
                 body,
                 latex_dual(prob, vectorizednames, matrixlimit),
                 latex_footer())

    print(f, tex)
    close(f)
end


text:
using Test
using Sockets
import Sockets: TCPServer, close

import ProtoBuf: call_method, write_bytes, read_bytes

include("testsvc_pb.jl")

# Our RpcController for the test does nothing as of now
mutable struct TestRpcController <: ProtoRpcController
    debug::Bool
end

debug_log(controller::TestRpcController, msg) = controller.debug && println(msg)
error_log(controller::TestRpcController, msg) = println(stderr, msg)

# RpcChannel implementation for our test protocol
# The protocol is to write and read:
# - messages as delimited bytes
# - header type SvcHeader to identify the method being called
mutable struct TestRpcChannel <: ProtoRpcChannel
    sock::TCPSocket
end
close(channel::TestRpcChannel) = close(channel.sock)

mutable struct SvcHeader <: ProtoType
    method::String
    SvcHeader() = (o=new(); fillunset(o); o)
end

function write_request(channel::TestRpcChannel, controller::TestRpcController, service::ServiceDescriptor, method::MethodDescriptor, request)
    io = channel.sock
    hdr = SvcHeader()
    hdr.method = method.name

    iob = IOBuffer()

    hdr_len = writeproto(iob, hdr)
    hdr_buff = take!(iob)

    data_len = writeproto(iob, request)
    data_buff = take!(iob)

    write_bytes(io, hdr_buff)
    #debug_log(controller, "req hdr  ==> $hdr_len bytes: $hdr_buff")
    write_bytes(io, data_buff)
    #debug_log(controller, "req data ==> $data_len bytes: $data_buff")
    hdr_len + data_len
end

function read_request(channel::TestRpcChannel, controller::TestRpcController, srvr)
    io = channel.sock
    hdr = SvcHeader()
    hdr_buff = read_bytes(io)
    #debug_log(controller, "req hdr  <== $(length(hdr_buff)) bytes: $hdr_buff")
    readproto(IOBuffer(hdr_buff), hdr)
    method = find_method(srvr, hdr.method)

    request_type = get_request_type(srvr, method)
    request = request_type()
    data_buff = read_bytes(io)
    #debug_log(controller, "req data <== $(length(data_buff)) bytes: $data_buff")
    readproto(IOBuffer(data_buff), request)
    method, request
end

function write_response(channel::TestRpcChannel, controller::TestRpcController, response)
    io = channel.sock
    iob = IOBuffer()
    data_len = writeproto(iob, response)
    data_buff = take!(iob)
    write_bytes(io, data_buff)
    #debug_log(controller, "resp ==> $(length(data_buff)) bytes: $data_buff")
    data_len
end

function read_response(channel::TestRpcChannel, controller::TestRpcController, response)
    io = channel.sock
    data_buff = read_bytes(io)
    #debug_log(controller, "resp <== $(length(data_buff)) bytes: $data_buff")
    readproto(IOBuffer(data_buff), response)
    response
end

function call_method(channel::TestRpcChannel, service::ServiceDescriptor, method::MethodDescriptor, controller::TestRpcController, request)
    write_request(channel, controller, service, method, request)
    response_type = get_response_type(method)
    response = response_type()
    read_response(channel, controller, response)
end

# Test server implementation on the RpcChannel
#
mutable struct TestServer
    srvr::TCPServer
    impl::ProtoService
    run::Bool
    debug::Bool
    TestServer(port::Integer) = new(listen(port), TestMath(Main), true, false)
end

function process(srvr::TestServer, channel::TestRpcChannel)
    controller = TestRpcController(srvr.debug)
    #debug_log(controller, "starting processing from channel")

    try
        while(!eof(channel.sock))
            method, request = read_request(channel, controller, srvr.impl)
            response = call_method(srvr.impl, method, controller, request)
            #debug_log(controller, "response: $response")
            write_response(channel, controller, response)
        end
    catch ex
        debug_log(controller, "channel stopped with exception $ex")
    end
    #debug_log(controller, "stopped processing channel")
end

# implementations of our test services
function Add(req::BinaryOpReq)
    resp = BinaryOpResp()
    resp.result = req.i1 + req.i2
    resp
end

function Mul(req::BinaryOpReq)
    resp = BinaryOpResp()
    resp.result = req.i1 * req.i2
    resp
end

# Utility methods for running the server and client
function run_server(srvr::TestServer)
    controller = TestRpcController(srvr.debug)
    try
        while(srvr.run)
            sock = accept(srvr.srvr)
            channel = TestRpcChannel(sock)
            @async process(srvr, channel)
        end
    catch ex
        debug_log(controller, "server stopped with exception $ex")
    end
    debug_log(controller, "stopped server")
end

global nresults = 0
function chk_results(out::BinaryOpResp, expected, channel=nothing)
    global nresults
    if channel != nothing
        close(channel)
        nresults += 1
    end
    @test out.result == expected
    nothing
end

function run_client(debug::Bool)
    global nresults
    controller = TestRpcController(debug)

    debug_log(controller, "testing services...")
    debug_log(controller, "testing blocking stub...")
    let channel=TestRpcChannel(connect(9999)), stub=TestMathBlockingStub(channel)
        for i in 1:10
            inp = BinaryOpReq()
            inp.i1 = Int64(rand(Int8))
            inp.i2 = Int64(rand(Int8))

            out = Add(stub, controller, inp)
            chk_results(out, inp.i1+inp.i2)

            out = Mul(stub, controller, inp)
            chk_results(out, inp.i1*inp.i2)
        end
        close(channel)
    end

    debug_log(controller, "testing non blocking stub...")
    for i in 1:10
        inp = BinaryOpReq()
        inp.i1 = Int64(rand(Int8))
        inp.i2 = Int64(rand(Int8))

        nresults -= 1
        let channel=TestRpcChannel(connect(9999)), stub=TestMathStub(channel), expected=inp.i1+inp.i2
            Add(stub, controller, inp, (out)->chk_results(out, expected, channel))
        end

        nresults -= 1
        let channel=TestRpcChannel(connect(9999)), stub=TestMathStub(channel), expected=inp.i1*inp.i2
            Mul(stub, controller, inp, (out)->chk_results(out, expected, channel))
        end
    end
    while nresults != 0
        debug_log(controller, "waiting for $(abs(nresults)) results")
        yield()
        sleep(1)
    end
end

debug = true
srvr = TestServer(9999)
srvr.debug = debug
@async run_server(srvr)
sleep(1)
run_client(debug)
srvr.run = false
close(srvr.srvr)


text:
#!/usr/bin/env julia
# -*- coding: utf-8 -*-
#=
Apply HANTS process on a multidimensional array
=#
cd(dirname(@__FILE__))
include("../HANTS.jl")
using .HANTS, Plots, DelimitedFiles

arr = readdlm("test_data.csv", ',', Float64, '\n')
arr = reshape(arr, 4, 5, 23) # 3d array, the 3rd dimension to apply HANTS

nbase = 365
nfreq = 3
tseries = 1:16:365
validrange = (-1, 1)
fet = 0.02
dod = 1
δ = 0.1

A, φ, arr_rec = hants(
    arr, fet, dod, δ; dims=3, nbase=nbase, nfreq=nfreq, validrange=validrange, tseries=tseries, outlier=nothing
)
arr_reconstruct = reconstruct(A, φ, nbase; dims=3)

A_Hi, φ_Hi, arr_rec_Hi = hants(
    arr, fet, dod, δ; dims=3, nbase=nbase, nfreq=nfreq, validrange=validrange, tseries=tseries, outlier=:Hi
)
arr_reconstruct_Hi = reconstruct(A_Hi, φ_Hi, nbase; dims=3)

A_Lo, φ_Lo, arr_rec_Lo = hants(
    arr, fet, dod, δ; dims=3, nbase=nbase, nfreq=nfreq, validrange=validrange, tseries=tseries, outlier=:Lo
)
arr_reconstruct_Lo = reconstruct(A_Lo, φ_Lo, nbase; dims=3)

fig1 = plot(tseries, arr[2, 4, :]; lc=:black, label="Original Data", shape=:circle, mc=:black)
plot!(tseries, arr_rec[2, 4, :]; lc=:blue, label="HANTS - None", shape=:circle, mc=:blue)
plot!(tseries, arr_rec_Hi[2, 4, :]; lc=:red, label="HANTS - Hi", shape=:circle, mc=:red)
plot!(tseries, arr_rec_Lo[2, 4, :]; lc=:green, label="HANTS - Lo", shape=:circle, mc=:green)
title!("Testing HANTS Algorithm")

fig2 = plot(tseries, arr[2, 4, :]; lc=:black, label="Original Data", shape=:circle, mc=:black)
plot!(1:365, arr_reconstruct[2, 4, :]; lc=:blue, label="HANTS - None")
plot!(1:365, arr_reconstruct_Hi[2, 4, :]; lc=:red, label="HANTS - Hi")
plot!(1:365, arr_reconstruct_Lo[2, 4, :]; lc=:green, label="HANTS - Lo")
title!("Testing HANTS Reconstruct")

plot(fig1, fig2; layout=(1, 2), size=(800,300))

max_stars_repo_name: shenrq/JuliaHANTS


text:
@testset "Counter" begin
    counter = EvaluationCounter()

    @test evaluations(counter) ≡ 0

    @test increment!(counter, -1) ≡ 0

    @test increment!(counter) ≡ evaluations(counter) ≡ 1

    restricted_type_test(counter)
end

@testset "Timer" begin
    timer = EvaluationTimer()

    @test Δt(timer) ≡ 0.0

    on!(timer)
    sleep(0.5)
    off!(timer)

    @test Δt(timer) > 0.5

    restricted_type_test(timer)
end

@testset "Profile" begin
    profile = BlockOptProfile()

    @test Δt(trs_timer(profile)) ≡ Δt(ghs_timer(profile)) ≡ 0.0

    @test evaluations(trs_counter(profile)) ≡ evaluations(ghs_counter(profile)) ≡ 0

    restricted_type_test(profile)
end

@testset "Weaver" begin
    weaver = Weaver()

    @test weave_level(weaver) ≡ WEAVE_LEVEL

    weave!(weaver, :f_vals, 0.0)

    @test pop!(f_vals(weaver)) ≡ 0.0

    weave!(weaver, :∇f_norms, 0.0)

    @test pop!(∇f_norms(weaver)) ≡ 0.0

    weave!(weaver, :p_norms, 0.0)

    @test pop!(p_norms(weaver)) ≡ 0.0

    weave!(weaver, :Δ_vals, 0.0)

    @test pop!(Δ_vals(weaver)) ≡ 0.0

    weave!(weaver, :ρ_vals, 0.0)

    @test pop!(ρ_vals(weaver)) ≡ 0.0
end

@testset "Trace" begin

    test_trace = BlockOptTrace(test_model, test_driver)

    @testset "constructor test" begin

        @test model(test_trace) ≡ test_model

        @test driver(test_trace) ≡ test_driver

        @test isa(profile(test_trace), BlockOptProfile)

        @test isa(trs_timer(test_trace), EvaluationTimer)

        @test isa(ghs_timer(test_trace), EvaluationTimer)

        @test isa(trs_counter(test_trace), EvaluationCounter)

        @test isa(ghs_counter(test_trace), EvaluationCounter)

        @test isa(weaver(test_trace), Weaver)

        @test map(
            isempty,
            (
                f_vals(test_trace),
                ∇f_norms(test_trace),
                Δ_vals(test_trace),
                p_norms(test_trace),
                ρ_vals(test_trace),
            ),
        ) ≡ (true, true, true, true, true)

        @test weave_level(test_trace) ≡ weave_level(test_driver)

        @test log_level(test_trace) ≡ log_level(test_driver)

        @test isfile(io(test_trace))

        restricted_type_test(test_trace)
    end

    @testset "generated logging function test" begin

        # TODO: Visual confrimation suggests working as expected

    end

end


text:
#################### Metropolis-Adjusted Langevin Algorithm ####################

#################### Types ####################

mutable struct MALATune <: SamplerTune
  logfgrad::Union{Function, Missing}
  epsilon::Float64
  SigmaL::Union{UniformScaling{Bool}, LowerTriangular{Float64}}

  MALATune() = new()

  MALATune(x::Vector, epsilon::Real) = new(missing, epsilon, I)

  MALATune(x::Vector, epsilon::Real, logfgrad::Function) =
    new(logfgrad, epsilon, I)

  MALATune(x::Vector, epsilon::Real, Sigma::Matrix{T}) where {T<:Real} =
    new(missing, epsilon, cholesky(Sigma).L)

  function MALATune(x::Vector, epsilon::Real, Sigma::AbstractMatrix{T},
                    logfgrad::Function) where {T<:Real}
    new(logfgrad, epsilon, cholesky(Sigma).L)
  end
end

const MALAVariate = SamplerVariate{MALATune}

validate(v::MALAVariate) = validate(v, v.tune.SigmaL)

validate(v::MALAVariate, SigmaL::UniformScaling) = v

function validate(v::MALAVariate, SigmaL::LowerTriangular)
  n = length(v)
  size(SigmaL, 1) == n ||
    throw(ArgumentError("Sigma dimension differs from variate length $n"))
  v
end

#################### Sampler Constructor ####################
"""
    MALA(params::ElementOrVector{Symbol}, epsilon::Real; args...)

Construct a `Sampler` object for MALA sampling. Parameters are assumed to be
continuous, but may be constrained or unconstrained.

Returns a `Sampler{MALATune}`` type object.

* `params`: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochastic `unlist()` function.

* `epsilon`: factor by which the drift and covariance matrix of the proposal distribution are scaled.

* `Sigma`: covariance matrix for the multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated. If omitted, the identity matrix is assumed.

* `dtype` : differentiation for gradient calculations. Options are
    * `:central` : central differencing
    * `:forward` : forward differencing.
"""
function MALA(params::ElementOrVector{Symbol}, epsilon::Real; args...)
  MALASampler(params, epsilon; args...)
end
"""
    MALA(params::ElementOrVector{Symbol}, epsilon::Real, Sigma::Matrix{T}; args...)

Construct a `Sampler` object for MALA sampling. Parameters are assumed to be
continuous, but may be constrained or unconstrained.

Returns a `Sampler{MALATune}`` type object.
"""
function MALA(params::ElementOrVector{Symbol}, epsilon::Real,
               Sigma::Matrix{T}; args...) where {T<:Real}
  MALASampler(params, epsilon, Sigma; args...)
end

function MALASampler(params, pargs...; dtype::Symbol=:forward)
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block, true)
    v = SamplerVariate(block, pargs...)
    sample!(v, x -> logpdfgrad!(block, x, dtype))
    relist(block, v)
  end
  Sampler(params, samplerfx, MALATune())
end

#################### Sampling Functions ####################

sample!(v::MALAVariate) = sample!(v, v.tune.logfgrad)
"""
    sample!(v::MALAVariate, logfgrad::Function)

Draw one sample from a target distribution using the MALA sampler. Parameters
are assumed to be continuous and unconstrained.

Returns `v` updated with simulated values and associated tuning parameters.
"""
function sample!(v::MALAVariate, logfgrad::Function)
  tune = v.tune

  L = sqrt(tune.epsilon) * tune.SigmaL
  Linv = inv(L)
  M2 = 0.5 * L * L'

  logf0, grad0 = logfgrad(v.value)
  y = v + M2 * grad0 + L * randn(length(v))
  logf1, grad1 = logfgrad(y)

  q0 = -0.5 * sum(abs2, Linv * (v - y - M2 * grad1))
  q1 = -0.5 * sum(abs2, Linv * (y - v - M2 * grad0))

  if rand() < exp((logf1 - q1) - (logf0 - q0))
    v[:] = y
  end

  v
end


text:
# Configuration for nearly neutral simuation
const dfe_s_list = [0.0,-0.01]
dfe_str = "dfe_neut_fixed"
const nn_simtype = 1    # nearly neutral infinite alleles
const mu_list_flag = false
@everywhere const N_list     = [25,50,100,200,400,800,1600,3200,6400]     # popsize
const N_mu_list= [0.5,1.0,2.0,4.0]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci
const ngens = 200000      # Generations after burn-in
const fix_minimum = 0.45  # Fraction of population needed to "fix"
const burn_in= 2.0    # generations of burn_in as a multiple of 1/mu

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
## deprecated functions..

function abc_mcmc_dbscan_full(dfcart::GaiaClustering.Df, params::abcfull)
    let
        println("# abc_mcmc_dbscan_full is deprecated ...")

        Random.seed!()
        println("## ABC/MCMC for DBSCAN FULL (parameters+weighting)...")
        mci = mcfull(zeros(Float64,0),zeros(Int32,0),zeros(Int32,0) , zeros(Float64,0), zeros(Int32,0),
        zeros(Float64,0), zeros(Float64,0), zeros(Float64,0))

        initial = true
        th = []
        mi = modelfull(0.0,0,0, 0., 0., 0.)
        micurrent = modelfull(0.0,0,0, 0., 0., 0.)
        probi = 0.
        probcurrent = 0.

        minimumQ = params.minQ
        minstars = params.minstars
        niter = params.niter
        nburn = params.nburnout

        maxstars= 5000

        println("### Minimum Q : $minimumQ")
        println("### Minimum nstars : $minstars")
        println("### Maximum nstars : $maxstars")
        minimumQ , minstars = check_qminqstar_full(dfcart, params, minimumQ, minstars)
        println("### Minimum Q : $minimumQ")
        println("### Minimum nstars : $minstars")
        if minstars < params.forcedminstars
            minstars = params.forcedminstars
            println("### Minimum nstars forced to : $minstars")
        end

        while initial
            mi, probi = theta_full(params)
            dfcartnorm = getDfcartnorm(dfcart , mi)
            qres , nstars = find_clusters(dfcartnorm, dfcart, mi)
            if qres > minimumQ && nstars >= minstars && nstars <= maxstars
                println("### init done ...")
                initial = false
                push!(mci.eps, mi.eps)
                push!(mci.mne, mi.min_nei)
                push!(mci.mcl, mi.min_cl)
                push!(mci.w3d, mi.w3d)
                push!(mci.wvel, mi.wvel)
                push!(mci.whrd, mi.whrd)
                push!(mci.qc,  qres)
                push!(mci.qn,  nstars)
            end
        end

        ministats_full(100, dfcart, mi , params)

        nchain = 1
        loopAgain = true
        burndone = false

        while loopAgain
            micurrent, probcurrent = thetaiter_full(mi , params)
            dfcartnorm = getDfcartnorm(dfcart , micurrent)
            qres , nstars = find_clusters(dfcartnorm, dfcart, micurrent)

            if qres > minimumQ && nstars >= minstars && nstars <= maxstars
            ### Metropolis-Hasting
                α = probcurrent / probi
                if α > rand()
                    mi = micurrent
                    probi = probcurrent
                    nchain += 1
                    if (nchain%1000 == 0) println("### chain:",nchain) end
                    if nchain > nburn && !burndone println("### burnout done...") ; nchain = 0 ; burndone = true end
                    if nchain > niter loopAgain = false end
                    if burndone
                        push!(mci.eps, mi.eps)
                        push!(mci.mne, mi.min_nei)
                        push!(mci.mcl, mi.min_cl)
                        push!(mci.w3d, mi.w3d)
                        push!(mci.wvel, mi.wvel)
                        push!(mci.whrd, mi.whrd)
                        push!(mci.qc , qres)
                        push!(mci.qn,  nstars)
                    end
                else
                    nchain += 1
                    if (nchain%1000 == 0) println("### chain:",nchain) end
                    if nchain > nburn && !burndone println("### burnout done...") ; nchain = 0 ; burndone = true end
                    if nchain > niter loopAgain = false end
                    if burndone
                        push!(mci.eps, mi.eps)
                        push!(mci.mne, mi.min_nei)
                        push!(mci.mcl, mi.min_cl)
                        push!(mci.w3d, mi.w3d)
                        push!(mci.wvel, mi.wvel)
                        push!(mci.whrd, mi.whrd)
                        push!(mci.qc , qres)
                        push!(mci.qn,  nstars)
                    end
                end
            end
        end
        println("## ABC/MCMC FULL done")
        println("##")
        return(mci)
    end
end

### function to adapt the minimum condition to evaluate the  solution
###
function check_qminqstar_full(dfcart::GaiaClustering.Df,
        params::abcfull, minimumQ, minstars)
    let
        new_minq = minimumQ
        new_minstars = minstars
        notfound = true

        ## standard sofar 50
        mingoodsolution = 10
        niter = 500

        println("#### Checking the minQ and minStars conditions...")
        println("#### Minimum good solutions $mingoodsolution")

        while notfound
            goodsolutions = 0
            for i in 1:niter
                mi, probi = theta_full(params)
                dfcartnorm = getDfcartnorm(dfcart , mi)
                qres , nstars = find_clusters(dfcartnorm, dfcart, mi)
                if qres > new_minq && nstars >= new_minstars
                    goodsolutions += 1
                end
                if goodsolutions > mingoodsolution
                    notfound = false
                    return(new_minq, new_minstars)
                end
            end
            if notfound
                new_minq *= 0.95
                new_minstars = trunc(Int, 0.95 * new_minstars)
                println("#### MinQ not reached yet... testing with $new_minq")
            end

            if new_minstars == 0
                println("#### Problem with minStars, return(0.5,5)")
                return(0.5,5)
            end
        end
        return(new_minq, new_minstars)
    end
end
function find_clusters(df::GaiaClustering.Df, dfcart::GaiaClustering.Df , m::GaiaClustering.model,
    aperture2d = 1.5, maxaperture2d = 15, aperturev = 3.0, maxaperturev = 20, nboot = 50 ,
    aperture3d = 3., maxaperture3d = 20)
    let
        println("Warning...find_clusters deprecated, find_clusters2 should be used instead")
        labels = clusters(df.data , m.eps , 20, m.min_nei, m.min_cl)
        nsol= length(labels)
        if nsol == 0
            return(0, 0)
        end
        if nsol > 1000
            println("### Warning... $nsol clusters found in find_clusters ")
        end

    ### metrics of the clusters
        q2d = metric2(dfcart, labels, "spatial2d" , aperture2d, maxaperture2d, nboot)
        q3d = metric2(dfcart, labels, "spatial3d" , aperture3d, maxaperture3d, nboot)     #### Added
        qv  = metric2(dfcart, labels, "velocity" , aperturev, maxaperturev, nboot)
        qp, qa = metric2(dfcart, labels, "HRD" )

        nlab = []
        for ilab in labels
            push!(nlab,length(ilab))
        end

    #### metric for the number of stars in the cluster
        qn = []
        for nl in nlab
            push!(qn,log10(nl))
        end

        qc = 0.
        qstar = 0
        for i in 1:length(nlab)
            k1 = q2d[i][1]
            k1bis = q3d[i][1]
            k2 = qv[i][1]
            k3 = qa[i][1]
            k4 = qn[i]
            ############### Composite metric ###
            qq = (2k1 + k1bis + 3k2 + k3 + k4) / 8.0
            # qq = (3k1 + k1bis + 3k2 + k4) / 8.0
            ###############
            if qq > qc
                qc = qq
                qstar = nlab[i]
            end
        end
        return(qc, qstar)
    end
end

## broadcasting modelfull ..
##
function find_clusters(df::GaiaClustering.Df, dfcart::GaiaClustering.Df , m::GaiaClustering.modelfull,
    aperture2d = 1.5, maxaperture2d = 15, aperturev = 3.0, maxaperturev = 20, nboot = 50,
    aperture3d = 3., maxaperture3d = 20)
    let
        println("Warning...find_clusters deprecated, find_clusters2 should be used instead")
        labels = clusters(df.data , m.eps , 20, m.min_nei, m.min_cl)
        nsol= length(labels)
        if nsol == 0 || nsol > 1000
            println("### Warning... $nsol clusters found in find_clusters2 ")
            return(0, 0)
        end

    ### metrics of the clusters
        q2d = metric2(dfcart, labels, "spatial2d" , aperture2d, maxaperture2d, nboot)
        q3d = metric2(dfcart, labels, "spatial3d" , aperture3d, maxaperture3d, nboot)
        qv  = metric2(dfcart, labels, "velocity" , aperturev, maxaperturev, nboot)
        qp, qa = metric2(dfcart, labels, "HRD" )

        nlab = []
        for ilab in labels
            push!(nlab,length(ilab))
        end

    #### metric for the number of stars in the cluster
        qn = []
        for nl in nlab
            push!(qn,log10(nl))
        end

        qc = 0.
        qstar = 0
        for i in 1:length(nlab)
            k1 = q2d[i][1]
            k1bis = q3d[i][1]
            k2 = qv[i][1]
            k3 = qa[i][1]
            k4 = qn[i]
            ############### Composite metric ###
            qq = (2k1 + k1bis + 3k2 + k3 + k4) / 8.0
            # qq = (3k1 + k1bis + 3k2 + k4) / 8.0
            ###############
            if qq > qc
                qc = qq
                qstar = nlab[i]
            end
        end
        return(qc, qstar)
    end
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
using Spring
using Test


text:
using TailRiskScenGen
using Distributions

# Inputs - change as necessary
dim = 20
gens = 30
samples = 100000
srand(1)

# Randomly generate distribution and cone
Z = Normal()
μ = rand(Z, dim)
A = rand(Z, (dim,dim))
Σ = A'A
dist = MvNormal(μ, Σ)
K = FiniteCone(rand(Z, (dim, gens)))

function test(dist::MvNormal, K::FiniteCone, samples::Int64)
    dim = K.dim
    for i in 1:samples
        project(K, rand(dim))
    end
end

# Small test to compile necessary functions
test(dist, K, 5)

# Real test
@time test(dist, K, samples)


text:
@testset "139.word-break.jl" begin
    tcase1 = ("leetcode", ["leet", "code"])
    tcase2 = ("applepenapple", ["apple", "pen"])
    tcase3 = ("bb", ["a", "b", "bbb", "bbbb"])
    fcase1 = ("catsandog", ["cats", "dog", "sand", "and", "cat"])
    fcase2 = ("a", ["b"])
    fcase3 = ("a"^150 * "b", ["a"^i for i in 1:10])
    @test word_break(tcase1...) && word_break_bfs(tcase1...)
    @test word_break(tcase2...) && word_break_bfs(tcase2...)
    @test word_break(tcase3...) && word_break_bfs(tcase3...)
    @test !word_break(fcase1...) && !word_break_bfs(fcase1...)
    @test !word_break(fcase2...) && !word_break_bfs(fcase2...)
    @test !word_break(fcase3...) && !word_break_bfs(fcase3...)
end