text:
# Copyright 2017 Stephen L. Smith and Frank Imeson
#
# 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.

"""
Select a removal and an insertion method using powers, and then perform
removal followed by insertion on tour.  Operation done in place.
"""
function remove_insert(current::Tour, best::Tour, dist::Array{Int64,2}, member::Array{Int64,1},
						setdist::Distsv, sets::Array{Any,1},
						powers, param::Dict{Symbol,Any}, phase::Symbol)
	# make a new tour to perform the insertion and deletion on
    trial = Tour(copy(current.tour), current.cost)
	pivot_tour!(trial.tour)
	num_removals = rand(param[:min_removals]:param[:max_removals])

	removal = power_select(powers["removals"], powers["removal_total"], phase)
	if removal.name == "distance"
		sets_to_insert = distance_removal!(trial.tour, dist, num_removals,
													member, removal.value)
    elseif removal.name == "worst"
		sets_to_insert = worst_removal!(trial.tour, dist, num_removals,
													member, removal.value)
	else
		sets_to_insert = segment_removal!(trial.tour, num_removals, member)
	end

	randomize_sets!(sets, sets_to_insert)
	# then perform insertion
	insertion = power_select(powers["insertions"], powers["insertion_total"], phase)
	noise = power_select(powers["noise"], powers["noise_total"], phase)
	if insertion.name == "cheapest"
		cheapest_insertion!(trial.tour, sets_to_insert, dist, setdist, sets)
	else
		randpdf_insertion!(trial.tour, sets_to_insert, dist, setdist, sets,
							insertion.value, noise)
	end
	rand() < param[:prob_reopt] && opt_cycle!(trial, dist, sets, member, param, setdist, "partial")

	# update power scores for remove and insert
	score = 100 * max(current.cost - trial.cost, 0)/current.cost
	insertion.scores[phase] += score
	insertion.count[phase] += 1
	removal.scores[phase] += score
	removal.count[phase] += 1
	noise.scores[phase] += score
	noise.count[phase] += 1
	return trial
end

"""
Select an integer between 1 and num according to
and exponential distribution with lambda = power
# goes from left of array if power is positive
# and right of array if it is negative
"""
function select_k(num::Int64, power::Float64)
	base = (1/2)^abs(power)
	# (1 - base^num)/(1 - base) is sum of geometric series
	rand_select = (1 - base^num)/(1 - base) * rand()
	bin = 1
	@inbounds for k = 1:num
		if rand_select < bin
			return (power >= 0 ? (num - k + 1) : k)
		end
		rand_select -= bin
		bin *= base
	end
	return (power >=0 ? num : 1)
end

"""
selecting a random k in 1 to length(weights) according to power
and then selecting the kth smallest element in weights
"""
function pdf_select(weights::Array{Int64,1}, power::Float64)
    power == 0.0 && return rand(1:length(weights))
    power > 9.0 && return rand_select(weights, maximum(weights))
    power < - 9.0 && return rand_select(weights, minimum(weights))

	# select kth smallest.  If 1 or length(weights), simply return
	k = select_k(length(weights), power)
	k == 1 && return rand_select(weights, minimum(weights))
	k == length(weights) && return rand_select(weights, maximum(weights))
	val = partialsort(weights, k)

	return rand_select(weights, val)
end

"""  choose set with pdf_select, and then insert in best place with noise  """
function randpdf_insertion!(tour::Array{Int64,1}, sets_to_insert::Array{Int64,1},
							dist::Array{Int64, 2}, setdist::Distsv,
							sets::Array{Any, 1}, power::Float64, noise::Power)

    mindist = [typemax(Int64) for i=1:length(sets_to_insert)]
    @inbounds for i = 1:length(sets_to_insert)
        set = sets_to_insert[i]
        for vertex in tour
            if setdist.min_sv[set, vertex] < mindist[i]
                mindist[i] = setdist.min_sv[set, vertex]
            end
        end
    end
    new_vertex_in_tour = 0

    @inbounds while length(sets_to_insert) > 0
        if new_vertex_in_tour != 0
            for i = 1:length(sets_to_insert)
                set = sets_to_insert[i]
                if setdist.min_sv[set, new_vertex_in_tour] < mindist[i]
                    mindist[i] = setdist.min_sv[set, new_vertex_in_tour]
                end
            end
        end
        set_index = pdf_select(mindist, power) # select set to insert from pdf
        # find the closest vertex and the best insertion in that vertex
        nearest_set = sets_to_insert[set_index]
		if noise.name == "subset"
			bestv, bestpos = insert_subset_lb(tour, dist, sets[nearest_set], nearest_set,
											  setdist, noise.value)
		else
			bestv, bestpos =
					insert_lb(tour, dist, sets[nearest_set], nearest_set, setdist, noise.value)
		end
        insert!(tour, bestpos, bestv)  # perform the insertion
        new_vertex_in_tour = bestv
        # remove the inserted set from data structures
        splice!(sets_to_insert, set_index)
        splice!(mindist, set_index)
    end
end

function cheapest_insertion!(tour::Array{Int64,1}, sets_to_insert::Array{Int64,1},
	dist::Array{Int64, 2}, setdist::Distsv, sets::Array{Any, 1})
    """
	choose vertex that can be inserted most cheaply, and insert it in that position
	"""
	while length(sets_to_insert) > 0
        best_cost = typemax(Int64)
        best_v = 0
        best_pos = 0
        best_set = 0
        for i = 1:length(sets_to_insert)
            set_ind = sets_to_insert[i]
            # find the best place to insert the vertex
            best_v, best_pos, cost = insert_cost_lb(tour, dist, sets[set_ind], set_ind, setdist,
										  best_v, best_pos, best_cost)
			if cost < best_cost
				best_set = i
				best_cost = cost
			end
        end

        # now, perform the insertion
        insert!(tour, best_pos, best_v)
        # remove the inserted set from data structures
        splice!(sets_to_insert, best_set)
    end
end

"""
Given a tour and a set, this function finds the vertex in the set with minimum
insertion cost, along with the position of this insertion in the tour.  If
best_position is i, then vertex should be inserted between tour[i-1] and tour[i].
"""
@inline function insert_lb(tour::Array{Int64,1}, dist::Array{Int64,2}, set::Array{Int64, 1},
							setind::Int, setdist::Distsv, noise::Float64)
	best_cost = typemax(Int64)
	bestv = 0
	bestpos = 0

	@inbounds for i = 1:length(tour)
		v1 = prev_tour(tour, i)
		lb = setdist.vert_set[v1, setind] + setdist.set_vert[setind, tour[i]] - dist[v1, tour[i]]
		lb > best_cost && continue

		for v in set
	        insert_cost = dist[v1, v] + dist[v, tour[i]] - dist[v1, tour[i]]
			noise > 0.0 && (insert_cost += round(Int64, noise * rand() * abs(insert_cost)))
			if insert_cost < best_cost
				best_cost = insert_cost
				bestv = v
				bestpos = i
	        end
		end
    end
    return bestv, bestpos
end

@inline function insert_subset_lb(tour::Array{Int64,1}, dist::Array{Int64,2}, set::Array{Int64, 1},
							setind::Int, setdist::Distsv, noise::Float64)
	best_cost = typemax(Int64)
	bestv = 0
	bestpos = 0
	tour_inds = collect(1:length(tour))

	@inbounds for i = 1:ceil(Int64, length(tour) * noise)
		i = incremental_shuffle!(tour_inds, i)
		v1 = prev_tour(tour, i)
		lb = setdist.vert_set[v1, setind] + setdist.set_vert[setind, tour[i]] - dist[v1, tour[i]]
		lb > best_cost && continue

		for v in set
	        insert_cost = dist[v1, v] + dist[v, tour[i]] - dist[v1, tour[i]]
	        if insert_cost < best_cost
				best_cost = insert_cost
				bestv = v
				bestpos = i
	        end
		end
    end
    return bestv, bestpos
end

############ Initial Tour Construction ##########################

"""build tour from scratch on a cold restart"""
function initial_tour!(lowest::Tour, dist::Array{Int64, 2}, sets::Array{Any, 1},
						setdist::Distsv, trial_num::Int64, param::Dict{Symbol,Any})
	sets_to_insert = collect(1:param[:num_sets])
	best = Tour(Int64[], typemax(Int64))

	# compute random initial tour only if past first trial
	# in this case, randomly choose between random and insertion tour.
	if param[:init_tour] == "warm" && param[:warm_tour] != nothing
 		best.tour = param[:warm_tour]
	elseif param[:init_tour] == "rand" && (trial_num > 1) && (rand() < 0.5)
		random_initial_tour!(best.tour, sets_to_insert, dist, sets)
	else
		random_insertion!(best.tour, sets_to_insert, dist, sets, setdist)
	end
	best.cost = tour_cost(best.tour, dist)
	lowest.cost > best.cost && (lowest = best)
	return best
end

"""
Randomly shuffle the sets, and then insert the best vertex from each set back into
the tour where sets are considered in shuffled order.
"""
function random_insertion!(tour::Array{Int64,1}, sets_to_insert::Array{Int64,1},
						   dist::Array{Int64, 2}, sets::Array{Any, 1}, setdist::Distsv)
    shuffle!(sets_to_insert)  # randomly permute the sets
    for set in sets_to_insert
        # only have to compute the insert cost for the changed portion of the tour
        if isempty(tour)
            best_vertex = rand(sets[set])
            best_position = 1
        else
            best_vertex, best_position = insert_lb(tour, dist, sets[set], set, setdist, 0.75)
        end
        # now, perform the insertion
        insert!(tour, best_position, best_vertex)
    end
end

"""
Randomly shuffle the sets, and then insert the best vertex from each set back into
the tour where sets are considered in shuffled order.
"""
function random_initial_tour!(tour::Array{Int64,1}, sets_to_insert::Array{Int64,1},
							  dist::Array{Int64, 2}, sets::Array{Any, 1})
    shuffle!(sets_to_insert)
    for set in sets_to_insert
		push!(tour, rand(sets[set]))
    end
end

######################### Removals ################################
"""
Remove the vertices randomly, but biased towards those that add the most length to the
tour.  Bias is based on the power input.  Vertices are then selected via pdf select.
"""
function worst_removal!(tour::Array{Int64,1}, dist::Array{Int64, 2},
							num_to_remove::Int64, member::Array{Int64,1}, power::Float64)
    deleted_sets = Array{Int}(undef, 0)
	while length(deleted_sets) < num_to_remove
		removal_costs = worst_vertices(tour, dist)
		ind = pdf_select(removal_costs, power)
		set_to_delete = member[tour[ind]]

        # perform the deletion
        push!(deleted_sets, set_to_delete)
		splice!(tour, ind)
	end
    return deleted_sets
end

""" removing a single continuos segment of the tour of size num_remove """
function segment_removal!(tour::Array{Int64, 1}, num_to_remove::Int64, member::Array{Int64,1})
	i = rand(1:length(tour))
	deleted_sets = Array{Int}(undef, 0)
	while length(deleted_sets) < num_to_remove
		i > length(tour) && (i = 1)
		push!(deleted_sets, member[tour[i]])
		splice!(tour, i)
	end
	return deleted_sets
end

"""  pick a random vertex, and delete its closest neighbors  """
function distance_removal!(tour::Array{Int64,1}, dist::Array{Int64, 2},
							   num_to_remove::Int64, member::Array{Int64,1}, power::Float64)
    deleted_sets = Array{Int}(undef, 0)
    deleted_vertices = Array{Int}(undef, 0)

    seed_index = rand(1:length(tour))
    push!(deleted_sets, member[tour[seed_index]])
    push!(deleted_vertices, tour[seed_index])
    splice!(tour, seed_index)

    while length(deleted_sets) < num_to_remove
        # pick a random vertex from the set of deleted vertices
        seed_vertex = rand(deleted_vertices)
        # find closest vertex to the seed vertex that's still in the tour
        mindist = zeros(Int64, length(tour))
        for i = 1:length(tour)
			mindist[i] = min(dist[seed_vertex, tour[i]], dist[tour[i], seed_vertex])
        end
        del_index = pdf_select(mindist, power)
        push!(deleted_sets, member[tour[del_index]])
        push!(deleted_vertices, tour[del_index])
        splice!(tour, del_index)
    end

    return deleted_sets
end

"""
determine the cost of removing each vertex from the tour, given that all others remain.
"""
function worst_vertices(tour::Array{Int64, 1}, dist::Array{Int64, 2})
    removal_cost = zeros(Int64, length(tour))
    @inbounds for i = 1:length(tour)
        if i == 1
            removal_cost[i] = dist[tour[end], tour[i]] +
								dist[tour[i], tour[i+1]] - dist[tour[end], tour[i+1]]
        elseif i == length(tour)
            removal_cost[i] = dist[tour[i-1], tour[i]] +
								dist[tour[i], tour[1]] - dist[tour[i-1], tour[1]]
        else
            removal_cost[i] = dist[tour[i-1], tour[i]] +
								dist[tour[i], tour[i+1]] - dist[tour[i-1], tour[i+1]]
		end
    end
    return removal_cost
end


text:
push!(LOAD_PATH,joinpath("..","src"))
using Multilane
using Base.Test
using GenerativeModels
using POMDPs

import Iterators.product
import Base.assert

##############
##UNIT TESTS##
##############

function assert(expr,val,fn::Function= ==,varname::AbstractString="")
	if !fn(expr,val)
    error("Assertion failed: $varname : expected $val, got $expr")
	end
end

function assert(expr,val,flag::Bool)
	if expr == val
		error("Assertion failed: expected NOT $val, got $expr")
	end
end

function Multilane.MLState(pos::Real, vel::Real, cars::Array{CarState,1}, x::Real=50.)
    #x = mdp.phys_param.lane_length/2.
    insert!(cars,1,CarState(x, pos, vel, 0, Multilane.NORMAL, 0))
    return MLState(0.0, 0.0, cars, nothing)
end
# function MLState(crashed::Bool,pos::Real,vel::Real,cars::Array{CarState,1},x::Real=50.)
#     insert!(cars,1,CarState(x, pos, vel, 0, NORMAL, 0))
#     return MLState(false, crashed, 0.0, 0.0, cars)
# end
function Multilane.MLState(crashed::Bool,x::Float64,t::Float64,pos::Real,vel::Real,cars::Array{CarState,1},x_ego::Real=50.)
    insert!(cars,1,CarState(x_ego, pos, vel, 0, Multilane.NORMAL, 0))
    if crashed
        return MLState(x, t, cars, :crashed)
    else
        return MLState(x, t, cars)
    end
end

function test_hashing(s::AbstractString,ps)
	##TODO: make a more comprehensive test that uses more of the created objects
	println("\t\tTesting $s Hashing")
	if length(ps) < 2
		error("Make More objects for testing $s hashing")
	end
	d1 = Dict{typeof(ps[1]),Int}()
	d1[ps[1]] = 1
	d1[ps[1]] = 1
	assert(length(d1),1)
	d2 = Dict{typeof(ps[1]),Int}(p=>3 for p in [ps[2];ps[2]])
	assert(length(d2),1)
	assert(get(d1,ps[1],0),1)
end

function test_equality(s::AbstractString,ps)
	println("\t\tTesting $s Equality")
	if length(ps) < 3
		error("Make More objects for testing $s equality")
	end
	for i = 1:length(ps)
		for j = 1:length(ps)
			if i == j
				assert(ps[i],deepcopy(ps[j]))
			else
				assert(ps[i],ps[j],true)
				assert(ps[j],ps[i],true)
			end
		end
	end
	#assert(ps[1] == ps[1])
	#assert(ps[2] != ps[3])
	#assert(ps[3] != ps[2])
end

include("test_idm.jl")
include("test_mobil.jl")
include("test_types.jl")
include("test_crash.jl")
include("test_model.jl")
include("test_actions.jl")

println("Running Unit Tests...")
test_idm()
test_mobil()
test_mdp_types()
test_crash()
test_model()
test_actions()

# println("Running an OriginalMDP example...")
# include("run_example.jl")

include("change_into.jl")

println("Running a NoCrashMDP example...")
include("run_nocrash.jl")
# write_tmp_gif(mdp, sim)
# println("Running a NoCrashMDP example with the robust solver...")
# include("run_with_robust_solver.jl")

ext: jl

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
using Base.Test
using celerite

omega = 2pi/12.203317
alpha_r = [1.0428542, 0.38361831, 0.30345984/2, 0.30345984/2]
alpha_i = [0. ,               0. , 0.1, 0.1]
beta = [complex(0.1229159,0.0),complex(0.48922908,0.0),complex(0.09086397,omega),complex(0.09086397,-omega)]
nt = 512
#t = collect(linspace(0,nt-1,nt))
t = sort!(rand(nt))
acf = zeros(nt)
p = length(alpha_r)

# First solve in traditional manner:
w = 0.03027 * ones(nt)
A = zeros(Float64,nt,nt)
for i=1:nt
  A[i,i] += w[i]
  for j=1:nt
    for k=1:p
      A[i,j] += alpha_r[k] * exp(-real(beta[k])*abs(t[j]-t[i]))*cos(imag(beta[k])*(t[j]-t[i]))
      A[i,j] += alpha_i[k] * exp(-real(beta[k])*abs(t[j]-t[i]))*sin(abs(imag(beta[k])*(t[j]-t[i])))
    end
  end
end

# Compute a realization of correlated noise with this covariance matrix.
y=randn(nt)
# First, do Cholesky decomposition:
sqrta = chol(A)
# Now make a realization of the correlated noise:
corrnoise=*(transpose(sqrta),y);

# Now, solve for A*y = x to get inverse of Kernel (A) times the correlated noise:
x2 = \(A,corrnoise)

# Take dot product with correlated noise, and compare with original noise realization:
@test_approx_eq dot(y,y) dot(corrnoise,x2)
println("Dot product of white noise:                          ",dot(y,y))
println("Dot product of correlated noise with inverse kernel: ",dot(corrnoise,x2))

# Now use Ambikarasan O(N) method:
y = corrnoise
n = nt
alpha_final = [1.0428542, 0.38361831, 0.30345984]
alpha_imag =  [0. ,  0. ,  0.2]
beta_real_final = [0.1229159,0.48922908,0.09086397]
beta_imag_final = [0.0,0.0,omega]
w0 = 0.03027
p_final = 3
p0_final = 2
nex_final = (4(p_final-p0_final)+2p0_final+1)*(n-1)+1
m1_final = 2(p_final-p0_final)+p0_final+2
if p0_final == p_final
  m1_final = p0_final + 1
end
width_final = 2m1_final+1
aex_final = zeros(Float64,width_final,nex_final)
al_small_final = zeros(Float64,m1_final,nex_final)
indx_final= collect(1:nex_final)
tic()
logdeta_final = compile_matrix_symm(alpha_final,alpha_imag,beta_real_final,beta_imag_final,w0,t,nex_final,aex_final,al_small_final,indx_final)
tic()
bex = zeros(Float64,nex_final)
log_like_final = compute_likelihood(p_final,p0_final,y,aex_final,al_small_final,indx_final,logdeta_final,bex)

# Now select solution:
#x1 = zeros(n)
x3 = zeros(n)
log_like3 = 0.0
for i=1:n
  log_like3 += x2[i]*y[i]
end
log_like3 = -0.5 * log_like3
log_like3 += -0.5*logdet(A) - 0.5*n*log(2*pi)

@test_approx_eq logdeta_final logdet(A)
@test_approx_eq log_like_final log_like3
println("Log Determinant: ",logdeta_final," ",logdet(A))
println("Log Likelihood:  ",log_like_final," ",log_like3)


text:
module HTTP2

import HttpCommon: Headers

# package code goes here
include("Frame.jl")
include("Session.jl")

## Below we try to fire a HTTP2 request. This request is meant to be tested
## against a local nghttp2 server running on port 9000.
##
## The server is started by the command: `nghttpd --verbose --no-tls 9000`
import HTTP2.Session

function request(dest, port, url)
    ## Fire a HTTP connection to port 9000
    buffer = connect(dest, port)

    ## Create a HTTPConnection object
    connection = Session.new_connection(buffer; isclient=true)

    ## Create a request with headers
    headers = Headers(":method" => "GET",
                      ":path" => url,
                      ":scheme" => "http",
                      ":authority" => "127.0.0.1:9000",
                      "accept" => "*/*",
                      "accept-encoding" => "gzip, deflate",
                      "user-agent" => "HTTP2.jl")

    Session.put_act!(connection, Session.ActSendHeaders(Session.next_free_stream_identifier(connection), headers, true))

    return (Session.take_evt!(connection).headers, Session.take_evt!(connection).data)
end

function handle_util_frames_until(connection)
    received = Session.recv_next(connection)
    while typeof(received) == Frame.SettingsFrame || typeof(received) == Frame.PriorityFrame
        # do nothing for now
        received = Session.recv_next(connection)
    end
    return received
end

function serve(port, body)
    server = listen(port)

    println("Server started.")
    while(true)
        buffer = accept(server)
        println("Processing a connection ...")

        connection = Session.new_connection(buffer; isclient=false)
        ## Recv the client preface, and send an empty SETTING frame.

        headers_evt = Session.take_evt!(connection)
        stream_identifier = headers_evt.stream_identifier

        sending_headers = Headers(":status" => "200",
                                  "server" => "HTTP2.jl",
                                  "date" => Dates.format(now(Dates.UTC), Dates.RFC1123Format),
                                  "content-type" => "text/html; charset=UTF-8")

        Session.put_act!(connection, Session.ActSendHeaders(stream_identifier, sending_headers, false))
        Session.put_act!(connection, Session.ActSendData(stream_identifier, body, true))

        ## We are done!
    end
end

end # module


text:
import Printf: @printf
using Distributed
using FromFile
@from "Core.jl" import SRThreaded, SRSerial, SRDistributed

function debug(verbosity, string...)
    if verbosity > 0
        println(string...)
    end
end

function debug_inline(verbosity, string...)
    if verbosity > 0
        print(string...)
    end
end

function getTime()::Int
    return round(Int, 1e3*(time()-1.6e9))
end

function check_numeric(n)
    return tryparse(Float64, n) !== nothing
end

function is_anonymous_function(op)
	op_string = string(nameof(op))
	return length(op_string) > 1 && op_string[1] == '#' && check_numeric(op_string[2:2])
end

function recursive_merge(x::AbstractVector...)
    cat(x...; dims=1)
end

function recursive_merge(x::AbstractDict...)
    merge(recursive_merge, x...)
end

function recursive_merge(x...)
    x[end]
end

isgood(x::T) where {T<:Number} = !(isnan(x) || !isfinite(x))
isgood(x) = true
isbad(x) = !isgood(x)

macro return_on_false(flag, retval)
    :(if !$(esc(flag))
          return ($(esc(retval)), false)
    end)
end

function next_worker(worker_assignment::Dict{Tuple{Int,Int}, Int}, procs::Vector{Int})::Int
    job_counts = Dict(proc=>0 for proc in procs)
    for (key, value) in worker_assignment
        @assert haskey(job_counts, value)
        job_counts[value] += 1
    end
    least_busy_worker = reduce(
        (proc1, proc2) -> (
            job_counts[proc1] <= job_counts[proc2] ? proc1 : proc2
        ),
        procs
    )
    return least_busy_worker
end

function next_worker(worker_assignment::Dict{Tuple{Int,Int}, Int}, procs::Nothing)::Int
    return 0
end

macro sr_spawner(parallel, p, expr)
    quote
        if $(esc(parallel)) == SRSerial
            $(esc(expr))
        elseif $(esc(parallel)) == SRDistributed
            @spawnat($(esc(p)), $(esc(expr)))
        else
            Threads.@spawn($(esc(expr)))
        end
    end
end


text:
using ProximalAlgorithms
using ProximalOperators
using LinearAlgebra
using Random
using Test

@testset "Nonconvex QP (tiny, $T)" for T in [Float64]
    Q = Matrix(Diagonal(T[-0.5, 1.0]))
    q = T[0.3, 0.5]
    low = T(-1.0)
    upp = T(+1.0)

    f = Quadratic(Q, q)
    g = IndBox(low, upp)

    n = 2

    Lip = maximum(diag(Q))
    gamma = T(0.95) / Lip

    @testset "PANOC" begin
        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.PANOC()
        x, it = solver(x0, f = f, g = g)
        z = min.(upp, max.(low, x .- gamma .* (Q * x + q)))
        @test norm(x - z, Inf) / gamma <= solver.tol
        @test x0 == x0_backup
    end

    @testset "PANOCplus" begin
        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.PANOCplus()
        x, it = solver(x0, f = f, g = g)
        z = min.(upp, max.(low, x .- gamma .* (Q * x + q)))
        @test norm(x - z, Inf) / gamma <= solver.tol
        @test x0 == x0_backup
    end

    @testset "ZeroFPR" begin
        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.ZeroFPR()
        x, it = solver(x0, f = f, g = g)
        z = min.(upp, max.(low, x .- gamma .* (Q * x + q)))
        @test norm(x - z, Inf) / gamma <= solver.tol
        @test x0 == x0_backup
    end

    @testset "LiLin" begin
        x0 = zeros(T, n)
        x0_backup = copy(x0)
        solver = ProximalAlgorithms.LiLin(gamma = gamma)
        x, it = solver(x0, f = f, g = g)
        z = min.(upp, max.(low, x .- gamma .* (Q * x + q)))
        @test norm(x - z, Inf) / gamma <= solver.tol
        @test x0 == x0_backup
    end
end

@testset "Nonconvex QP (small, $T)" for T in [Float64]
    @testset "Random problem $k" for k = 1:5
        Random.seed!(k)

        n = 100
        A = randn(T, n, n)
        U, R = qr(A)
        eigenvalues = T(2) .* rand(T, n) .- T(1)
        Q = U * Diagonal(eigenvalues) * U'
        Q = 0.5 * (Q + Q')
        q = randn(T, n)

        low = T(-1.0)
        upp = T(+1.0)

        f = Quadratic(Q, q)
        g = IndBox(low, upp)

        Lip = maximum(abs.(eigenvalues))
        gamma = T(0.95) / Lip

        TOL = 1e-4

        @testset "PANOC" begin
            x0 = zeros(T, n)
            x0_backup = copy(x0)
            solver = ProximalAlgorithms.PANOC(tol = TOL)
            x, it = solver(x0, f = f, g = g)
            z = min.(upp, max.(low, x .- gamma .* (Q * x + q)))
            @test norm(x - z, Inf) / gamma <= solver.tol
            @test x0 == x0_backup
        end

        @testset "PANOCplus" begin
            x0 = zeros(T, n)
            x0_backup = copy(x0)
            solver = ProximalAlgorithms.PANOCplus(tol = TOL)
            x, it = solver(x0, f = f, g = g)
            z = min.(upp, max.(low, x .- gamma .* (Q * x + q)))
            @test norm(x - z, Inf) / gamma <= solver.tol
            @test x0 == x0_backup
        end

        @testset "ZeroFPR" begin
            x0 = zeros(T, n)
            x0_backup = copy(x0)
            solver = ProximalAlgorithms.ZeroFPR(tol = TOL)
            x, it = solver(x0, f = f, g = g)
            z = min.(upp, max.(low, x .- gamma .* (Q * x + q)))
            @test norm(x - z, Inf) / gamma <= solver.tol
            @test x0 == x0_backup
        end

        @testset "LiLin" begin
            x0 = zeros(T, n)
            x0_backup = copy(x0)
            solver = ProximalAlgorithms.LiLin(gamma = gamma, tol = TOL)
            x, it = solver(x0, f = f, g = g)
            z = min.(upp, max.(low, x .- gamma .* (Q * x + q)))
            @test norm(x - z, Inf) / gamma <= solver.tol
            @test x0 == x0_backup
        end
    end
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-DEV
#     language: julia
#     name: julia-1.7
# ---

# %%
using BenchmarkTools
using DistributedArrays
struct Foo{T, N} d::DArray{T, N} end
struct Bar{T, N, A} d::DArray{T, N, A} end

# %%
d = drand(4, 3, 2)
l = localpart(d)
foo = Foo(d)
bar = Bar(d);

# %%
f(foobar, i, j, k) = localpart(foobar.d)[i, j, k]
g(d, i, j, k) = localpart(d)[i, j, k]

# %%
@btime f($foo, 3, 2, 1)

# %%
@btime f($bar, 3, 2, 1)

# %%
@btime g($d, 3, 2, 1)

# %%
@btime getindex($l, 3, 2, 1)

# %%
@code_warntype f(foo, 3, 2, 1)

# %%
@code_warntype f(bar, 3, 2, 1)

# %%
isconcretetype(Array{Float64, 3})

# %%
isconcretetype(DArray{Float64, 3})

# %%
typeof(d)

# %%
isconcretetype(DArray{Float64, 3, Array{Float64, 3}})

# %%

max_stars_repo_path: 0004/DArray{Float64, 3} is not a concrete type.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using LogCompose, Test, Logging, SyslogLogging, SyslogLoggingCompose

function test()
    config = joinpath(@__DIR__, "testapp.toml")

    let logger = LogCompose.logger(config, "syslog.testapp"; section="loggers")
        with_logger(logger) do
            @info("testsyslog")
        end
    end
end

test()


text:
"""
    activate(; shell=true)
    activate(config::Config, args...; shell=true)
    activate(env::Environment; shell=true)

Modifies the current environment to operate within a specific playground environment.
When `shell=true` a new shell environment will be created.
However, when `shell=false` the existing julia REPL will be modifed and
`deactivate()` must be called to restore the REPL state.
"""
activate(; shell=true) = activate(Environment(); shell=shell)

function activate(config::Config, args...; shell=true)
    activate(Environment(config, args...); shell=shell)
end

function activate(env::Environment; shell=true)
    debug(logger, "Activating playground $(name(env))...")

    if shell
        withenv(env) do
            sh = getshell(env.config.default_shell)
            run(sh, env)
        end
    else
        repl = Julia()
        prompt = getprompt(repl, env)

        old = Dict{Symbol, Any}()
        old[:ENV] = set!(env, getenvs(env)...)
        try
            old[:PROMPT] = input_prompt()
            input_prompt!(prompt, :magenta)
        catch e
            warn(logger, "Failed to set the julia prompt to $prompt ($e)")
        finally
            push!(cache, old)
        end
    end
end

"""
    deactivate()

Deactivates the active environment and restores the original julia environment.
"""
function deactivate()
    if !isempty(cache)
        debug(logger, "Deactivating playground ...")
        old = pop!(cache)

        try
            input_prompt!(old[:PROMPT])
        catch _
            warn(logger, string("Failed to restore the julia prompt."))
        end

        restore!(old[:ENV])
    else
        warn(logger, "There are no cached environment settings to restore.")
    end
end


text:
using Test
using LinearAlgebra
using DataFlowTasks
using DataFlowTasks: DAG, num_nodes, num_edges, addnode!, addedge!,
                    has_edge, outneighbors, inneighbors, isconnected,
                    adjacency_matrix, addedge_transitive!, remove_node!, R,W,RW, TaskGraph

@testset "DAG{Int}" begin
    dag = DAG{Int}()
    @test num_nodes(dag) == 0
    addnode!(dag,1)
    @test num_nodes(dag) == 1
    @test num_edges(dag) == 0
    addnode!(dag, 2)
    @test num_nodes(dag) == 2
    addedge!(dag, 1, 2)
    @test num_edges(dag) == 1
    @test has_edge(dag, 1, 2) == true
    @test has_edge(dag, 2, 1) == false
    @test outneighbors(dag, 1) == Set(2)
    @test inneighbors(dag, 1) == Set()
    @test outneighbors(dag, 2) == Set()
    @test inneighbors(dag, 2) == Set(1)
    addnode!(dag, 3)
    addnode!(dag, 4)
    addedge!(dag,2,4)
    @test isconnected(dag, 1, 3) == false
    @test isconnected(dag, 2, 4) == true
    @test isconnected(dag, 1, 4) == true # 1 --> 2 --> 4
    @test adjacency_matrix(dag) == [0 1 0 0;
                                    0 0 0 1;
                                    0 0 0 0;
                                    0 0 0 0]
    @test num_edges(addedge_transitive!(dag, 1, 4)) == 2 # vtx 1 is already connected to vtx through vtx2
end

@testset "TaskGraph" begin
    dag = TaskGraph()
    C = rand(3,3)
    cl1 = @dtask rmul!(C,2) (C,) (RW,)
    @test num_nodes(dag) == 0
    addnode!(dag,cl1)
    @test num_nodes(dag) == 1
    @test num_edges(dag) == 0
    cl2 = @dtask rmul!(C,2) (C,) (RW,)
    addnode!(dag, cl2)
    @test num_nodes(dag) == 2
    addedge!(dag, cl1, cl2)
    @test num_edges(dag) == 1
    @test has_edge(dag, cl1, cl2) == true
    @test has_edge(dag, cl2, cl1) == false
    @test outneighbors(dag, cl1) == Set((cl2,))
    @test inneighbors(dag, cl1) == Set()
    @test outneighbors(dag, cl2) == Set()
    @test inneighbors(dag, cl2) == Set([cl1])
    cl3 = @dtask rmul!(C,2) (C,) (RW,)
    cl4 = @dtask rmul!(C,2) (C,) (RW,)
    addnode!(dag, cl3)
    addnode!(dag, cl4)
    addedge!(dag,cl2,cl4)
    @test isconnected(dag, cl1, cl3) == false
    @test isconnected(dag, cl2, cl4) == true
    @test isconnected(dag, cl1, cl4) == true # 1 --> 2 --> 4
    @test adjacency_matrix(dag) == [0 1 0 0;
                                    0 0 0 1;
                                    0 0 0 0;
                                    0 0 0 0]
    @test num_edges(addedge_transitive!(dag, cl1, cl4)) == 2 # vtx 1 is already connected to vtx through vtx2
end

@testset "Buffering" begin
    sz = 10
    dag = DAG{Int}(sz)
    for i in 1:sz
        addnode!(dag, i)
    end
    @test num_nodes(dag) == sz
    t = @async addnode!(dag, 11)
    sleep(0.1)
    @test num_nodes(dag) == sz
    remove_node!(dag,1)
    sleep(0.1)
    @test num_nodes(dag) == sz
end


text:
type GaussianProcess
    noise::Float64
    kernel::CovarianceKernel
end

type GaussianProcessPrediction
    mean::Vector{Float64}
    sd::Vector{Float64}
end

function predict(gp::GaussianProcess,
                 x_in::Vector{Float64},
                 y_in::Vector{Float64},
                 x_out::Vector{Float64})

    n_in = size(x_in, 1)
    n_out = size(x_out, 1)

    mu_in = zeros(n_in)
    mu_out = zeros(n_out)

    sigma_in = gpcov(gp.kernel, x_in, x_in)
    sigma_in += gp.noise ^ 2 * eye(n_in)

    sigma_out = gpcov(gp.kernel, x_out, x_out)

    sigma_in_out = gpcov(gp.kernel, x_in, x_out)
    sigma_out_in = sigma_in_out'

    pred_mean = mu_out
    pred_mean += sigma_out_in * inv(sigma_in) * y_in

    pred_var = sigma_out
    pred_var -= sigma_out_in * inv(sigma_in) * sigma_in_out

    pred_sd = sqrt(diag(pred_var))

    GaussianProcessPrediction(pred_mean, pred_sd)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# BEGIN 1.0 deprecations

# END 1.0 deprecations

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# this file setup is less than ideal. We will refactor this with packing macros
# to automatically generate the converters and packed data types.
# The eval functions are also currently spread out, they should be concentrated here as well.
# The Pose2D and Pose3D types will most likely be packaged with the RoME package in the future.

# Examples costraint functions which can be used, however,
# \/\/\/these won't pack into ProtoBuf -- need to write special converters
# for Neo4j DataBase storage of complicated types. See constraints hereafter
# for more standard types.
type UniPriorPose2D <: Singleton
  Z::Distributions.MvNormal
end
type GMMPriorPose2D <: Singleton
  Z::Array{Distributions.MvNormal,1}
  W::Array{Float64,1}
end
type KDEPriorPoint2D <: Singleton
  Z::BallTreeDensity
end
type KDERangePoint2D <: Pairwise
  Z::BallTreeDensity
end
# ^^^will work on these later

# Active constraint types listed below
# -------------

# define the pose group
type Odo <: Pairwise
    Zij::Array{Float64,2} # 0rotations, 1translation in each column
    Cov::Array{Float64,2}
    W::Array{Float64,1}
    Odo() = new()
    Odo(x...) = new(x[1], x[2], x[3])
end
type PackedOdo
    vecZij::Array{Float64,1} # 0rotations, 1translation in each column
    dimz::Int64
    vecCov::Array{Float64,1}
    dimc::Int64
    W::Array{Float64,1}
    PackedOdo() = new()
    PackedOdo(x...) = new(x[1], x[2], x[3], x[4], x[5])
end
function convert(::Type{Odo}, d::PackedOdo)
  Zij = reshapeVec2Mat(d.vecZij,d.dimz)
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Odo(Zij, Cov, d.W)
end
function convert(::Type{PackedOdo}, d::Odo)
  v1 = d.Zij[:];
  v2 = d.Cov[:];
  return PackedOdo(v1,length(v1),
                    v2,length(v2),
                    d.W)
end
function convert(::Type{FunctionNodeData{PackedOdo}}, d::FunctionNodeData{Odo})
  return FunctionNodeData{PackedOdo}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(PackedOdo, d.fnc))
end
function convert(::Type{FunctionNodeData{Odo}}, d::FunctionNodeData{PackedOdo})
  return FunctionNodeData{Odo}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(Odo, d.fnc))
end
function FNDencode(d::FunctionNodeData{Odo})
  return convert(FunctionNodeData{PackedOdo}, d)
end
function FNDdecode(d::FunctionNodeData{PackedOdo})
  return convert(FunctionNodeData{Odo}, d)
end

type OdoMM <: Pairwise
    Zij::Array{Float64,2} # 0rotations, 1translation in each column
    Cov::Array{Float64,2}
    W::Array{Float64,1}
end

type Ranged <: Pairwise
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    Ranged() = new()
    Ranged(x...) = new(x[1], x[2], x[3])
end
function FNDencode(d::FunctionNodeData{Ranged})
  return d
end
function FNDdecode(d::FunctionNodeData{Ranged})
  return d
end

type GenericMarginal <: Pairwise
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    GenericMarginal() = new()
    GenericMarginal(a,b,c) = new(a,b,c)
end
function FNDencode(d::FunctionNodeData{GenericMarginal})
  return d
end
function FNDdecode(d::FunctionNodeData{GenericMarginal})
  return d
end

# ------------------------------------------------------------

type Obsv2 <: Singleton
    pts::Array{Float64,2}
    bws::Array{Float64,2}
    W::Array{Float64,1}
    Obsv2() = new()
    Obsv2(x...) = new(x[1], x[2], x[3])
end
type PackedObsv2
    vecZij::Array{Float64,1} # 0rotations, 1translation in each column
    dimz::Int64
    vecCov::Array{Float64,1}
    dimc::Int64
    W::Array{Float64,1}
    PackedObsv2() = new()
    PackedObsv2(x...) = new(x[1],x[2],x[3],x[4],x[5])
end
function convert(::Type{Obsv2}, d::PackedObsv2)
  Zij = reshapeVec2Mat(d.vecZij,d.dimz)
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Obsv2(Zij, Cov, d.W)
end
function convert(::Type{PackedObsv2}, d::Obsv2)
  v1 = d.pts[:];
  v2 = d.bws[:];
  return PackedObsv2(v1,size(d.pts,1),
                    v2,size(d.bws,1),
                    d.W)
end
function convert(::Type{FunctionNodeData{PackedObsv2}}, d::FunctionNodeData{Obsv2})
  return FunctionNodeData{PackedObsv2}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(PackedObsv2, d.fnc))
end
function convert(::Type{FunctionNodeData{Obsv2}}, d::FunctionNodeData{PackedObsv2})
  return FunctionNodeData{Obsv2}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(Obsv2, d.fnc))
end
function FNDencode(d::FunctionNodeData{Obsv2})
  return convert(FunctionNodeData{PackedObsv2}, d)
end
function FNDdecode(d::FunctionNodeData{PackedObsv2})
  return convert(FunctionNodeData{Obsv2}, d)
end

# ------------------------------------------------------

type PriorPose2 <: Singleton
    Zi::Array{Float64,2}
    Cov::Array{Float64,2}
    W::Array{Float64,1}
    PriorPose2() = new()
    PriorPose2(x...) = new(x[1], x[2], x[3])
end
type PackedPriorPose2
    vecZij::Array{Float64,1} # 0rotations, 1translation in each column
    dimz::Int64
    vecCov::Array{Float64,1}
    dimc::Int64
    W::Array{Float64,1}
    PackedPriorPose2() = new()
    PackedPriorPose2(x...) = new(x[1], x[2], x[3], x[4], x[5])
end
function convert(::Type{PriorPose2}, d::PackedPriorPose2)
  Zij = reshapeVec2Mat(d.vecZij,d.dimz)
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return PriorPose2(Zij, Cov, d.W)
end
function convert(::Type{PackedPriorPose2}, d::PriorPose2)
  v1 = d.Zij[:];
  v2 = d.Cov[:];
  return PackedPriorPose2(v1,size(d.Zij,1),
                          v2,size(d.Cov,1),
                          d.W)
end
function convert(::Type{FunctionNodeData{PackedPriorPose2}}, d::FunctionNodeData{PriorPose2})
  return FunctionNodeData{PackedPriorPose2}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(PackedPriorPose2, d.fnc))
end
function convert(::Type{FunctionNodeData{PriorPose2}}, d::FunctionNodeData{PackedPriorPose2})
  return FunctionNodeData{PriorPose2}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(PriorPose2, d.fnc))
end
function FNDencode(d::FunctionNodeData{PriorPose2})
  return convert(FunctionNodeData{PackedPriorPose2}, d)
end
function FNDdecode(d::FunctionNodeData{PackedPriorPose2})
  return convert(FunctionNodeData{PriorPose2}, d)
end

# ------------------------------------

type Pose2Pose2 <: Pairwise
    Zij::Array{Float64,2} # 2translations, 1rotation
    Cov::Array{Float64,2}
    W::Array{Float64,1}
    Pose2Pose2() = new()
    Pose2Pose2(x...) = new(x[1], x[2], x[3])
end
type PackedPose2Pose2
  vecZij::Array{Float64,1} # 2translations, 1rotation
  dimz::Int64
  vecCov::Array{Float64,1}
  dimc::Int64
  W::Array{Float64,1}
  PackedPose2Pose2() = new()
  PackedPose2Pose2(x...) = new(x[1], x[2], x[3], x[4], x[5])
end
function convert(::Type{Pose2Pose2}, d::PackedPose2Pose2)
  return Pose2Pose2(reshapeVec2Mat(d.vecZij,d.dimz),
                    reshapeVec2Mat(d.vecCov, d.dimc), W)
end
function convert(::Type{PackedPose2Pose2}, d::Pose2Pose2)
  v1 = d.Zij[:];
  v2 = d.Cov[:];
  return PackedPose2Pose2(v1,size(d.Zij,1),
                          v2,size(d.Cov,1),
                          d.W)
end
function convert(::Type{FunctionNodeData{PackedPose2Pose2}}, d::FunctionNodeData{Pose2Pose2})
  return FunctionNodeData{PackedPose2Pose2}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(PackedPose2Pose2, d.fnc))
end
function convert(::Type{FunctionNodeData{Pose2Pose2}}, d::FunctionNodeData{PackedPose2Pose2})
  return FunctionNodeData{Pose2Pose2}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(Pose2Pose2, d.fnc))
end
function FNDencode(d::FunctionNodeData{Pose2Pose2})
  return convert(FunctionNodeData{PackedPose2Pose2}, d)
end
function FNDdecode(d::FunctionNodeData{PackedPose2Pose2})
  return convert(FunctionNodeData{Pose2Pose2}, d)
end

# --------------------------------------------

type Pose2DPoint2DBearingRange <: Pairwise
    Zij::Array{Float64,2} # bearing and range hypotheses as columns
    Cov::Array{Float64,2}
    W::Array{Float64,1}
    Pose2DPoint2DBearingRange() = new()
    Pose2DPoint2DBearingRange(x...) = new(x[1],x[2],x[3])
end
type PackedPose2DPoint2DBearingRange
    vecZij::Array{Float64,1} # 0rotations, 1translation in each column
    dimz::Int64
    vecCov::Array{Float64,1}
    dimc::Int64
    W::Array{Float64,1}
    PackedPose2DPoint2DBearingRange() = new()
    PackedPose2DPoint2DBearingRange(x...) = new(x[1], x[2], x[3], x[4], x[5])
end
function convert(::Type{Pose2DPoint2DBearingRange}, d::PackedPose2DPoint2DBearingRange)
  Zij = reshapeVec2Mat(d.vecZij,d.dimz)
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Pose2DPoint2DBearingRange(Zij, Cov, d.W)
end
function convert(::Type{PackedPose2DPoint2DBearingRange}, d::Pose2DPoint2DBearingRange)
  v1 = d.Zij[:];
  v2 = d.Cov[:];
  return PackedPose2DPoint2DBearingRange(v1,size(d.Zij,1),
                                         v2,size(d.Cov,1),
                                         d.W)
end
function convert(::Type{FunctionNodeData{PackedPose2DPoint2DBearingRange}}, d::FunctionNodeData{Pose2DPoint2DBearingRange})
  return FunctionNodeData{PackedPose2DPoint2DBearingRange}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(PackedPose2DPoint2DBearingRange, d.fnc))
end
function convert(::Type{FunctionNodeData{Pose2DPoint2DBearingRange}}, d::FunctionNodeData{PackedPose2DPoint2DBearingRange})
  return FunctionNodeData{Pose2DPoint2DBearingRange}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          convert(Pose2DPoint2DBearingRange, d.fnc))
end
function FNDencode(d::FunctionNodeData{Pose2DPoint2DBearingRange})
  return convert(FunctionNodeData{PackedPose2DPoint2DBearingRange}, d)
end
function FNDdecode(d::FunctionNodeData{PackedPose2DPoint2DBearingRange})
  return convert(FunctionNodeData{Pose2DPoint2DBearingRange}, d)
end

# ------------------------------------------------------

type PriorPose3 <: Singleton
    Zi::SE3
    Cov::Array{Float64,2}
    PriorPose3() = new()
    PriorPose3(st::FloatInt, sr::Float64) = new(SE3(0), [[st*eye(3);zeros(3,3)];[zeros(3);sr*eye(3)]])
    PriorPose3(s::SE3, c::Array{Float64,2}) = new(s,c)
end
# type PackedPriorPose3
#     vecZij::Array{Float64,1} # 0rotations, 1translation in each column
#     dimz::Int64
#     vecCov::Array{Float64,1}
#     dimc::Int64
#     W::Array{Float64,1}
#     PackedPriorPose3() = new()
#     PackedPriorPose3(x...) = new(x[1], x[2], x[3], x[4], x[5])
# end
# function convert(::Type{PriorPose3}, d::PackedPriorPose3)
#   Zij = reshapeVec2Mat(d.vecZij,d.dimz)
#   Cov = reshapeVec2Mat(d.vecCov, d.dimc)
#   return PriorPose3(Zij, Cov, d.W)
# end
# function convert(::Type{PackedPriorPose3}, d::PriorPose3)
#   v1 = d.Zij[:];
#   v2 = d.Cov[:];
#   return PackedPriorPose3(v1,size(d.Zij,1),
#                           v2,size(d.Cov,1),
#                           d.W)
# end
# function convert(::Type{FunctionNodeData{PackedPriorPose3}}, d::FunctionNodeData{PriorPose3})
#   return FunctionNodeData{PackedPriorPose3}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           convert(PackedPriorPose3, d.fnc))
# end
# function convert(::Type{FunctionNodeData{PriorPose3}}, d::FunctionNodeData{PackedPriorPose3})
#   return FunctionNodeData{PriorPose3}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           convert(PriorPose3, d.fnc))
# end
# function FNDencode(d::FunctionNodeData{PriorPose3})
#   return convert(FunctionNodeData{PackedPriorPose3}, d)
# end
# function FNDdecode(d::FunctionNodeData{PackedPriorPose3})
#   return convert(FunctionNodeData{PriorPose3}, d)
# end

# ------------------------------------

type Pose3Pose3 <: Pairwise
    Zij::SE3 # 3translations, 3exponential param rotation
    Cov::Array{Float64,2}
    Pose3Pose3() = new()
    Pose3Pose3(st::FloatInt, sr::Float64) = new(SE3(0), [[st*eye(3);zeros(3,3)];[zeros(3);sr*eye(3)]])
    Pose3Pose3(s::SE3, c::Array{Float64,2}) = new(s,c)
end
# type PackedPose3Pose3
#   vecZij::Array{Float64,1} # 2translations, 1rotation
#   dimz::Int64
#   vecCov::Array{Float64,1}
#   dimc::Int64
#   W::Array{Float64,1}
#   PackedPose2Pose2() = new()
#   PackedPose2Pose2(x...) = new(x[1], x[2], x[3], x[4], x[5])
# end
# function convert(::Type{Pose3Pose3}, d::PackedPose3Pose3)
#   return Pose3Pose3(reshapeVec2Mat(d.vecZij,d.dimz),
#                     reshapeVec2Mat(d.vecCov, d.dimc), W)
# end
# function convert(::Type{PackedPose3Pose3}, d::Pose3Pose3)
#   v1 = d.Zij[:];
#   v2 = d.Cov[:];
#   return PackedPose3Pose3(v1,size(d.Zij,1),
#                           v2,size(d.Cov,1),
#                           d.W)
# end
# function convert(::Type{FunctionNodeData{PackedPose3Pose3}}, d::FunctionNodeData{Pose3Pose3})
#   return FunctionNodeData{PackedPose3Pose3}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           convert(PackedPose3Pose3, d.fnc))
# end
# function convert(::Type{FunctionNodeData{Pose3Pose3}}, d::FunctionNodeData{PackedPose3Pose3})
#   return FunctionNodeData{Pose3Pose3}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           convert(Pose3Pose3, d.fnc))
# end
# function FNDencode(d::FunctionNodeData{Pose3Pose3})
#   return convert(FunctionNodeData{PackedPose3Pose3}, d)
# end
# function FNDdecode(d::FunctionNodeData{PackedPose3Pose3})
#   return convert(FunctionNodeData{Pose3Pose3}, d)
# end

# -----------------------


text:
#Merge methods together using metaprogramming?
function xgridlines!(ec::EChart; show::Bool = true,
								 interval::Union{Int, String, JSFunction, Nothing} = "auto",
								 kwargs...)

	#Set if not present
    ec.xAxis[1].splitLine == nothing ? ec.xAxis[1].splitLine = SplitLine() : nothing

	#Set top-level properties
    ec.xAxis[1].splitLine.show = show
	ec.xAxis[1].splitLine.interval = interval

	#kwargs all the lineStyle arguments
	length(kwargs) > 0 ? for (k, v) in kwargs setfield!(ec.xAxis[1].splitLine.lineStyle, k, v) end : nothing

    return ec

end

function ygridlines!(ec::EChart; show::Bool = true,
								 interval::Union{Int, String, JSFunction, Nothing} = "auto",
								 kwargs...)

	#Set if not present
    ec.yAxis[1].splitLine == nothing ? ec.yAxis[1].splitLine = SplitLine() : nothing

	#Set top-level properties
    ec.yAxis[1].splitLine.show = show
	ec.yAxis[1].splitLine.interval = interval

	#kwargs all the lineStyle arguments
	length(kwargs) > 0 ? for (k, v) in kwargs setfield!(ec.yAxis[1].splitLine.lineStyle, k, v) end : nothing

    return ec

end


text:
"""
    sparsity_pattern(M::AbstractArray, ratio::Real, p::Real, min_per_row=0::Integer, min_per_col=0::Integer)

Compute a pattern matrix `M_pat` which is a p-norm sparsity pattern for `M`.

`M_pat` is a sparse `SparseMatrixCSC{Int64, Int64}` of the same shape as `M` and
contains only 0 or 1.

`ratio` is a measure for how sparse the pattern is supposed to be and must be in ``[0,1]``.
`p` determines which norm is used and should be in ``(0, \\infty ]``. `min_per_row` is the minimum number
of non-zeros needed per row and `min_per_col` is the minimum number of non-zeros needed per column,
Defaults are `0`.

See also: [`sparsify`](@ref), [`bin_sparse_matrix!`](@ref).

# Examples
```julia-repl
julia> sparsity_pattern([4 1; 0.1 17], 0.6, 2)
2×2 SparseMatrixCSC{Int64, Int64} with 2 stored entries:
 1  ⋅
 ⋅  1
```
"""
function sparsity_pattern(M::AbstractArray, ratio::Real, p::Real, min_per_row=0::Integer, min_per_col=0::Integer)
    0 <= ratio <= 1 || throw(DomainError(ratio, "ratio must be in [0,1]"))
    p > 0 || throw(DomainError(p, "p must be greater than zero"))
    min_per_row >= 0|| throw(DomainError("min_per_row must be greater equal zero"))
    min_per_col >= 0|| throw(DomainError("min_per_col must be greater equal zero"))

    # todo make almost symmetric matrices symmetric
    m = size(M, 1)
    n = size(M, 2)

    M_pat = spzeros(Int, m, n)
    @views for i in 1:m
        M_pat[i, :] = p_norm_sparsity_vector(M[i, :], ratio, p, min_per_row)
    end
    @views for j in 1:n
        M_pat[:, j] = max.(p_norm_sparsity_vector(M[:, j], ratio, p, min_per_col), M_pat[:, j])
    end
    return M_pat
end

function p_norm_sparsity_vector(v::AbstractVector, ratio::Real, p::Real, min_num_nnz::Integer)
    0 <= ratio <= 1 || throw(DomainError(ratio, "ratio must be in [0,1]"))
    p > 0 || throw(DomainError(p, "p must be greater than zero"))
    min_num_nnz >= 0|| throw(DomainError("min_num_nnz must be greater equal zero"))

    n = length(v)

    v = abs.(v)

    nnz_v = sum(v .> 0)

    if nnz_v == 0
        v_pat = spzeros(Int, n)
        return v_pat
    end

    nnz_v >= min_num_nnz || throw(ArgumentError("minimum number for non zeros entries bigger than non zero entries in v"))

    idx = sortperm(v)
    v = v[idx]

    if p == Inf
        if min_num_nnz == nnz_v
            num_needed = min_num_nnz
        else
            discarded_bound = (1-ratio)*v[end]
            num_needed = n - upper_bound(v[1:end-min_num_nnz], discarded_bound) + 1
        end
    else
        v = v.^p
        csum = cumsum(v)
        sum_threshold = csum[end] * (1 - ratio) ^ max(1, p)

        ub = upper_bound(csum, sum_threshold)
        num_needed = max(min_num_nnz, n - ub + 1)
    end

    data = ones(Int, num_needed)
    v_pat = sparsevec(idx[n - num_needed + 1:n], data, n)
    return v_pat
end

function upper_bound(v::AbstractVector, value::Real)
    n = length(v)

    if n == 0
        low = 0
    else
        low = 1
        high = n
        while low <= high
            mid = floor(Int, (low + high) / 2)
            if value < v[mid]
                high = mid - 1
            else
                low = mid + 1
            end
        end
    end
    return low
end


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

using Markdown
using InteractiveUtils

# ╔═╡ eb033006-f788-11ea-086f-357e2353381f
using Pkg, DrWatson

# ╔═╡ eb03721e-f788-11ea-10b4-07a4a4032d76
begin
	@quickactivate "SR2TuringPluto"
	using StatisticalRethinking
end

# ╔═╡ c4b79432-f788-11ea-1092-df50f9dca4bf
md"## Fig3.5t.jl"

# ╔═╡ eb03e064-f788-11ea-274c-cbe1c1978108
begin
  N = 100000
  d = rand(Binomial(9, 0.7), N);
  histogram(d; normalize=:probability,
    bins=-0.5:1:9.5, leg=false, xticks=0:9, bar_width=0.2)
end

# ╔═╡ eb0aecd8-f788-11ea-14da-fdfc08b9a986
md"## End of Fig3.5t.jl"

# ╔═╡ Cell order:
# ╟─c4b79432-f788-11ea-1092-df50f9dca4bf
# ╠═eb033006-f788-11ea-086f-357e2353381f
# ╠═eb03721e-f788-11ea-10b4-07a4a4032d76
# ╠═eb03e064-f788-11ea-274c-cbe1c1978108
# ╟─eb0aecd8-f788-11ea-14da-fdfc08b9a986


text:
# using QMTK

# op = Linear(3, 4)
# display(op)


text:
module Getopts

export getopts

done(instr, idx) = iterate(instr, idx)===nothing

"""
```
function skipcwhile(instr, startidx=1, mchars=" \\t")
```
Consume and skip characters **WHILE** they match characters specified in
mchars. Return index of first character following consumed characters
in string.
"""
function skipcwhile(instr, startidx=1, mchars=" \t")
  idx = startidx
  while !done(instr, idx) && instr[idx] in mchars
    idx = iterate(instr, idx)[2]
  end
  idx
end

"""
```
function getcuntil(instr, startidx=1, mchars=" \\t")
```
Consume and return characters in input string **UNTIL** a character
specified in mchars (match characters) occurs. Also return index
of first character following consumed characters in string.
"""
function getcuntil(instr, startidx=1, mchars=" \t")
  idx = startidx
  outstr = ""
  while !done(instr, idx) && !(instr[idx] in mchars)
    outstr = string(outstr, instr[idx])
    idx = iterate(instr, idx)[2]
  end
  outstr, idx
end

"""
```
function getcuntilunless(instr, startidx=1, mchars=" \\t", uchars="-")
```
Consume and return characters in input string **UNTIL** a character
specified in mchars (match characters) occurs, **UNLESS** the input string
starts with a character specified in uchars. Also return index
of first character following consumed characters in string.
"""
function getcuntilunless(instr, startidx=1, mchars=" \t", uchars="-")
  idx = skipcwhile(instr, startidx, mchars)
  outstr, idx = done(instr, idx) || (instr[idx] in uchars) ? ("", idx) : getcuntil(instr, idx, " \t")
end

getopt = getcuntil
getarg = getcuntilunless

function getopts(instr::AbstractString)
  opts = Dict{AbstractString,Array{AbstractString,1}}()
  argv = Array{AbstractString,1}()
  idx = 1
  while !done(instr, idx)
    while true
      arg, idx = getarg(instr, idx, " \t", "-")
      arg != "" && push!(argv, arg)
      (arg == "" || done(instr, idx)) && break
    end
    if !done(instr, idx)
      key, idx = getopt(instr, idx, " \t")
      !(key in keys(opts)) && (opts[key] = Array{AbstractString,1}())
      val, idx = getarg(instr, idx, " \t", "-")
      push!(opts[key], val)
    end
  end
  opts, argv
end

function getopts(inarray::AbstractArray{S,1}) where S<:AbstractString
  getopts(foldl(string,map(x->string(" ",x),inarray)))
end

end

max_stars_repo_name: msravi/Getopts

max_issues_repo_name: msravi/Getopts

max_forks_repo_name: msravi/Getopts


text:
typeof(2) #> Int64
typeof(Int64) #> DataType
typeof(Complex{Int64}) #> DataType
typeof(Any) #> DataType
typeof(DataType) #> DataType

super(Int64) #> Signed
super(Signed) #> Integer
super(Integer) #> Real
super(Real) #> Number
super(Number) #> Any
super(Any) #> Any

subtypes(Integer) #> 5-element Array{Any,1}:
# BigInt
# Bool
# Char
# Signed
# Unsigned
subtypes(Signed) #> 5-element Array{Any,1}:
# Int128
# Int16
# Int32
# Int64
# Int8
subtypes(Int64) #> 0-element Array{Any,1}
typeof(None) #> UnionType

Bool <: Integer #> true
Bool <: Any #> true
Bool <: Char #> false
issubtype(Bool, Integer) #> true
issubtype(Float64, Integer) #> false

typeof(5) #> Int64
isa(5, Number) #> true

max_stars_repo_path: Module 1/Chapter06/type_hierarchy.jl

max_issues_repo_path: Module 1/Chapter06/type_hierarchy.jl

max_forks_repo_path: Module 1/Chapter06/type_hierarchy.jl


text:
function makevis(x, name)
    m = Automa.compile(x)
    write("$name.dot", Automa.machine2dot(m))
    run(`dot -Tpng -o $name.png $name.dot`)
end

function appendfrom!(dst, dpos, src, spos, n)
    if length(dst) < dpos + n - 1
        resize!(dst, dpos + n - 1)
    end
    copyto!(dst, dpos, src, spos, n)
    return dst
end

newline = let
    lf = re"\n"
    lf.actions[:enter] = [:countline]

    re.cat(re.opt('\r'), lf)
end

segment = let
    hspace = re"[ \t\v]"

    seqid = re"[!-~]+"
    seqid.actions[:enter] = [:pos]
    seqid.actions[:exit] = [:seg_id]

    letters = re"[A-Za-z=.]+"
    letters.actions[:enter] = [:filledseq]

    emptyseq = re"\*"
    emptyseq.actions[:enter] = [:emptyseqs]

    sequence = letters | emptyseq
    sequence.actions[:enter] = [:pos]
    sequence.actions[:exit] = [:seg_seq]

    tag = let
        key = re"[A-Za-z][A-Za-z0-9_]"
        val = re.alt(
            re"A:[!-~]",
            re"i:[-+]?[0-9]+",
            re"f:[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?",
            re"Z:[ !-~]*",
            re"H:([0-9A-F][0-9A-F])*",
            re"B:[cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+")

        re.cat(key, ':', val)
    end
    tag.actions[:enter] = [:pos]
    tag.actions[:exit] = [:seg_tags]

    re.cat('S', hspace, seqid, hspace, sequence, re.rep(re.cat(hspace, tag)))
end

segment.actions[:enter] = [:mark]
segment.actions[:exit] = [:segment]

file = re.rep(re.cat(segment, newline))

segment_actions = Dict(
    :mark => :(@mark),
    :pos => :(pos = @relpos(p)),
    :countline => :(linenum += 1),
    :seg_id => :(segname = pos:@relpos(p - 1)),
    :seg_seq => :(seq = pos:@relpos(p - 1)),
    :emptyseq => :(println("Emptysequence!")),
    :seg_tags => :(tags = pos:@relpos(p - 1)),
    :segment => quote
        println("Segment!")
        appendfrom!(recorddata, 1, data, @markpos, p-@markpos)
        println(recorddata[segname])
        println(recorddata[seq])
        println(recorddata[tags])
    end
)

context = Automa.CodeGenContext(generator=:goto)
Automa.Stream.generate_reader(
    :gosegments,
    machine,
    actions = segment_actions,
    context = context
) |> eval

link = let
    hspace = re"[ \t\v]"

    fromid = re"[!-~]+"
    fromid.actions[:enter] = [:pos]
    fromid.actions[:exit] = [:from_id]

    fromsgn = re"[-+]"
    fromsgn.actions[:enter] = [:pos]
    fromsgn.exit[:exit] = [:from_sgn]

    toid = re"[!-~]+"
    toid.actions[:enter] = [:pos]
    toid.actions[:exit] = [:to_id]

    tosgn = re"[-+]"
    tosgn.actions[:enter] = [:pos]
    tosgn.actions[:exit] = [:to_sgn]

end

link.actions[:enter] = [:mark]
link.actions[:exit] = [:link]


text:
# Embedders
show(io::IO, embdr::E) where {E<:AbstractEmbedder} = begin
    local _vecs, _embedder, _indim, _outdim
    if E <: WordVectorsEmbedder
        _indim = size(embdr.embeddings)[1]
        _outdim = dimensionality(embdr)
        if E<:BOEEmbedder
            _embedder = "BOE"
        elseif E<:SIFEmbedder
            _embedder = "SIF"
        elseif E<:BOREPEmbedder
            _embedder = "BOREP"
        elseif E<:CPMeanEmbedder
            _embedder = "CPMean"
        elseif E<:DisCEmbedder
            _embedder = "DisC"
        else
            _embedder = "?"
        end
        L = typeof(embdr.embeddings)
        if L <: Word2Vec.WordVectors
            _vecs = "Word2Vec"
        elseif L <: Glowe.WordVectors
            _vecs = "GloVe"
        elseif L <: ConceptnetNumberbatch.ConceptNet
            _vecs = "Conceptnet"
        elseif L <: EmbeddingsAnalysis.CompressedWordVectors
            _vecs = "Compressed"
        else
            _vecs = "?"
        end
    elseif E<: DTVEmbedder
        _vecs = "DTV($(embdr.config.vectors))"
        _indim = length(embdr.model.vocab)
        _outdim = _indim
        L = typeof(embdr.model)
        if L <: StringAnalysis.LSAModel
            _embedder = "LSA"
            _outdim = dimensionality(embdr)
        elseif L <: StringAnalysis.RPModel && embdr.config.vectors_transform==:rp
            _embedder = "RP"
            _outdim = dimensionality(embdr)
        else
            _embedder = "-"
        end
    end
    printstyled(io, "$_embedder embedder, $_indim to $_outdim, $_vecs vectors")
end

# Searcher
show(io::IO, srcher::Searcher{T,E,I}) where {T,E,I} = begin
    _status = ifelse(isenabled(srcher), "enabled", "disabled")
    _status_color = ifelse(isenabled(srcher), :light_green, :light_black)
    printstyled(io, "[$_status] ", color=_status_color, bold=true)
    printstyled(io, "Searcher $(id(srcher))/")
    printstyled(io, "$(srcher.config.id_aggregation)", bold=true)
    printstyled(io, ", ")

    # Get search index type string
    if I <: NaiveIndex
        _index_type = "Naive index"
    elseif I <: BruteTreeIndex
        _index_type = "BruteTree index"
    elseif I<: KDTreeIndex
        _index_type = "KDTree index"
    elseif I <: HNSWIndex
        _index_type = "HNSW index"
    elseif I <: IVFIndex
        _index_type = "IVFADC index"
    elseif I <: NoopIndex
        _index_type = "Noop index"
    else
        _index_type = "<Unknown index>"
    end
    printstyled(io, "$_index_type", bold=true)
    #printstyled(io, "$(description(srcher))", color=:normal)
    printstyled(io, ", $(length(srcher.index)) $T embedded documents, ")
    # Get embeddings type string
    _embsstr = ifelse(srcher.data_embedder[] === srcher.input_embedder[],
                      "one embedder", "two embedders")
    printstyled(io, "$_embsstr")
end

# SearchResult
show(io::IO, result::SearchResult) = begin
    n = length(result.query_matches)
    nm = length(result.needle_matches)
    ns = length(result.suggestions)
    printstyled(io, "Search results for $(result.id): ")
    printstyled(io, " $n hits, $nm query terms, $ns suggestions.", bold=true)
end

# InternalRequest
show(io::IO, request::T) where {T<:InternalRequest} = begin
    _field_lengths = Dict(:query => 50)
    itstr = (uppercase(string(field)) * "=" *
               chop_to_length(repr(getproperty(request, field)),
                              get(_field_lengths, field, 10))
             for field in fieldnames(T))
    print(io, "InternalRequest: ", join(itstr, " | "))
end

# SearchEnv
Base.show(io::IO, env::SearchEnv{T}) where {T} = begin
    print(io, "SearchEnv{$T} with:\n")
    printstyled(io, "`-dbdata = ")
    buf = IOBuffer();
    print(buf, env.dbdata);
    seekstart(buf);
    dbstr = readuntil(buf, ':')
    printstyled(io, "$(dbstr)\n", bold=true)
    printstyled(io, "  id_key = ")
    printstyled(io, "$(env.id_key)\n", bold=true)
    printstyled(io, "  sampler = ")
    printstyled(io, "$(repr(env.sampler))\n", bold=true)
    printstyled(io, "  embedders = [\n")
    printstyled(io, "$(join(map(o->"    "*repr(o), env.embedders), "\n"))\n", bold=true)
    printstyled(io, "  ]\n")
    printstyled(io, "  searchers = [\n")
    printstyled(io, "$(join(map(o->"    "*repr(o), env.searchers), "\n"))\n", bold=true)
    printstyled(io, "  ]\n")
    printstyled(io, "  config_path = ")
    printstyled(io, "$(env.config_path)", bold=true)
end


text:
export PowerDivergenceTest, ChisqTest, MultinomialLRT

typealias Levels{T} @compat(Tuple{UnitRange{T},UnitRange{T}})

function boundproportion{T<:Real}(x::T)
  max(min(convert(Float64,x),1.0),0.0)
end

immutable PowerDivergenceTest <: HypothesisTest
  lambda::Float64
  theta0::Vector{Float64}
  stat::Float64
  df::Int64
  observed::Matrix{Int64}
  n::Int64
  thetahat::Vector{Float64}

  expected::Matrix{Float64}
  residuals::Matrix{Float64}
  stdresiduals::Matrix{Float64}
end

#test name
function testname(x::PowerDivergenceTest)
  if x.lambda == 1
    return("Pearson's Chi-square Test")
  elseif x.lambda == 0
    return("Multinomial Likelihood Ratio Test")
  elseif x.lambda == -1
    return("Minimum Discrimination Information Test")
  elseif x.lambda == -2
    return("Neyman Modified Chi-square Test")
  elseif x.lambda == -.5
    return("Freeman-Tukey Test")
  else
    return("Power Divergence Test")
  end
end

# parameter of interest: name, value under h0, point estimate
population_param_of_interest(x::PowerDivergenceTest) = ("Multinomial Probabilities", x.theta0, x.thetahat)

pvalue(x::PowerDivergenceTest; tail=:right) = pvalue(Chisq(x.df),x.stat; tail=tail)

function ci(x::PowerDivergenceTest, alpha::Float64=0.05; tail::Symbol=:both, method::Symbol=:sison_glaz, correct::Bool=true, bootstrap_iters::Int64=10000, GC::Bool=true)
  check_alpha(alpha)

  m  = length(x.thetahat)

  if tail == :left
    i = ci(x, alpha*2,method=method, GC=GC)
    @compat Tuple{Float64,Float64}[ (0.0, i[j][2]) for j in 1:m]
  elseif tail == :right
    i = ci(x, alpha*2,method=method, GC=GC)
    @compat Tuple{Float64,Float64}[ (i[j][1], 1.0) for j in 1:m]
  elseif tail == :both
    if method == :gold
      ci_gold(x,alpha,correct=correct,GC=GC)
    elseif method == :sison_glaz
      ci_sison_glaz(x,alpha, skew_correct=correct)
    elseif method == :quesenberry_hurst
      ci_quesenberry_hurst(x,alpha,GC=GC)
    elseif method == :bootstrap
      ci_bootstrap(x,alpha,bootstrap_iters)
    else
      throw(ArgumentError("method=$(method) is invalid or not implemented yet"))
    end
  else
    throw(ArgumentError("tail=$(tail) is invalid"))
  end
end

# Bootstrap
function ci_bootstrap(x::PowerDivergenceTest,alpha::Float64, iters::Int64)
  m = mapslices(x -> quantile(x,[alpha/2,1-alpha/2]), rand(Multinomial(x.n,convert(Vector{Float64},x.thetahat)),iters)/x.n, 2)
  @compat Tuple{Float64,Float64}[ ( boundproportion(m[i,1]), boundproportion(m[i,2])) for i in 1:length(x.thetahat)]
end

# Quesenberry, Hurst (1964)
function ci_quesenberry_hurst(x::PowerDivergenceTest,alpha::Float64; GC::Bool=true)
  m  = length(x.thetahat)
  cv = GC ? quantile(Chisq(1),1-alpha/m) : quantile(Chisq(m-1), 1-alpha)

  u = (cv + 2*x.observed + sqrt(  cv * (cv + 4 * x.observed .* (x.n - x.observed)/x.n)))/( 2*(x.n + cv))
  l = (cv + 2*x.observed - sqrt(  cv * (cv + 4 * x.observed .* (x.n - x.observed)/x.n)))/( 2*(x.n + cv))
  @compat Tuple{Float64,Float64}[ (boundproportion(l[j]), boundproportion(u[j])) for j in 1:m]
end

# asymptotic simultaneous confidence interval
# Gold (1963)
function ci_gold(x::PowerDivergenceTest, alpha::Float64; correct::Bool=true, GC::Bool=true)
  m  = length(x.thetahat)
  cv = GC ? quantile(Chisq(1), 1-alpha/(2*m)) : quantile(Chisq(m-1), 1-alpha/2)

  u = [ x.thetahat[j] + sqrt(cv * x.thetahat[j]*(1-x.thetahat[j])/x.n) + (correct ? 1/(2*x.n) : 0 ) for j in 1:m]
  l = [ x.thetahat[j] - sqrt(cv * x.thetahat[j]*(1-x.thetahat[j])/x.n) - (correct ? 1/(2*x.n) : 0 ) for j in 1:m]
  @compat Tuple{Float64,Float64}[ (boundproportion(l[j]), boundproportion(u[j])) for j in 1:m]
end

# Simultaneous confidence interval
# Sison, Glaz (1995)
# adapted from SAS macro by May, Johnson (2000)
function ci_sison_glaz(x::PowerDivergenceTest, alpha::Float64; skew_correct::Bool=true)
  k = length(x.thetahat)
  probn = 1/pdf( Poisson(x.n),x.n)

  c = 0
  p_old = 0

  for c in 1:x.n
    #run truncpoi
    m = zeros(k,5)
    for i in 1:k
      lambda = x.observed[i]
      #run moments
      a = lambda + c
      b = max(lambda - c, 0)
      poislama = lambda > 0 ? cdf(Poisson(lambda),a) : 1
      poislamb = lambda > 0 ? cdf(Poisson(lambda),b-1) : 1
      den = b > 0 ? poislama-poislamb : poislama

      mu  = zeros(4,1)
      for r in 1:4
        poisA = 0
        poisB = 0
        plar = lambda > 0 ? cdf(Poisson(lambda),a-r) : 1
        plbr = lambda > 0 ? cdf(Poisson(lambda),b-r-1): 1
        poisA = a - r >= 0 ? poislama-plar : poislama
        if  b - r - 1 >= 0
          poisB = poislamb - plbr
        end
        if b - r - 1 < 0 && b-1 >= 0
          poisB = poislamb
        end
        if b - r - 1 < 0 && b-1 < 0
          poisB = 0
        end
        mu[r] = (lambda^r) * (1-(poisA-poisB)/den)
      end
      # end of moments
      m[i,1] = mu[1]
      m[i,2] = mu[2] + mu[1] - mu[1]^2
      m[i,3] = mu[3]+mu[2]*(3-3*mu[1])+(mu[1]-3*mu[1]^2+2*mu[1]^3)
      m[i,4] = mu[4] + mu[3] * (6-4*mu[1]) + mu[2] * (7-12*mu[1]+6*mu[1]^2) + mu[1]-4*mu[1]^2+6*mu[1]^3-3*mu[1]^4
      m[i,5] = den
    end
    [ m[i,4] = m[i,4] - 3*m[i,2]^2 for i in 1:k]
    s1,s2,s3,s4 = mapslices(sum,m,1)
    z  = (x.n-s1)/sqrt(s2)
    g1 = s3/(s2^(3/2))
    g2 = s4*(s2^2)

    poly = 1 + g1 * (z^3 - 3*z)/6 + g2*(z^4-6*z^2+3)/24 + g1^2 * (z^6 - 15*z^4 + 45*z^2 -15)/72
    f = poly*exp(-z^2/2)/(sqrt(2*pi))
    probx = 1

    [ probx *= probx*m[i,5] for i in 1:k]

    p = probn*probx*f/sqrt(s2)
    # end of truncpoi

    p > 1-alpha && p_old < 1-alpha ? break : nothing
    p_old = p
  end

  delta = (1-alpha - p_old)/(p_old)
  out = zeros(k,5)
  num = zeros(k,1)

  c = c-1

  vol1 = 1
  vol2 = 1
  for i in 1:k
    num[i,1] = i
    cn = c/x.n
    onen = 1/x.n
    out[i,1] = x.thetahat[i]

    out[i,2] = max(x.thetahat[i]-cn,0)
    out[i,3] = min(x.thetahat[i]+cn+2*delta/x.n,1)

    out[i,4] = max(x.thetahat[i] - cn - onen,0)
    out[i,5] = min(x.thetahat[i] + cn + onen,1)
  end
  if skew_correct
    return( @compat Tuple{Float64,Float64}[ (boundproportion(out[i,2]),boundproportion(out[i,3])) for i in 1:k] )
  else
    return( @compat Tuple{Float64,Float64}[ (boundproportion(out[i,4]),boundproportion(out[i,5])) for i in 1:k] )
  end
end

# power divergence statistic for testing goodness of fit
# Cressie and Read 1984; Read and Cressie 1988
# lambda  =  1: Pearson's Chi-square statistic
# lambda ->  0: Converges to Likelihood Ratio test stat
# lambda -> -1: Converges to minimum discrimination information statistic (Gokhale, Kullback 1978)
# lambda =  -2: Neyman Modified chi-squared statistic (Neyman 1949)
# lambda = -.5: Freeman-Tukey statistic (Freeman, Tukey 1950)

# Under regularity conditions, their asymptotic distributions are all the same (Drost 1989)
# Chi-squared null approximation works best for lambda near 2/3
function PowerDivergenceTest{T<:Integer,U<:AbstractFloat}(x::AbstractMatrix{T}; lambda::U=1.0, theta0::Vector{U} = ones(length(x))/length(x))

  nrows, ncols = size(x)
  n = sum(x)

  #validate date
  any( x .< 0) || any( !isfinite(x)) ? error("all entries must be nonnegative and finite") : nothing
  n == 0 ? error("at least one entry must be positive") : nothing
  isfinite(nrows) && isfinite(ncols) && isfinite(nrows*ncols) ? nothing : error("Invalid number of rows or columns")

  df = thetahat = xhat = V = 0

  if nrows > 1 && ncols > 1
    rowsums = mapslices(sum, x, 2)
    colsums = mapslices(sum, x, 1)
    df = (nrows-1)*(ncols-1)
    thetahat = [ x[i]/n for i in 1:length(x)]
    xhat = rowsums * colsums / n
    theta0 = (xhat/n)[:]
    V = Float64[ colsums[j]*rowsums[i]*(n-rowsums[i])*(n-colsums[j])/n^3 for i in 1:nrows, j in 1:ncols]
  elseif nrows == 1 || ncols == 1
    df = length(x) - 1
    xhat = reshape(n * theta0, size(x))
    thetahat = x/n
    V = reshape(n * theta0 .* (1-theta0), size(x))

    abs( 1 - sum(theta0)) <= sqrt(eps()) ? nothing : error("Probabilities must sum to one")
  else
    error("Number of rows or columns must be at least 1")
  end

  stat = 0
  if lambda == 0
    for i in 1:length(x)
      stat += x[i]*(log(x[i]) - log(xhat[i]))
    end
    stat *= 2
  elseif lambda == -1
    for i in 1:length(x)
      stat += xhat[i]*(log(xhat[i]) - log(x[i]))
    end
    stat *= 2
  else
    for i in 1:length(x)
      stat += x[i]*( (x[i]/(xhat[i]))^lambda - 1)
    end
    stat *= 2/(lambda*(lambda+1))
  end

  PowerDivergenceTest(lambda, theta0, stat, df, x, n, thetahat[:], xhat, (x - xhat)./sqrt(xhat), (x-xhat)./sqrt(V))
end

#convenience functions

#PDT
function PowerDivergenceTest{T<:Integer,U<:AbstractFloat}(x::AbstractVector{T}, y::AbstractVector{T}, levels::Levels{T}; lambda::U=1.0)
  d = counts(x,y,levels)
  PowerDivergenceTest(d,lambda=lambda)
end

function PowerDivergenceTest{T<:Integer,U<:AbstractFloat}(x::AbstractVector{T}, y::AbstractVector{T}, k::T; lambda::U=1.0)
  d = counts(x,y,k)
  PowerDivergenceTest(d,lambda=lambda)
end

PowerDivergenceTest{T<:Integer,U<:AbstractFloat}(x::AbstractVector{T}; lambda::U=1.0, theta0::Vector{U} = ones(length(x))/length(x)) =
  PowerDivergenceTest(reshape(x,length(x),1),lambda=lambda,theta0=theta0)

#ChisqTest
function ChisqTest{T<:Integer}(x::AbstractMatrix{T})
  PowerDivergenceTest(x,lambda=1.0)
end

function ChisqTest{T<:Integer}(x::AbstractVector{T}, y::AbstractVector{T}, levels::Levels{T})
  d = counts(x,y,levels)
  PowerDivergenceTest(d,lambda=1.0)
end

function ChisqTest{T<:Integer}(x::AbstractVector{T}, y::AbstractVector{T}, k::T)
  d = counts(x,y,k)
  PowerDivergenceTest(d,lambda=1.0)
end

ChisqTest{T<:Integer,U<:AbstractFloat}(x::AbstractVector{T}, theta0::Vector{U} = ones(length(x))/length(x)) =
  PowerDivergenceTest(reshape(x,length(x),1), lambda=1.0, theta0 = theta0)

#MultinomialLRT
function MultinomialLRT{T<:Integer}(x::AbstractMatrix{T})
  PowerDivergenceTest(x,lambda=0.0)
end

function MultinomialLRT{T<:Integer}(x::AbstractVector{T}, y::AbstractVector{T}, levels::Levels{T})
  d = counts(x,y,levels)
  PowerDivergenceTest(d,lambda=0.0)
end

function MultinomialLRT{T<:Integer}(x::AbstractVector{T}, y::AbstractVector{T}, k::T)
  d = counts(x,y,k)
  PowerDivergenceTest(d,lambda=0.0)
end

MultinomialLRT{T<:Integer,U<:AbstractFloat}(x::AbstractVector{T}, theta0::Vector{U} = ones(length(x))/length(x)) =
  PowerDivergenceTest(reshape(x,length(x),1), lambda=0.0, theta0 = theta0)

function show_params(io::IO, x::PowerDivergenceTest, ident="")
  println(io,ident, "Sample size:        $(x.n)")
  println(io,ident, "statistic:          $(x.stat)")
  println(io,ident, "degrees of freedom: $(x.df)")
  println(io,ident, "residuals:          $(x.residuals[:])")
  println(io,ident, "std. residuals:     $(x.stdresiduals[:])")
end


text:
module titanic

using CSV
using DataFrames
using JuMP
using LinearAlgebra
using OSQP
using Parameters
using Statistics

include("dataset.jl")
include("io.jl")
include("logistic.jl")
include("adaboost.jl")
include("svm.jl")

export loadtrain
export loadtest
export savepredictions

export regression
export boost
export svm
export classify
export classificationerror

end


text:
using ApproxFunRational, ApproxFunFourier, ApproxFunBase, ApproxFun
using Test

## TODO: Test SumFuns

@testset "ApproxFunRational.jl: Fourier transform" begin
    f = x -> 1. + (x+x^2)*exp(-x^2+1im*x^2)
    α = 2.0
    L = 2.0
    F = Fun(cai(f,1.0),OscLaurent(α,L))
    𝓕 = FourierTransform(3.0)
    𝓗 = FourierTransform(-3.0)
    FF = 𝓗*(𝓕*F)
    @test F(.1) ≈ FF(.1)

    f = x -> (x+x^2)*exp(-x^2+1im*x^2)
    α = 2.0
    L = 2.0
    F = Fun(zai(f),OscLaurent(α,L))
    𝓕 = FourierTransform(1.0)
    Fhat = 𝓕*F
    @test Fhat(.1) ≈ -0.1979914917932347 + 0.2718859895675688im

    f = x -> (x+x^2)*exp(-x^2+1im*x^2)
    α = 2.0
    L = 2.0
    F = Fun(zai(f),OscLaurent(α,L))
    𝓕 = FourierTransform(-1.0)
    Fhat = 𝓕*F
    @test Fhat(.1) ≈ -0.20287825133470003 + 0.22704533835040858im

    f = x -> 1. + (x+x^2)*exp(-x^2+1im*x^2)
    α = -2.0
    L = 2.0
    F = Fun(cai(f,1.0),OscLaurent(α,L),4)
    𝓕 = FourierTransform(-1.0)
    𝓗 = FourierTransform(1.0)
    FF = 𝓗*(𝓕*F)
    @test F(.1)-FF(.1) ≈ 0
end

@testset "ApproxFunRational.jl: Matrix-vector function product" begin
    L = 1.; α = -2.; β = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    f = z -> exp(-z^2+z)
    h = z -> 1.0 + 1/(z^2+1)
    F = z -> [cai(h,1.0)(z) zai(g)(z); zai(f)(z) cai(h,1.0)(z)]
    G = z -> [zai(g)(z), zai(f)(z)]
    FF = Fun(F,OscLaurent(0.0))
    GG = Fun(G,OscLaurent(0.0))
    J = z -> [zai(g)(z)^2 + zai(f)(z)^2]
    JJ = transpose(GG)*GG
    HH = FF*GG
    @test JJ(.1) ≈ J(.1)
    @test HH(.1) ≈ F(.1)*G(.1)
    FF = Fun(F,OscLaurent(β))
    GG = Fun(G,OscLaurent(0.0)) + Fun(G,OscLaurent(α))
    @test (FF*GG)(.1) ≈ F(.1)*G(.1)*exp(1im*β*(.1)) + F(.1)*G(.1)*exp(1im*(α+β)*(.1))
end

@testset "ApproxFunRational.jl: Vector-valued functions" begin
    L = 1.; α = -2.; β = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    f = z -> exp(-z^2+z)
    h = z -> 1.0 + 1/(z^2+1)
    F = z-> [zai(g)(z), zai(f)(z), cai(h,1.0)(z)]
    FF = Fun(F,OscLaurent(α))
    @test F(.1)*exp(1im*α*(.1)) ≈ FF(.1)
end

@testset "ApproxFunRational.jl: Oscillatory Cauchy integrals" begin
    L = 1.; α = -2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    G1 = Fun(zai(g), OscLaurent(dom,α), 200)
    GCP = Cauchy(1)*G1
    GCM = Cauchy(-1)*G1
    @test GCP(.1) - GCM(.1) ≈ G1(.1)
    a1 = 0.027382689548799077 + 0.0017439616653601026im
    @test GCP(.1) ≈ a1
    a2 = -0.9478038907592183 + 0.1994240679809568im
    @test GCM(.1) ≈ a2

    L = 1.; α = -2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> exp(-z^2+z)
    G2 = Fun(zai(g), OscLaurent(dom,α), 1000)
    GCP = Cauchy(1)*G2
    GCM = Cauchy(-1)*G2
    @test GCP(.1) - GCM(.1) ≈ G2(.1)
    b1 = 0.030786661098655133 - 0.09357539334812659im
    @test GCP(.1) ≈ b1
    b2 = -1.0415769846942182 + 0.12380347936876926im
    @test GCM(.1) ≈ b2

    L = 1.; α = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> exp(-z^2+z)
    G3 = Fun(zai(g), OscLaurent(dom,α), 1000)
    GCP = Cauchy(1)*G3
    GCM = Cauchy(-1)*G3
    @test GCP(.1) - GCM(.1) ≈ G3(.1)
    c1 = 1.0415769846942182 + 0.12380347936876926im
    @test GCP(.1) ≈ c1
    c2 = -0.030786661098655133 - 0.09357539334812659im
    @test GCM(.1) ≈ c2

    GCP = ApproxFunRational.CauchyP(G1 + G3)
    @test GCP(.1) ≈ a1 + c1

    L = 2.; α = -2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    G = Fun(zai(g), OscLaurent(dom,α), 300)
    GCP = Cauchy(1)*G
    GCM = Cauchy(-1)*G
    @test GCP(.1) - GCM(.1) ≈ G(.1)
    @test GCP(.1) ≈ 0.027382689548799077 + 0.0017439616653601026im
    @test GCM(.1) ≈ -0.9478038907592183 + 0.1994240679809568im

    L = 0.5; α = -2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    G = Fun(zai(g), OscLaurent(dom,α), 400)
    GCP = Cauchy(1)*G
    GCM = Cauchy(-1)*G
    @test GCP(.1) - GCM(.1) ≈ G(.1)
    @test GCP(.1) ≈ 0.027382689548799077 + 0.0017439616653601026im
    @test GCM(.1) ≈ -0.9478038907592183 + 0.1994240679809568im
end

@testset "ApproxFunRational.jl: Evaluation and multiplication" begin
    L = 1.; α = -2.; β = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    G = Fun(zai(g), OscLaurent(dom,α))
    println(length(G.coefficients))
    @test G(.1) ≈ g(.1)*exp(α*im*.1) #test adaptivity

    f = z -> exp(-z^2+z)
    G = Fun(zai(g), OscLaurent(dom,α), 200)
    @test G(.1) ≈ g(.1)*exp(α*im*.1)
    F = Fun(zai(f), OscLaurent(dom,β), 200)
    @test F(.1) ≈ f(.1)*exp(β*im*.1)
    H = F*G
    @test H(.1) ≈ F(.1)*G(.1)

    L = 2.; α = -2.; β = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    f = z -> exp(-z^2+z)
    G = Fun(zai(g), OscLaurent(dom,α), 500)
    @test G(.1) ≈ g(.1)*exp(α*im*.1)
    F = Fun(zai(f), OscLaurent(dom,β), 200)
    @test F(.1) ≈ f(.1)*exp(β*im*.1)
    H = F*G
    @test H(.1) ≈ F(.1)*G(.1)

    L = 0.5; α = -2.; β = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    f = z -> exp(-z^2+z)
    G = Fun(zai(g), OscLaurent(dom,α), 400)
    @test G(.1) ≈ g(.1)*exp(α*im*.1)
    F = Fun(zai(f), OscLaurent(dom,β), 400)
    @test F(.1) ≈ f(.1)*exp(β*im*.1)
    H = F*G
    @test H(.1) ≈ F(.1)*G(.1)

    f = z -> exp(-z^2+z*1im)*(1.0+1im)
    F = Fun(zai(f), OscLaurent(dom,β), 200)
    Fbar = conj(F)
    @test Fbar(.1) ≈ conj(F(.1))
end

@testset "ApproxFunRational.jl: Differentiation and integration" begin
    L = 1.; α = -2.; β = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    f = z -> exp(-z^2+z)
    df = z-> (1-2z)*exp(-z^2+z)
    G = Fun(zai(g), OscLaurent(dom,α), 200)
    F = Fun(zai(f), OscLaurent(dom,β), 200)
    H = F*G
    dF = F'
    dF_2 = Derivative()*F
    @test dF(.1) ≈ df(.1)*exp(β*im*.1) + (β*im)*f(.1)*exp(β*im*.1)
    @test dF_2(.1) ≈ df(.1)*exp(β*im*.1) + (β*im)*f(.1)*exp(β*im*.1)
    @test sum(H) ≈ sqrt(pi)
    @test sum(G) ≈ pi*sech(1.0*pi)
    @test sum(F) ≈ exp(-3/4 + im)*sqrt(pi)

    L = 2.; α = -2.; β = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    f = z -> exp(-z^2+z)
    df = z-> (1-2z)*exp(-z^2+z)
    G = Fun(zai(g), OscLaurent(dom,α), 400)
    F = Fun(zai(f), OscLaurent(dom,β), 200)
    H = F*G
    dF = F'
    dF_2 = Derivative()*F
    @test dF(.1) ≈ df(.1)*exp(β*im*.1) + (β*im)*f(.1)*exp(β*im*.1)
    @test dF_2(.1) ≈ df(.1)*exp(β*im*.1) + (β*im)*f(.1)*exp(β*im*.1)
    @test sum(H) ≈ sqrt(pi)
    @test sum(G) ≈ pi*sech(1.0*pi)
    @test sum(F) ≈ exp(-3/4 + im)*sqrt(pi)

    L = 0.5; α = -2.; β = 2.;
    dom = PeriodicLine{false,Float64}(0.,L)
    g = z -> sech(z)
    f = z -> exp(-z^2+z)
    df = z-> (1-2z)*exp(-z^2+z)
    G = Fun(zai(g), OscLaurent(dom,α), 400)
    F = Fun(zai(f), OscLaurent(dom,β), 400)
    H = F*G
    dF = F'
    dF_2 = Derivative()*F
    @test dF(.1) ≈ df(.1)*exp(β*im*.1) + (β*im)*f(.1)*exp(β*im*.1)
    @test dF_2(.1) ≈ df(.1)*exp(β*im*.1) + (β*im)*f(.1)*exp(β*im*.1)
    @test sum(H) ≈ sqrt(pi)
    @test sum(G) ≈ pi*sech(1.0*pi)
    @test sum(F) ≈ exp(-3/4 + im)*sqrt(pi)
end


text:
module TestExamplesUsefulPatterns
using LiterateTest.AssertAsTest: @assert
include("../examples/useful_patterns.jl")
end  # module


text:
module DrugInteractions

export drug_interactions_app

using Gtk
using RxNav

const _apps = GtkWindow[]
const _apps_should_persist = [true]

"""
    drug_interactions_app(title = "Drug Interaction Queries", rlabel = "Results")

Create a `Gtk` widget with entries for one or more substances to check for interactions.
"""
function drug_interactions_app(title = "Drug Interaction Queries", rlabel = "Results")
    label = GtkLabel("Drug(s) to Check:  ")
    substances = GtkEntry()
    topbox = GtkBox(:h)
    push!(topbox, label, substances)
    highonly = GtkCheckButton("Only Search For High Severity")
    set_gtk_property!(highonly, :active, true)
    resultbutton = GtkButton(rlabel)
    win = GtkWindow(title, 300, 100) |> (GtkFrame() |> (vbox = GtkButtonBox(:v)))
    push!(_apps, win)
    push!(vbox, topbox, highonly, resultbutton)
    set_gtk_property!(substances, :expand, true)

    function queryRxNav(w)
        high = get_gtk_property(highonly, :active, Bool)
        drugs = split(something(get_gtk_property(substances, :text, String), ""), r"\s+")
        if !isempty(drugs)
            tuples = length(drugs) == 1 ? interaction(first(drugs); ONCHigh = high) :
                                          interaction_within_list(string.(drugs))
            if !isempty(tuples)
                lines = "Substance 1     Substance 2    Severity    Description\n" * "-"^160 * "\n"
                for t in tuples
                    line = rpad(t[1], 20) * " " * rpad(t[2], 20) * " " * rpad(t[3], 8) * " " * t[4]
                    lines *= line * "\n"
                end
                info_dialog(lines)
            else
                info_dialog("No results found")
            end
        end
    end

    signal_connect(queryRxNav, resultbutton, :clicked)

    if _apps_should_persist[1]
        cond = Condition()
        endit(w) = notify(cond)
        signal_connect(endit, win, :destroy)
        showall(win)
        wait(cond)
    else
        sleep(5)
    end
end

end # module


text:
module VoronoiGraph

using LinearAlgebra
using NearestNeighbors
using Polyhedra
using ProgressMeter
using SparseArrays
using StaticArrays
using RecipesBase

include("voronoi.jl")
include("raycast.jl")
include("volume.jl")
include("plot.jl")
include("montecarlo.jl")

export voronoi, voronoi_random
export volumes
export mc_volumes, mc_integrate

end


text:
@testset "Map" begin
  @testset "Existence of (pre)image" begin
    G = DiagonalGroup([3, 5, 0])
    H, mH = sub(G, [G[1]])
    a = G[1]
    b = G[3]
    bb, c = @inferred haspreimage(mH, a)
    @test bb
    @test mH(c) == a
    bb, c = @inferred haspreimage(mH, b)
    @test !bb
    # TODO: Test for hasimage missing
  end

  @testset "Homomorphisms" begin
    G = DiagonalGroup([4, 4, 4])
    H = DiagonalGroup([4, 4, 4])
    h = @inferred hom(gens(G), gens(H))
    @test h(G[1]) == H[1]
    @test h(G[2]) == H[2]
    @test h(G[3]) == H[3]

    h = @inferred hom(G, [2*h for h in gens(H)])
    @test h(G[1]) == 2*H[1]
    @test h(G[2]) == 2*H[2]
    @test h(G[3]) == 2*H[3]
  end

  @testset "Kernel" begin
    G = DiagonalGroup([4, 4, 4])
    H = DiagonalGroup([4, 4, 4])
    h = @inferred hom(G, [2*h for h in gens(H)])
    @test h(G[1]) == 2*H[1]
    @test h(G[2]) == 2*H[2]
    @test h(G[3]) == 2*H[3]

    K, mK = @inferred kernel(h)
    @test all(iszero(h(mK(k))) for k in K)
    @test order(K) == 8
  end

  @testset "Image" begin
    G = DiagonalGroup([4, 4, 4])
    H = DiagonalGroup([4, 4, 4])
    h = @inferred hom(G, [2*h for h in gens(H)])
    @test h(G[1]) == 2*H[1]
    @test h(G[2]) == 2*H[2]
    @test h(G[3]) == 2*H[3]

    I, mI = @inferred image(h)
    @test all(haspreimage(h, mI(i))[1] for i in I)
    @test order(I) == 8
  end

  @testset "Injectivity" begin
    G = DiagonalGroup([4, 4, 4])
    H = DiagonalGroup([4, 4, 4])

    h = @inferred hom(G, [2*h for h in gens(H)])
    b = @inferred isinjective(h)
    @test !b

    h = @inferred hom(G, [3*h for h in gens(H)])
    b = @inferred isinjective(h)
    @test b
  end

  @testset "Surjectivity" begin
    G = DiagonalGroup([4, 4, 4])
    H = DiagonalGroup([4, 4, 4])

    h = @inferred hom(G, [2*h for h in gens(H)])
    b = @inferred issurjective(h)
    @test !b

    h = @inferred hom(G, [3*h for h in gens(H)])
    b = @inferred issurjective(h)
    @test b
  end

  @testset "Bijectivity" begin
    G = DiagonalGroup([4, 4, 4])
    H = DiagonalGroup([4, 4, 4])

    h = @inferred hom(G, [2*h for h in gens(H)])
    b = @inferred isbijective(h)
    @test !b

    h = @inferred hom(G, [3*h for h in gens(H)])
    b = @inferred isbijective(h)
    @test b
  end
end


text:
@testset "Testing Main deconvolution code" begin
    Random.seed!(42)

    img = [0.5560268761463861 0.29948409035891055 0.46860588216767457; 0.444383357109696 1.7778610980573246 0.15614346264074028; 0.027155338009193845 1.14490153172882 2.641991008076796]
    psf = [1.0033099014594844 0.5181487878771377 0.8862052960481365; 1.0823812056084292 1.4913791170403063 0.6845647041648603; 0.18702790710363 0.3675627461748204 1.590579974922555]

    res = [2.4393835034275493 0.013696697097634842 0.0002833052222499294; 0.07541628019133978 1.0066536888249171 0.02222160874466724; 0.0004945773667781262 0.008547708184955495 3.717245734531717]

    #@show deconvolution(img, psf, λ=0.01)[1]
    @test all(≈(res, deconvolution(img, psf, λ=0.01)[1], rtol=0.1))
    @test all(≈(res, deconvolution(img, psf, λ=0.01)[1], rtol=0.1))
    @test all(≈(res, deconvolution(img, psf, λ=0.01, iterations=20)[1], rtol=0.1))

    # testing regularizer
    res2 = [4.188270526990536 5.999388400251461e-10 2.8299849680327642e-8; 1.725273124171714e-7 2.54195544512864 2.0216187854619135e-9; 9.594324085846738e-10 1.2000166997002865e-8 0.7863126081711094]
    @test all(≈(res2, deconvolution(img, psf, regularizer=nothing, padding=0.0)[1], rtol=0.1))

    # testing padding
    img = [8.21306764808666 10.041589152470781 86.74936458947307; 17.126996611046078 4.324960146596254 11.39657297820361; 21.019754207225656 14.128485444028698 28.441178191470662]

    psf = [0.37240087577993225 0.562668812321259 0.6810849274435286; 0.36901028455183293 0.10686911035365092 1.3391251213773154; 0.007612980079313577 0.5694584949295476 0.23828371819888622]

    #@show deconvolution(img, psf, regularizer=nothing, padding=0.1)[1]
    res3 = [3.700346192004195 16.163004402781457 0.0005317269571170434; 0.3852386808335988 14.378882057906575 0.04691706650167405; 1.5059695704739982 16.94303953714345 22.72731111751148]
    @test all(≈(res3, deconvolution(img, psf, regularizer=nothing, padding=0.1)[1], rtol=1e-2))

    # test without mapping
    res4 = [-13.085096066984729 34.39174935297922 -11.353417020874208; -13.01079706662783 43.76596781851713 -7.565144283495296; 16.740081837985805 -7.488374542587605 37.666978022259336]
    #= @show deconvolution(img, psf, regularizer=nothing, padding=0.1, mapping=nothing)[1] =#
    @test all(≈(res4, deconvolution(img, psf, regularizer=nothing, padding=0.1, mapping=nothing)[1], rtol=0.1))

    # test OptimPackNextGen optimizers  (which is currently not officially released yet)
    res5 = [0.49633      16.4936   0.00862045; 0.0139499    15.8587   2.31716; 0.000227487   9.263   19.5289]
    @test all(≈(res5, deconvolution(img, psf, opt=vmlmb!, opt_options=(mem=20, lower=0, lnsrch=OptimPackNextGen.LineSearches.MoreThuenteLineSearch()),
    regularizer=nothing, padding=0.1, mapping=nothing, opt_package=Opt_OptimPackNextGen)[1], rtol=0.1))

    # test broadcasting with image having more dimensions
    img = zeros((3, 3, 2))
    imgc = abs.(randn((3, 3, 1)))
    img[:, :, 1] = imgc
    img[:, :, 2] = imgc
    psf = zeros((3, 3))
    psf[1,1] = 1

    res = deconvolution(img, psf, regularizer=GR(num_dims=3, sum_dims=[1,2]))[1]
    @test all(res[:, :, 1] .≈ res[:, :, 2])

    img = zeros((3, 3, 2, 1))
    imgc = abs.(randn((3, 3, 1, 1)))
    img[:, :, 1, 1] = imgc
    img[:, :, 2, 1] = imgc
    psf = zeros((3, 3))
    psf[1,1] = 1

    res = deconvolution(img, psf, regularizer=GR(num_dims=4, sum_dims=[1,2]))[1]
    @test all(≈(res[:, :, 1, :], res[:, :, 2, :], rtol=0.1))

end

@testset "Compare optimization with iterative lucy richardson scheme" begin

    img = Float32.(testimage("resolution_test_512"));
    psf = Float32.(generate_psf(size(img), 30));
    img_b = conv(img, psf);
    img_n = poisson(img_b, 300);

    reg = GR()
    res = richardson_lucy_iterative(img_n, psf, regularizer=reg, iterations=200);
    res2, o = deconvolution(img_n, psf, regularizer=reg, iterations=30);
    @test ≈(res .+ 1, res2 .+ 1, rtol=0.003)

    reg = TV()
    res = richardson_lucy_iterative(img_n, psf, iterations=500, λ=0.005, regularizer=reg);
    res2, o = deconvolution(img_n, psf, iterations=35, regularizer=reg, λ=0.005);
    @test ≈(res2 .+1, res .+ 1, rtol=0.02)

    reg = nothing
    res = richardson_lucy_iterative(img_n, psf, regularizer=reg, iterations=400);
    res2, o = deconvolution(img_n, psf, regularizer=reg, iterations=40);
    @test ≈(res2 .+1, res .+ 1, rtol=0.02)

end