text:
struct Params{FT}
    θ0 :: FT
    π0 :: FT
    Π :: FT
    S :: FT
    f :: FT
end

function Params(T = Float64;
    θ0 = 285, θT = 340, π0 = 1000, πT = 500, f = 4π/86400*sind(45)
)
    Π = π0 - πT
    S = (θT - θ0)/Π
    θ0 = θ0/(S*Π)
    π0 = π0/Π
    return Params(T(θ0), T(π0), T(Π), T(S), T(f))
end

get_Rd(p::Params) = sqrt(p.S*p.Π^2/p.f^2)
float_type(p::Params{T}) where {T} = T

struct Domain{R,I,FT,A}
    x :: R
    y :: R
    z :: R
    nx :: I
    ny :: I
    nz :: I
    hx :: I
    hy :: I
    hz :: I
    Δx :: FT
    Δy :: FT
    Δz :: FT
    xc :: A
    yc :: A
    zc :: A
end

function Domain(p::Params{T};
    x = (-3, 3), y = (-3, 3),
    size = (32, 32, 32), halo = (1,1,1)
) where {T}

    π0 = p.π0
    x = T.(x)
    y = T.(y)
    z = (π0, π0 - 1)
    nx, ny, nz = size
    hx, hy, hz = halo
    Δx = (x[2] - x[1])/nx
    Δy = (y[2] - y[1])/ny
    Δz = (z[2] - z[1])/nz
    xc = range(x[1] + (0.5 - hx)*Δx, x[2] - (0.5 - hx)*Δx, length = nx + 2*hx)
    xc = OffsetArray(xc, -hx)
    yc = range(y[1] + (0.5 - hy)*Δy, y[2] - (0.5 - hy)*Δy, length = ny + 2*hy)
    yc = OffsetArray(yc, -hy)
    zc = range(z[1] + (0.5 - hz)*Δz, z[2] - (0.5 - hz)*Δz, length = nz + 2*hz)
    zc = OffsetArray(zc, -hz)
    return Domain(x, y, z, nx, ny, nz, hx, hy, hz, T(Δx), T(Δy), T(Δz), xc, yc, zc)

end

halo(d::Domain) = (d.hx, d.hy, d.hz)
Base.size(d::Domain) = (d.nx, d.ny, d.nz)
Base.IteratorsMD.CartesianIndices(d::Domain) = CartesianIndices((d.nx, d.ny, d.nz))
Base.LinearIndices(d::Domain) = LinearIndices((d.nx, d.ny, d.nz))
Base.getindex(d::Domain, I::CartesianIndex{3}) = (d.xc[I[1]], d.yc[I[2]], d.zc[I[3]])
Base.getindex(d::Domain, i, j, k) = (d.xc[i], d.yc[j], d.zc[k])
const dx = CartesianIndex(1,0,0)
const dy = CartesianIndex(0,1,0)
const dz = CartesianIndex(0,0,1)
function new_field(d::Domain)
    underlying_data = zeros(d.nx + 2d.hx, d.ny + 2d.hy, d.nz + 2d.hz)
    return OffsetArray(underlying_data, -d.hx, -d.hy, -d.hz)
end

function new_rhs(d::Domain)
    return zeros(prod(size(d)))
end
function new_linearized_rhs(d::Domain)
    return zeros(2*prod(size(d)))
end
function field_from_rhs!(f, rhs, d::Domain)
    l = LinearIndices(d)
    c = CartesianIndices(d)
    for I = c
        f[I] = rhs[l[I]]
    end
    return f
end
function field_from_rhs(rhs, d::Domain)
    f = new_field(d)
    return field_from_rhs!(f, rhs, d)
end
function fields_from_linearized_rhs!(ϕ, ψ, rhs, d::Domain)
    c = CartesianIndices(d)
    l = LinearIndices(d)
    ol = last(l)
    xϕ = @view rhs[1:ol]
    xψ = @view rhs[ol+1:end]
    for I = c
        ϕ[I] = xϕ[l[I]]
        ψ[I] = xψ[l[I]]
    end
    return (ϕ, ψ)
end
function fields_from_linearized_rhs(rhs, d::Domain)
    ϕ = new_field(d)
    ψ = new_field(d)
    return fields_from_linearized_rhs!(ϕ, ψ, rhs, d)
end
function rhs_from_field!(rhs, f, d::Domain)
    l = LinearIndices(d)
    c = CartesianIndices(d)
    for I = c
        rhs[l[I]] = f[I]
    end
    return rhs
end
function rhs_from_field(f, d::Domain)
    rhs = new_rhs(d)
    return rhs_from_field!(rhs, f, d)
end

mutable struct Convergence{FT,I}
    Δ :: FT
    n :: I
end

struct Solver{FT}
    atol :: FT
    convergence :: Convergence{FT}
    ω :: FT
end

function Solver(p::Params; atolϕ = 1e0, ω = 0.7)
    S = p.S
    Π = p.Π
    atol = atolϕ/(S*Π^2)
    T = float_type(p)
    return Solver(T(atol), Convergence(T(0), 0), T(ω))
end

function is_converged(s::Solver)
    return s.convergence.Δ ≤ s.atol
end

function set_Δ!(s::Solver, Δ)
    s.convergence.Δ = Δ
    return s
end

function increment!(s::Solver, n = 1)
    s.convergence.n += n
    return s
end

function print_convergence_info(s::Solver)
    n = s.convergence.n
    Δ = s.convergence.Δ
    tol = s.atol
    print(n, "\tΔ = ", Δ, "; ", Δ/tol, "x tolerance\n")
end


text:
using CompScienceMeshes

coarse = meshsphere(1.0, 0.45)
fine = CompScienceMeshes.lineofsight_refinement(m)
vtoc, nc = CompScienceMeshes.vertextocellmap(fine)

colors = zeros(numcells(fine))
for i in 1:numvertices(coarse)
    c = rand()
    for j in 1 : (nc[i])
        colors[vtoc[i,j]] = c
    end
end

import PlotlyJS
PlotlyJS.plot(patch(fine.mesh, colors))

ext: jl


text:
mutable struct Canvas
	x::Int64
	y::Int64
	clusters::Vector{Cluster}
	num_cluster::Int64
end

MAX_ITERATION = 100000

Canvas(x::Int64,y::Int64) = Canvas(x, y, [], 0)

function initialize(canvas::Canvas, num_clusters::Int64, num_waypoints::Int64, trasmission_range::Int64)
	waypoints_per_cluster = random_array_fixed_sum(num_waypoints,num_clusters)
	canvas.num_cluster = num_clusters
	# For each cluster
	for i = 1:num_clusters
		push!(canvas.clusters,define_cluster(canvas, waypoints_per_cluster[i], trasmission_range))
	end
end

function define_cluster(canvas::Canvas, waypoints_per_cluster::Int, trasmission_range)
	cluster = Cluster(trasmission_range,canvas.x,canvas.y)
	groups = random_array_fixed_sum(waypoints_per_cluster,rand(1:waypoints_per_cluster))
	last_plotted_point = Tuple{Int64,Int64}((0,0))

	for i in 1:length(groups)
		group_points = Tuple{Int64,Int64}[]
		# println("Plotting $i group")
		# First we plot the first point in the cluster, which has a different rule
		if length(cluster.points) == 0  ##If I'm plotting the very first point, plot it with an uniform distribution on the canvas
			last_plotted_point = generate_point(canvas)
			# println("First Point $i group : $last_plotted_point")
			add_point_to_cluster(cluster,last_plotted_point)
		else #The first point of the other groups are plotted with a Y/4 - Y/3 distance from the last plotted point
			last_plotted_point = generate_group_first_point(canvas,last_plotted_point,cluster)
			# println("Point $i group : $last_plotted_point")
			add_point_to_cluster(cluster,last_plotted_point)
		end
		push!(group_points,last_plotted_point)

		for j in 2:groups[i] #Plot the other ones within 0.1R from the last plotted point
			last_plotted_point = generate_group_internal_point(canvas,trasmission_range,cluster,group_points)
			push!(group_points,last_plotted_point)
			add_point_to_cluster(cluster,last_plotted_point)
		end
	end
	return cluster
end

#TODO: Collapse point generation point functions into one, with acceptance funcion as parameter
function generate_group_internal_point(canvas::Canvas,R::Int64,cluster::Cluster,group_points::Array{Tuple{Int64,Int64}})
	t = ceil(Int,0.1 * R)

	current_iterations = 0
	while true
		println(current_iterations)
		last_point = group_points[rand(1:length(group_points))]
		xl = (last_point[1] - t) >= 0 ? last_point[1] - t : 0
		xr = (last_point[1] + t) <= canvas.x ? last_point[1] + t : canvas.x
		yb = (last_point[2] - t) >= 0 ? last_point[2] - t : 0
		yt = (last_point[2] + t) <= canvas.y ? last_point[2] + t : canvas.y

		x = rand(xl : xr) # A point inside the radious
		y = rand(yb : yt)

		d = evaluate(Euclidean(),collect(last_point),[x,y])
		# println("Point ($x,$y) -- distance $d -- t: $t")
		if d <= t && !((x,y) in cluster.points) && !is_inside_another_cluster(canvas,(x,y))
			return (x,y)
		end
		current_iterations +=1
		if current_iterations > MAX_ITERATION
			throw(ErrorException("Unable to plot point, the graph is too dense.
						Try reducing trasmission range or reduce waypoints number"))
		end
	end
end

function generate_group_first_point(canvas::Canvas,last_point::Tuple{Int64,Int64},cluster::Cluster)
	# Y = round(Int,(2 * canvas.x * canvas.y) / canvas.num_cluster)
	Y = round(Int,(canvas.x ) / canvas.num_cluster)

	xl = (last_point[1] - Y/3) >= 0 ? last_point[1] - Y/3 : 0
	xr = (last_point[1] + Y/3) <= canvas.x ? last_point[1] + Y/3 : canvas.x
	yb = (last_point[2] - Y/3) >= 0 ? last_point[2] - Y/3 : 0
	yt = (last_point[2] + Y/3) <= canvas.y ? last_point[2] + Y/3 : canvas.y

	current_iterations = 0
	while true
		println(current_iterations)
		x = ceil(Int,rand(xl : xr)) # A point inside the radious
		y = ceil(Int,rand(yb : yt))

		d = evaluate(Euclidean(),collect(last_point),[x,y])
		# println("Point ($x,$y) -- distance $d -- Y/3:$(Y/3) -- Y/4:$(Y/4) -- Y: $Y")
		if d >= Y/4 && d <= Y/3 && !((x,y) in cluster.points) && !is_inside_another_cluster(canvas,(x,y))
			return (x,y)
		end
		current_iterations +=1
		if current_iterations > MAX_ITERATION
			throw(ErrorException("Unable to plot point, the graph is too dense.
						Try reducing trasmission range or reduce waypoints number"))
		end
	end
end

# Get a poin in the Canvas

function generate_point(canvas::Canvas)
	current_iterations = 0
	while true
		println(current_iterations)
		p = (rand(0:canvas.x) , rand(0:canvas.y))
		if !is_inside_another_cluster(canvas,p)
			return p
		end
		current_iterations +=1
		if current_iterations > MAX_ITERATION
			throw(ErrorException("Unable to plot point, the graph is too dense.
						Try reducing trasmission range or reduce waypoints number"))
		end
	end

end

function is_inside_another_cluster(canvas::Canvas, point::Tuple)
	if length(canvas.clusters) == 0
		return false
	else
		for i in 1:length(canvas.clusters)
			# println("Checking point $point inside $(canvas.clusters[i].boundary)")
			if is_inside_cluster(canvas.clusters[i], point)
				# print("Is inside another is_inside_another_cluster")
				return true
			end
		end
		return false
	end
end


text:
# Code for load-balancing optimization

using Random, Distributions

function setup_loadbalancer_channels(;chin_size=200, chout_size=200)
    chin = ConcurrentCircularBuffer{ProfilerOutput}(chin_size)
    chout = ConcurrentCircularBuffer{@NamedTuple{worker::Int,p::Int}}(chout_size)
    chin, chout
end

"""

For each worker, return a lower bound on the probability of this worker being the fastest of all
workers, i.e., when `nwait < nworkers`, `rv[i]` is a lower bound on `log(ls[i])`.
"""
function fastest_lower_bound!(rv, dzs, dys)
    length(rv) == length(dzs) || throw(DimensionMismatch("rv has dimension $(length(rv)), but dzs has dimension $(length(dzs))"))
    length(rv) == length(dys) || throw(DimensionMismatch("rv has dimension $(length(rv)), but dys has dimension $(length(dys))"))
    n = length(rv)

    # compute the midpoint of the means
    cz = mean(mean, dzs)
    cy = mean(mean, dys)

    # for each i, prob. of (zi <= cz and all others > cz) and (yi <= cy and all others > cy)
    pz, py = 0.0, 0.0
    for i in 1:n
        pz += logccdf(dzs[i], cz)
        py += logccdf(dys[i], cy)
    end
    for i in 1:n
        vz = logccdf(dzs[i], cz)
        vy = logccdf(dys[i], cy)
        pz -= vz
        py -= vy
        rv[i] = max(rv[i], pz + logcdf(dzs[i], cz) + py + logcdf(dys[i], cy))
        pz += vz
        py += vy
    end
    rv
end

fastest_lower_bound(dzs, dys) = fastest_lower_bound!(fill(-Inf, length(dzs)), dzs, dys)

"""

For each worker, return a lower bound on the probability of this worker being the slowest of all
workers.
"""
function slowest_lower_bound!(rv, dzs, dys)
    length(rv) == length(dzs) || throw(DimensionMismatch("rv has dimension $(length(rv)), but dzs has dimension $(length(dzs))"))
    length(rv) == length(dys) || throw(DimensionMismatch("rv has dimension $(length(rv)), but dys has dimension $(length(dys))"))
    n = length(rv)

    # compute the midpoint of the means
    cz = mean(mean, dzs)
    cy = mean(mean, dys)

    # for each i, prob. of (zi > cz and all others <= cz) and (yi > cy and all others <= cy)
    pz, py = 0.0, 0.0
    for i in 1:n
        pz += logcdf(dzs[i], cz)
        py += logcdf(dys[i], cy)
    end
    for i in 1:n
        vz = logcdf(dzs[i], cz)
        vy = logcdf(dys[i], cy)
        pz -= vz
        py -= vy
        rv[i] = min(rv[i], pz + logccdf(dzs[i], cz) + py + logccdf(dys[i], cy))
        pz += vz
        py += vy
    end
    rv
end

slowest_lower_bound(dzs, dys) = slowest_lower_bound!(fill(Inf, length(dzs)), dzs, dys)

"""

Setup the simulator and simulate the fraction of iterations each worker participates in.
"""
function simulate!(ls, ps; sim, θs, comp_mcs, comp_vcs, simulation_nsamples, simulation_niterations)
    nworkers = length(ls)

    # setup compute latency distributions
    for i in 1:nworkers
        m = comp_mcs[i] * θs[i] / ps[i]
        v = comp_vcs[i] * (θs[i] / ps[i])^2
        sim.comp_distributions[i] = CodedComputing.distribution_from_mean_variance(Gamma, m, v)
    end

    # run nsamples simulations, each consisting of niterations steps
    latency = 0.0
    ls .= 0
    for _ in 1:simulation_nsamples
        empty!(sim)
        step!(sim, simulation_niterations)
        latency += sim.time / simulation_niterations
        ls .+= sim.nfresh ./ simulation_niterations
    end
    latency /= simulation_nsamples
    ls ./= simulation_nsamples
    ls .= log.(ls)
    fastest_lower_bound!(ls, sim.comp_distributions, sim.comm_distributions)
    latency, ls
end

"""

Return the expected latency of the `i`-th worker.
"""
expected_worker_latency(i; θs, ps, comp_mcs, comm_mcs) = comp_mcs[i] * θs[i] / ps[i] + comm_mcs[i]

function min_max_expected_worker_latency(;θs, ps, comp_mcs, comm_mcs)
    nworkers = length(ps)
    vmin, vmax = Inf, -Inf
    for i in 1:nworkers
        v = expected_worker_latency(i; θs, ps, comp_mcs, comm_mcs)
        vmin = min(vmin, v)
        vmax = max(vmax, v)
    end
    vmin, vmax
end

function optimize!(ps::AbstractVector, ps_prev::AbstractVector, sim::EventDrivenSimulator; ls=zeros(length(ps)), contribs=zeros(length(ps)), θs, comp_mcs, comp_vcs, comm_mcs, comm_vcs, min_contribution::Real=(sim.nwait / length(ps_prev) / mean(ps_prev)), min_comp_fraction::Real=0.1, simulation_niterations::Integer=50, simulation_nsamples::Integer=3)
    nworkers = length(ps)
    0 < min_contribution <= 1 || throw(ArgumentError("min_contribution is $min_contribution"))
    length(ps_prev) == nworkers || throw(DimensionMismatch("ps_prev has dimension $(length(ps_prev)), but nworkers is $nworkers"))
    length(θs) == nworkers || throw(DimensionMismatch("θs has dimension $(length(θs)), but nworkers is $nworkers"))
    length(comp_mcs) == nworkers || throw(DimensionMismatch("comp_mcs has dimension $(length(comp_mcs)), but nworkers is $nworkers"))
    length(comp_vcs) == nworkers || throw(DimensionMismatch("comp_vcs has dimension $(length(comp_vcs)), but nworkers is $nworkers"))
    length(comm_mcs) == nworkers || throw(DimensionMismatch("comm_mcs has dimension $(length(comm_mcs)), but nworkers is $nworkers"))
    length(comm_vcs) == nworkers || throw(DimensionMismatch("comm_vcs has dimension $(length(comm_vcs)), but nworkers is $nworkers"))
    0 < min_comp_fraction < 1 || throw(ArgumentError("min_comp_fraction is $min_comp_fraction"))

    # setup communication latency distributions
    for i in 1:nworkers
        m = comm_mcs[i]
        v = comm_vcs[i]
        sim.comm_distributions[i] = CodedComputing.distribution_from_mean_variance(Gamma, m, v)
    end

    # loss for previous solution for reference
    latency0, ls = simulate!(ls, ps_prev; sim, θs, comp_mcs, comp_vcs, simulation_nsamples=2*simulation_nsamples, simulation_niterations=2*simulation_nsamples)
    contribs .= ls .+ log.(θs) .- log.(ps_prev)
    contrib0 = sum(exp, contribs)
    vmin, vmax = min_max_expected_worker_latency(;θs, ps=ps_prev, comp_mcs, comm_mcs)
    loss0 = vmax / vmin

    # equalize latency between workers
    ## find the slowest worker
    i = 0
    v = -Inf
    for j in 1:nworkers
        comp_latency = comp_mcs[j] * θs[j] / ps[j]
        if comp_latency + comm_mcs[j] > v
            i = j
            v = comp_latency + comm_mcs[j]
        end
    end

    ## increase the workload of all other workers, so that what was the slowest worker becomes the fastest
    for j in 1:nworkers
        if j == i
            continue
        end
        ps[j] = max(1.0, round(comp_mcs[j] * θs[j] / (v - comm_mcs[j])))
    end

    # latency and contribution of the current solution
    latency, ls = simulate!(ls, ps; sim, θs, comp_mcs, comp_vcs, simulation_nsamples=2*simulation_nsamples, simulation_niterations=2*simulation_niterations)
    contribs .= ls .+ log.(θs) .- log.(ps)
    contrib = sum(exp, contribs)

    # while below the contribution constraint, assign more work to the fastest workers
    while contrib < min_contribution

        # find the fastest worker with at least 2 sub-partitions
        i = 0
        v = Inf
        for j in 1:nworkers
            comp_latency = comp_mcs[j] * θs[j] / ps[j]
            if ps[j] >= 2 && (comp_latency + comm_mcs[j]) < v
                i = j
                v = comp_latency + comm_mcs[j]
            end
        end
        if iszero(i)
            break
        end
        δ = max(1, round(ps[i] * 0.01))
        ps[i] -= δ

        latency, ls = simulate!(ls, ps; sim, θs, comp_mcs, comp_vcs, simulation_nsamples, simulation_niterations)
        contribs .= ls .+ log.(θs) .- log.(ps)
        contrib = sum(exp, contribs)

        # double check the exit condition
        if contrib >= min_contribution
            latency, ls = simulate!(ls, ps; sim, θs, comp_mcs, comp_vcs, simulation_nsamples=2*simulation_nsamples, simulation_niterations=2*simulation_niterations)
            contribs .= ls .+ log.(θs) .- log.(ps)
            contrib = sum(exp, contribs)
        end
    end

    # while within the contribution constraint, reduce the workload of slow workers
    while contrib >= min_contribution * 0.99

        # find the slowest worker with a comp. latency that accounts for at least
        # min_comp_fraction of the overall latency of the worker
        # (to ensure that the load-balancer doesn't push comp. latency to zero)
        i = 0
        v = -Inf
        for j in 1:nworkers
            comp_latency = comp_mcs[j] * θs[j] / ps[j]
            if comp_latency / (comp_latency + comm_mcs[j]) >= min_comp_fraction && (comp_latency + comm_mcs[j]) > v
                i = j
                v = comp_latency + comm_mcs[j]
            end
        end
        if iszero(i)
            break
        end
        δ = max(1, round(ps[i] * 0.01))
        ps[i] += δ
        latency, ls = simulate!(ls, ps; sim, θs, comp_mcs, comp_vcs, simulation_nsamples, simulation_niterations)
        contribs .= ls .+ log.(θs) .- log.(ps)
        contrib = sum(exp, contribs)

        # double check the exit condition
        if contrib < min_contribution * 0.99
            latency, ls = simulate!(ls, ps; sim, θs, comp_mcs, comp_vcs, simulation_nsamples=2*simulation_nsamples, simulation_niterations=2*simulation_niterations)
            contribs .= ls .+ log.(θs) .- log.(ps)
            contrib = sum(exp, contribs)
        end
    end

    vmin, vmax = min_max_expected_worker_latency(;θs, ps, comp_mcs, comm_mcs)
    loss = vmax / vmin

    ps, latency0, contrib0, loss0, latency, contrib, loss
end

function load_balancer(chin::ConcurrentCircularBuffer, chout::ConcurrentCircularBuffer; nwait::Integer, nworkers::Integer, nsubpartitions::Union{Integer,<:AbstractVector}, min_improvement::Real=0.9)
    0 < nworkers || throw(ArgumentError("nworkers is $nworkers"))
    0 < nwait <= nworkers || throw(ArgumentError("nwait is $nwait, but nworkers is $nworkers"))
    if typeof(nsubpartitions) <: Integer
        0 < nsubpartitions || throw(ArgumentError("nsubpartitions must be positive, but is $nsubpartitions"))
    else
        length(nsubpartitions) == nworkers || throw(DimensionMismatch("nsubpartitions has dimension $(length(nsubpartitions)), but nworkers is $nworkers"))
    end
    @info "load_balancer task started on thread $(Threads.threadid())"

    # fraction of dataset stored by and fraction of local data processed per iteration for each worker
    θs = fill(1/nworkers, nworkers)
    ps = typeof(nsubpartitions) <: Integer ? fill(float(nsubpartitions), nworkers) : float.(nsubpartitions)
    ps_prev = copy(ps)

    # minimum contribution set to that of the initial solution
    min_contribution = nwait / length(ps_prev) / mean(ps_prev)

    # mean and variance coefficients for each worker
    # (dummy values used to force pre-compilation of the optimizer)
    comp_mcs = ones(nworkers)
    comp_vcs = ones(nworkers)
    comm_mcs = ones(nworkers)
    comm_vcs = ones(nworkers)

    # buffers used by the optimizer
    ls = zeros(nworkers)
    contribs = zeros(nworkers)

    # reusable simulator
    comp_distributions = [Gamma() for _ in 1:nworkers]
    comm_distributions = [Gamma() for _ in 1:nworkers]
    sim = EventDrivenSimulator(;nwait, nworkers, comp_distributions, comm_distributions)

    # run the optimizer once to force pre-compilation
    optimize!(ps, ps_prev, sim; ls, contribs, θs, min_contribution, comp_mcs, comp_vcs, comm_mcs, comm_vcs)

    # reset changes made by the optimizer
    ps .= nsubpartitions
    ps_prev .= nsubpartitions
    comp_mcs .= NaN
    comp_vcs .= NaN
    comm_mcs .= NaN
    comm_vcs .= NaN

    # helper to check if there is any missing latency data
    # (in which case we shouldn't run the load-balancer)
    all_populated = false
    function check_populated()
        rv = iszero(count(isnan, comp_mcs))
        if rv
            rv = rv && iszero(count(isnan, comp_vcs))
        end
        if rv
            rv = rv && iszero(count(isnan, comm_mcs))
        end
        if rv
            rv = rv && iszero(count(isnan, comm_vcs))
        end
        if rv
            rv = rv && iszero(count(isnan, θs))
        end
        if rv
            rv = rv && iszero(count(isnan, ps))
        end
        rv
    end

    # consume incoming ProfilerOutput samples and run the optimizer
    while isopen(chin)

        # consume all values currently in the channel
        try
            vin = take!(chin)
            isnan(vin.comp_mc) || (comp_mcs[vin.worker] = vin.comp_mc)
            isnan(vin.comp_vc) || (comp_vcs[vin.worker] = vin.comp_vc)
            isnan(vin.comm_mc) || (comm_mcs[vin.worker] = vin.comm_mc)
            isnan(vin.comm_vc) || (comm_vcs[vin.worker] = vin.comm_vc)
        catch e
            if e isa InvalidStateException
                @info "error taking value from input channel" e
                break
            else
                rethrow()
            end
        end
        while isready(chin)
            try
                vin = take!(chin)
                isnan(vin.comp_mc) || (comp_mcs[vin.worker] = vin.comp_mc)
                isnan(vin.comp_vc) || (comp_vcs[vin.worker] = vin.comp_vc)
                isnan(vin.comm_mc) || (comm_mcs[vin.worker] = vin.comm_mc)
                isnan(vin.comm_vc) || (comm_vcs[vin.worker] = vin.comm_vc)
            catch e
                if e isa InvalidStateException
                    @info "error taking value from input channel" e
                    break
                else
                    rethrow()
                end
            end
        end

        # verify that we have complete latency information for all workers
        if !all_populated
            all_populated = check_populated()
            if all_populated
                @info "load_balancer has latency information for all workers; starting to optimize"
            end
        end
        if !all_populated
            continue
        end

        try
            # @info "load_balancer optimization started with ps: $ps, θs: $θs, comp_mcs: $comp_mcs, comp_vcs: $comp_vcs, comm_mcs: $comm_mcs, comm_vcs: $comm_vcs"
            # @info "load_balancer optimization started with ps = $ps, θs = $θs, comp_mcs = $comp_mcs, comp_vcs = $comp_vcs, comm_mcs = $comm_mcs, comm_vcs = $comm_vcs"
            # @info "ps_prev: $ps_prev, ps_baseline: $ps_baseline"
            # @info "load_balancer optimization started"
            t = @elapsed begin
                ps, latency0, contrib0, loss0, latency, contrib, loss = optimize!(ps, ps_prev, sim; ls, contribs, θs, min_contribution, comp_mcs, comp_vcs, comm_mcs, comm_vcs)
            end

            # compare the initial and new solutions, and continue if the change isn't large enough
            if isnan(loss) || isinf(loss) || (loss / loss0) > min_improvement
                # @info "load-balancer finished in $(t) seconds with loss $loss and loss0 $loss0; continuing. ps: $ps"
                continue
            end
            # @info "load-balancer finished in $(t) seconds with loss $loss and loss0 $loss0; accepting it. ps: $ps"
            ps_prev .= ps

            # push any changes into the output channel
            for i in 1:nworkers
                p = max(1, round(Int, ps[i]))
                vout = @NamedTuple{worker::Int,p::Int}((i, p))
                try
                    push!(chout, vout)
                catch e
                    if e isa InvalidStateException
                        @info "error pushing value into output channel" e
                        break
                    else
                        rethrow()
                    end
                end
            end
        catch e
            if e isa ArgumentError
                @error "load-balancer failed; trying again" e
            else
                rethrow()
            end
        end
    end
    @info "load_balancer task finished"
    return
end


text:
using CfGRIB
using FileIO

#  Goal is for the backends to be optional, so only include the files if the
#  user has the required backend installed
BACKEND_PROVIDERS = [:AxisArrays, :DimensionalData]
AVAILABLE_BACKENDS = []
DEFAULT_BACKEND = DataSet

"""
Backends are wrappers around named array data types, contain information on the
dimensions, attributes, and encoding of the data, as well as the datasets
themselves.

# See also

[`AxisArrayWrapper`](@ref AxisArrayWrapper),
[`DimensionalArrayWrapper`](@ref DimensionalArrayWrapper)
"""
abstract type Backend end

"""
Abstract type for backends which also wrap data types provided by external
packages, instead of just providing a 'vanilla' implementation of the data type.

Python's `xarray` has two fundamental data structures, a `DataArray` which holds
the data for a **single** multi-dimensional variable and its coordinates, as
well as a `Dataset` which holds **multiple DataArrays** that (typically) share
the same coordinates.

# See also

[`AxisArrayWrapper`](@ref AxisArrayWrapper),
[`DimensionalArrayWrapper`](@ref DimensionalArrayWrapper)
"""
abstract type ArrayWrapper <: Backend end

try
    using DimensionalData
    include("./backends/dimensionaldata.jl")
    global DEFAULT_BACKEND = DimensionalArrayWrapper
    append!(AVAILABLE_BACKENDS, [:DimensionalArray])
catch e
    if !(e isa ArgumentError)
        throw(e)
    end
end

try
    using AxisArrays
    include("./backends/axisarrays.jl")
    global DEFAULT_BACKEND = AxisArrayWrapper
    append!(AVAILABLE_BACKENDS, [:AxisArray])
catch e
    if !(e isa ArgumentError)
        throw(e)
    end
end#

#  If no backends could be loaded print a warning
if DEFAULT_BACKEND == DataSet
    @warn """No backends could be loaded, setting default to bare `DataSet`.
    Please include one of: $BACKEND_PROVIDERS.
    Run `import Pkg; Pkg.add("BACKEND")` to install one of the array backends"""
end

# function load(f::String; backend=DEFAULT_BACKEND, kwargs...)
#     ds = DataSet(f, kwargs...)

#     return convert(backend, ds)
# end

function fileio_load(f::File{format"GRIB"}; backend=DEFAULT_BACKEND, kwargs...)
    """ Load command for FileIO integration. Following line must be in the
    FileIO Registry for this to function properly:

        add_format(format"GRIB", UInt8[0x47, 0x52, 0x49, 0x42], ".grib", [:CfGRIB])
    """
    ds = DataSet(f.filename, kwargs...)
    return convert(backend, ds)
end


text:
module TestFit

using Test
using Dates
using LightGBM

# test fixtures
train_matrix = rand(10000,70) # create random dataset
train_labels = rand([0, 1], 10000)
train_dataset = LightGBM.LGBM_DatasetCreateFromMat(train_matrix, "")

test_matrix = rand(5000,70) # create random dataset
test_labels = rand([0, 1], 5000)
test_dataset = LightGBM.LGBM_DatasetCreateFromMat(test_matrix, "", train_dataset)

test2_matrix = rand(2500,70) # create second random dataset
test2_labels = rand([0, 1], 2500)
test2_dataset = LightGBM.LGBM_DatasetCreateFromMat(test2_matrix, "", train_dataset)

@testset "test fit! with dataset -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(objective = "binary", num_class = 1)
    LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)
    LightGBM.LGBM_DatasetSetField(test_dataset, "label", test_labels)

    # Act
    LightGBM.fit!(estimator, train_dataset, test_dataset, verbosity = -1);
    p_binary = LightGBM.predict(estimator, train_matrix, verbosity = -1)
    binary_classes = LightGBM.predict_classes(estimator, train_matrix, verbosity = -1)

    # Assert
    # check also that binary outputs and classes are invariant in prediction output shapes, i.e. they're matrices
    @test length(size(p_binary)) == 2
    @test length(size(binary_classes)) == 2

end

@testset "test fit! with dataset without testset -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(objective = "binary", num_class = 1)
    LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)

    # Act
    LightGBM.fit!(estimator, train_dataset, verbosity = -1);
    p_binary = LightGBM.predict(estimator, train_matrix, verbosity = -1)
    binary_classes = LightGBM.predict_classes(estimator, train_matrix, verbosity = -1)

    # Assert
    # check also that binary outputs and classes are invariant in prediction output shapes, i.e. they're matrices
    @test length(size(p_binary)) == 2
    @test length(size(binary_classes)) == 2

end

@testset "test train! single test set -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        objective = "binary",
        num_class = 1,
        is_training_metric = true,
        metric = ["auc"],
    )

    LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)
    LightGBM.LGBM_DatasetSetField(test_dataset, "label", test_labels)

    bst_parameters = LightGBM.stringifyparams(estimator, LightGBM.BOOSTERPARAMS) * " verbosity=-1"
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)

    # Act
    results = LightGBM.train!(estimator, 10, ["test_set_with_whatever_label"], -1, DateTime(2020, 1, 11, 11, 00));

    # Assert
    @test length(results["training"]["auc"]) > 0
    @test length(results["test_set_with_whatever_label"]["auc"]) > 0

end

@testset "test train! multiple test set multiple metrics -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        objective = "binary",
        num_class = 1,
        is_training_metric = true,
        metric = ["auc", "l2"],
    )

    LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)
    LightGBM.LGBM_DatasetSetField(test_dataset, "label", test_labels)
    LightGBM.LGBM_DatasetSetField(test2_dataset, "label", test2_labels)

    bst_parameters = LightGBM.stringifyparams(estimator, LightGBM.BOOSTERPARAMS) * " verbosity=-1"
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test2_dataset)

    # Act
    results = LightGBM.train!(estimator, 10, ["test_1", "test_2"], -1, DateTime(2020, 1, 11, 11, 00));

    # Assert
    @test length(results["training"]["auc"]) > 0
    @test length(results["training"]["l2"]) > 0
    @test length(results["test_1"]["auc"]) > 0
    @test length(results["test_1"]["l2"]) > 0
    @test length(results["test_2"]["auc"]) > 0
    @test length(results["test_2"]["l2"]) > 0

end

@testset "test train! multiple test set multiple metrics -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(
        objective = "binary",
        num_class = 1,
        is_training_metric = true,
        metric = ["auc", "l2"],
    )

    LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)
    LightGBM.LGBM_DatasetSetField(test_dataset, "label", test_labels)
    LightGBM.LGBM_DatasetSetField(test2_dataset, "label", test2_labels)

    bst_parameters = LightGBM.stringifyparams(estimator, LightGBM.BOOSTERPARAMS) * " verbosity=-1"
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test2_dataset)

    # Act
    results = LightGBM.train!(estimator, 10, ["test_1", "test_2"], -1, DateTime(2020, 1, 11, 11, 00));

    # Assert
    @test length(results["training"]["auc"]) > 0
    @test length(results["training"]["l2"]) > 0
    @test length(results["test_1"]["auc"]) > 0
    @test length(results["test_1"]["l2"]) > 0
    @test length(results["test_2"]["auc"]) > 0
    @test length(results["test_2"]["l2"]) > 0

end

@testset "test store_scores! adds results" begin
    # Arrange
    results_fixture = Dict{String,Dict{String,Vector{Float64}}}()

    # Act
    LightGBM.store_scores!(results_fixture, "whatever_dataset","some_metric",0.5)
    LightGBM.store_scores!(results_fixture, "whatever_dataset","some_metric",0.55)

    LightGBM.store_scores!(results_fixture, "whatever_dataset","some_extra_metric",0.6)
    LightGBM.store_scores!(results_fixture, "whatever_dataset","some_extra_metric",0.65)

    LightGBM.store_scores!(results_fixture, "whatever_second_dataset","some_metric",0.7)
    LightGBM.store_scores!(results_fixture, "whatever_second_dataset","some_metric",0.75)

    # Assert
    @test results_fixture["whatever_dataset"]["some_metric"] == [0.5, 0.55]
    @test results_fixture["whatever_dataset"]["some_extra_metric"] == [0.6, 0.65]
    @test results_fixture["whatever_second_dataset"]["some_metric"] == [0.7, 0.75]

end

@testset "test eval_metrics! early stop disabled" begin
    # Arrange
    results_fixture = Dict{String,Dict{String,Vector{Float64}}}()
    estimator = LightGBM.LGBMClassification(
        num_iterations = 10, objective = "binary", num_class = 1,
        is_training_metric = false, metric = ["auc"],
        early_stopping_round = 0 # default value, but stating explicitly to test!
    )
    LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)
    LightGBM.LGBM_DatasetSetField(test_dataset, "label", test_labels)

    bst_parameters = LightGBM.stringifyparams(estimator, LightGBM.BOOSTERPARAMS) * " verbosity=-1"
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)

    bigger_is_better = Dict{String,Float64}("auc" => 1.0)
    best_score = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))
    best_iter = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))

    # Act and assert each round returns false (despite iter 1 have best score always)
    for iter in 1:10
        LightGBM.LGBM_BoosterUpdateOneIter(estimator.booster)
        output = LightGBM.eval_metrics!(
            results_fixture, estimator, ["test_bla"], iter, -1,
            bigger_is_better, best_score, best_iter, ["auc"]
        )

        @test output == false

    end
end

#=
This test assumes:
* early_stopping_round set to 5
* total iterations = 10
* best iter is ALWAYS round 1 of metric score 1 (auc)
Criteria: early_stopping should kick in on round 6
=#
@testset "test eval_metrics! stops correctly" begin
    # Arrange
    results_fixture = Dict{String,Dict{String,Vector{Float64}}}()
    estimator = LightGBM.LGBMClassification(
        num_iterations = 10, objective = "binary", num_class = 1,
        is_training_metric = false, metric = ["auc"],
        early_stopping_round = 5
    )
    LightGBM.LGBM_DatasetSetField(train_dataset, "label", train_labels)
    LightGBM.LGBM_DatasetSetField(test_dataset, "label", test_labels)

    bst_parameters = LightGBM.stringifyparams(estimator, LightGBM.BOOSTERPARAMS) * " verbosity=-1"
    estimator.booster = LightGBM.LGBM_BoosterCreate(train_dataset, bst_parameters)
    LightGBM.LGBM_BoosterAddValidData(estimator.booster, test_dataset)

    bigger_is_better = Dict{String,Float64}("auc" => 1.0)
    best_score = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))
    best_iter = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))

    # Act and assert each round returns expected output
    for iter in 1:10
        LightGBM.LGBM_BoosterUpdateOneIter(estimator.booster)
        output = LightGBM.eval_metrics!(
            results_fixture, estimator, ["test_bla"], iter, -1,
            bigger_is_better, best_score, best_iter, ["auc"]
        )

        # reset scores to round 1 being best
        best_score = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))
        best_iter = Dict{String,Dict{String,Real}}("auc" => Dict("test_bla" => 1))

        if iter < 6
            @test output == false
        else
            @test output == true
        end
    end
end

end # module


text:
"""
            Functions to read a snapshot and store it in an object.

            !!!! BROKEN AND NOT IN USE !!!!

            only head_to_obj can be used!

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2018-12-12

"""

# function snap_to_obj(filename)
#
#     head = head_to_obj(filename)
#
#     if head.snap_format == 2
#
#        P = snap_2_p(filename, head)
#
#     else
#
#        P = snap_1_p(filename, head)
#
#     end
#
#     return P
#
# end

# function snap_2_p(filename, head)
#
#         f = open(filename)
#
#         seek(f, 296)
#
#         N = sum(head.npart)
#         skipsize = read(f, Int32, 1)[1]
#         bit_size = skipsize/(3*N)
#
#         if Int(bit_size) == 4
#             P = Array{part_single,1}(0)
#             # set up array of particles
#
#             for i = 1:6
#                  if head.npart[i] != Int32(0)
#                      for k = 1:head.npart[i]
#                          push!(P,part_single(i-1))
#                      end
#                  end
#             end
#
#         elseif Int(bit_size) == 8
#
#                 P = Array{part_double,1}(0)
#                 for i = 1:6
#                      if head.npart[i] != Int32(0)
#                          for k = 1:head.npart[i]
#                              push!(P,part_double(i-1))
#                          end
#                      end
#                 end
#         else
#             println("read error! neither 32 nor 64 bits data!")
#             return -1
#         end
#
#         # read positions
#         if Int(bit_size) == 4
#             for i = 1:sum(head.npart)
#                  P[i].pos = read(f, Float32, 3)
#             end
#         else
#             for i = 1:sum(head.npart)
#                  P[i].pos = read(f, Float64, 3)
#             end
#         end
#
#
#         p = position(f)
#
#         # skip identifiers
#         seek(f, p+24)
#
#         if Int(bit_size) == 4
#             dummy = read(f, Float32, (3,sum(head.npart)))
#         else
#             dummy = read(f, Float64, (3,sum(head.npart)))
#         end
#
#         # read Velocities
#         for i = 1:sum(head.npart)
#
#             P[i].vel = dummy[:,i]
#
#         end
#
#         p = position(f)
#
#         # skip identifiers
#         seek(f, p+24)
#
#         for i = 1:sum(head.npart)
#
#             P[i].id = read(f, UInt32)
#
#         end
#
#         p = position(f)
#
#         # skip identifiers
#         seek(f, p+24)
#
#         # read masses
#         for k = 1:6
#             if head.massarr[k] == Int32(0)
#                  for i = 1:head.npart[k]
#                      if bit_size == 4
#                         P[i].m = read(f, Float32)
#                      else
#                         P[i].m = read(f, Float64)
#                      end
#                  end
#             else
#                  for i = 1:head.npart[k]
#                     P[i].m = head.massarr[k]
#                  end
#             end
#         end
#
#
#         if head.npart[1] != Int32(0)
#
#             # skip identifiers
#             p = position(f)
#             seek(f, p+24)
#
#             # read U
#             for i = 1:head.npart[1]
#
#                if bit_size == 4
#                   P[i].U = read(f, Float32)
#                else
#                   P[i].U = read(f, Float64)
#                end
#
#            end
#
#             # skip identifiers
#             p = position(f)
#             seek(f, p+24)
#
#             # Read Density
#             for i = 1:head.npart[1]
#
#                if bit_size == 4
#                   P[i].rho = read(f, Float32)
#                else
#                   P[i].rho = read(f, Float64)
#                end
#
#            end
#
#             # skip identifiers
#             p = position(f)
#             seek(f, p+24)
#
#             # read smoothing length
#             for i = 1:head.npart[1]
#
#                if bit_size == 4
#                   P[i].sml = read(f, Float32)
#                else
#                   P[i].sml = read(f, Float64)
#                end
#
#            end
#
#
#         end
#
#         close(f)
#
#         return P
#
# end
#
# function snap_1_p(filename, head)
#
#     f = open(filename)
#
#     seek(f, 264)
#
#     N = sum(head.npart)
#     skipsize = read(f, Int32, 1)[1]
#     bit_size = skipsize/(3*N)
#
#     if Int(bit_size) == 4
#         P = Array{part_single,1}(0)
#         # set up array of particles
#
#         for i = 1:6
#              if head.npart[i] != Int32(0)
#                  for k = 1:head.npart[i]
#                      push!(P,part_single(i-1))
#                  end
#              end
#         end
#
#     elseif Int(bit_size) == 8
#
#             P = Array{part_double,1}(0)
#             for i = 1:6
#                  if head.npart[i] != Int32(0)
#                      for k = 1:head.npart[i]
#                          push!(P,part_double(i-1))
#                      end
#                  end
#             end
#     else
#         println("read error! neither 32 nor 64 bits data!")
#         return -1
#     end
#
#     # read positions
#     if Int(bit_size) == 4
#         for i = 1:sum(head.npart)
#              P[i].pos = read(f, Float32, 3)
#         end
#     else
#         for i = 1:sum(head.npart)
#              P[i].pos = read(f, Float64, 3)
#         end
#     end
#
#
#     p = position(f)
#
#     # skip identifiers
#     seek(f, p+8)
#
#     if Int(bit_size) == 4
#         dummy = read(f, Float32, (3,sum(head.npart)))
#     else
#         dummy = read(f, Float64, (3,sum(head.npart)))
#     end
#
#     # read Velocities
#     for i = 1:sum(head.npart)
#
#         P[i].vel = dummy[:,i]
#
#     end
#
#     p = position(f)
#
#     # skip identifiers
#     seek(f, p+8)
#
#     for i = 1:sum(head.npart)
#
#         P[i].id = read(f, UInt32)
#
#     end
#
#     p = position(f)
#
#     # skip identifiers
#     seek(f, p+8)
#
#     # read masses
#     for k = 1:6
#         if head.massarr[k] == Int32(0)
#              for i = 1:head.npart[k]
#                  if bit_size == 4
#                     P[i].m = read(f, Float32)
#                  else
#                     P[i].m = read(f, Float64)
#                  end
#              end
#         else
#              for i = 1:head.npart[k]
#                 P[i].m = head.massarr[k]
#              end
#         end
#     end
#
#
#     if head.npart[1] != Int32(0)
#
#         # skip identifiers
#         p = position(f)
#         seek(f, p+8)
#
#         # read U
#         for i = 1:head.npart[1]
#
#            if bit_size == 4
#               P[i].U = read(f, Float32)
#            else
#               P[i].U = read(f, Float64)
#            end
#
#        end
#
#         # skip identifiers
#         p = position(f)
#         seek(f, p+8)
#
#         # Read Density
#         for i = 1:head.npart[1]
#
#            if bit_size == 4
#               P[i].rho = read(f, Float32)
#            else
#               P[i].rho = read(f, Float64)
#            end
#
#        end
#
#         # skip identifiers
#         p = position(f)
#         seek(f, p+8)
#
#         # read smoothing length
#         for i = 1:head.npart[1]
#
#            if bit_size == 4
#               P[i].sml = read(f, Float32)
#            else
#               P[i].sml = read(f, Float64)
#            end
#
#        end
#
#
#     end
#
#     close(f)
#
#     return P
#
# end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
module BQCESubroutine

export broutine!, @broutine, log2dim, log2i

using MLStyle
using Expronicon
using ExprTools
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, AbstractLocations, merge_locations
using LinearAlgebra
using LoopVectorization
using StrideArrays
using StrideArrays: StrideArray
using Polyester: @batch
using ArrayInterface

include("utils2.jl")
include("schedule.jl")
include("subspace/bit.jl")
include("mul/mul.jl")
include("routine/bit.jl")
include("threading.jl")

# include("utils.jl")
# using .Utilities

include("codegen/utils.jl")
include("codegen/broutine.jl")
# include("mul.jl")
include("statevector.jl")

end


text:
##An ABCDistance type defines which distance to use in an ABC algorithm.
##This file defines several concrete subtypes and two methods:
##  *  "init" initialises the distance based on simulated summary statistics.
##     These are sometimes used to initialise weights.
##     Alternatively this method can have no effect.
##     (nb infinite summary statistic values should be ignored here. They represent pre-rejected simulations.)
##  *  "evaldist" evaluates the distance for particular vector of simulated summaries.
##An ABCDistance object must be created for a particular vector of observed summaries.
##Ths is computationally convenient for some distance functions.
abstract type ABCDistance end

#################################
##Subtypes and methods of ABCDistance
#################################
"
`ABCDistance` objects are constructed specifying the observed summaries `sobs`.
Some types require further information which is initialised during an ABC algorithm.
This is set up by the `init` method.
Its arguments are

* `x` The `ABCDistance` object
* `sumstats` All the simulated summaries from the current iteration (including rejections). This is a `Array{Float64, 2}` whose columns are the simulations.
* `parameters` The parameters corresponding to the `sumstats` simulations. This is also a `Array{Float64, 2}` whose columns are the simulations.

`abcRejection` always calls this method.
So does `abcPMC` on the first iteration, or all iterations for `adaptive=true`.
`abcPMC_comparison` calls this on the first iteration if `initialise_dist=true` or otherwise not at all.
"
function init(x::ABCDistance, sumstats::Array{Float64, 2})
    ##Default is to return original distance as no initialisation required
    x
end

function init(x::ABCDistance, sumstats::Array{Float64, 2}, parameters::Array{Float64, 2})
    ##Default is to ignore parameters, as these are often not needed
    init(x, sumstats)
end

type Lp <: ABCDistance
    sobs::Array{Float64,1}
    p::Float64
end

function evaldist(x::Lp, s::Array{Float64,1})
    absdiff = abs.(x.sobs - s)
    norm(absdiff, x.p)
end

"
Euclidean distance. Its constructor is `Euclidean(sobs)`.
"
type Euclidean <: ABCDistance
    sobs::Array{Float64,1}
end

function evaldist(x::Euclidean, s::Array{Float64,1})
    absdiff = abs.(x.sobs - s)
    norm(absdiff, 2.0)
end

type Logdist <: ABCDistance
    sobs::Array{Float64,1}
end

function evaldist(x::Logdist, s::Array{Float64,1})
    absdiff = abs.(x.sobs - s)
    exp(sum(log(absdiff)))
end

"""
This is weighted Euclidean distance with summary statistics weighted by scalars - equation (3) in the paper.
It has several constructors.
  * `WeightedEuclidean(sobs, w, scale_type)`
  Here `w` is the weights for each summary and `scale_type` is how to initialise the scale estimate. Possible values for the latter include `"sd"` and `"MAD"` (standard deviation and median absolute deviation).
  * `WeightedEuclidean(sobs, scale_type)`
  Here the weights are left undefined until initialisation of the distance takes place.
  * `WeightedEuclidean(sobs)`
  Equivalent to `WeightedEuclidean(sons, "MAD")`
"""
type WeightedEuclidean <: ABCDistance
    sobs::Array{Float64,1}
    w::Array{Float64,1} ##Weights for each summary statistic - square root of estimated precisions
    scale_type::AbstractString ##Whether to initialise scale using "MAD", "sd" or "ADO"

    function WeightedEuclidean(sobs::Array{Float64,1}, w::Array{Float64,1}, scale_type::AbstractString)
        if (scale_type ∉ ("MAD", "sd", "sdreg", "ADO"))
            error("scale_type must be MAD sd sdreg or ADO")
        end
        new(sobs, w, scale_type)
    end
end

function WeightedEuclidean(sobs::Array{Float64, 1}, scale_type::AbstractString)
    WeightedEuclidean(sobs, Array{Float64}(0), scale_type)
end

function WeightedEuclidean(sobs::Array{Float64, 1})
    WeightedEuclidean(sobs, Array{Float64}(0), "MAD")
end

function init(x::WeightedEuclidean, sumstats::Array{Float64, 2}, parameters::Array{Float64, 2})
    (nstats, nsims) = size(sumstats)
    if (nsims == 0)
        sig = ones(nstats)
    elseif x.scale_type=="MAD"
        sig = Float64[MAD(sumstats[i,:]) for i in 1:nstats]
    elseif x.scale_type=="sd"
        sig = Float64[std(sumstats[i,:]) for i in 1:nstats]
    elseif x.scale_type=="sdreg"
        P = parameters'
        sig = Array{Float64}(nstats)
        for i in 1:nstats
            y = sumstats[i,:]
            beta = linreg(P, y)
            res = y - beta[1] - P * beta[2:end]
            sig[i] = std(res)
        end
    elseif x.scale_type=="ADO"
        sig = [ADO(sumstats[i,:], x.sobs[i]) for i in 1:nstats]
    end
    return WeightedEuclidean(x.sobs, 1.0./sig, x.scale_type)
end

##Median absolute deviation
function MAD(x::Array{Float64, 1})
    median(abs.(x - median(x)))
end

##Absolute deviation to observations
function ADO(x::Array{Float64, 1}, obs::Float64)
    median(abs.(x - obs))
end

function evaldist(x::WeightedEuclidean, s::Array{Float64, 1})
    absdiff = abs.(x.sobs - s)
    wdiff = absdiff .* x.w
    wdiff[absdiff .== 0.0] = 0.0 # Correctly handles case where absdiff zero and w infinite
    norm(wdiff, 2.0)
end

"
This is weighted Euclidean distance with summary statistics weighted by a matrix i.e. d(x,y)=(x-y)' W (x-y). (This can often be viewed as an empirical estimate of Mahalanobis distance.)
It has two constructors.
  * `MahalanobisEmp(sobs, Ω)`
  Where `Ω` is used as the matrix W.
  * `MahalanobisEmp(sobs)`
  Here the weights are left undefined until initialisation of the distance takes place.
"
type MahalanobisEmp <: ABCDistance
    sobs::Array{Float64,1}
    Ω::Array{Float64,2} ##Weights matrix. Empirical estimate of precision used.
end

function MahalanobisEmp(sobs::Array{Float64, 1})
    MahalanobisEmp(sobs, eye(1))
end

function init(x::MahalanobisEmp, sumstats::Array{Float64, 2})
    (nstats, nsims) = size(sumstats)
    if (nsims == 0)
        Ω = eye(nstats)
    else
        Ω = inv(cov(sumstats, 2)) ##TO DO: need to deal with case where empirical covariance not invertible
    end
    return MahalanobisEmp(x.sobs, Ω)
end

function evaldist(x::MahalanobisEmp, s::Array{Float64, 1})
    absdiff = abs.(x.sobs - s)
    (absdiff' * x.Ω * absdiff)[1]
end


text:
"""
Package de fonctions qui définissent les fonctions Theta utilisées ici

liste des fonctions :
theta(x::Vector,r::Float64)

invtheta(x::Vector,r::Float64)

dtheta(x::Vector,r::Float64)

AlphaThetaMax(x::Float64,dx::Float64,y::Float64,dy::Float64,r::Float64,s::Float64,t::Float64)
AlphaThetaMaxOneLeaf(x::Float64,dx::Float64,y::Float64,dy::Float64,r::Float64,s::Float64,t::Float64)
"""

module Thetafunc

"""
La fonction theta
theta_r(x)
In : x :: vecteur ; r :: float
"""
function theta(x::Any,r::Float64)
 y=zeros(length(x))
 for i=1:length(x)
  y[i]= x[i]>=0 ? x[i]/(x[i]+r) : x[i]/r-x[i]^2/(r^2)
 end
 return y
end
theta(x::Float64,r::Float64)= x>=0 ? x/(x+r) : x/r-x^2/(r^2)

"""
La fonction inverse de theta
invtheta_r(x)
In : x :: vecteur ; r :: float
"""
function invtheta(x::Any,r::Float64)
 y=zeros(length(x))
 for i=1:length(x)
  y[i]= x[i]>=0 ? -r*x[i]/(x[i]-1) : 1/2*(r-r*sqrt(1-4*x[i]))
 end
 return y
end
invtheta(x::Float64,r::Float64)= x>=0 ? -r*x/(x-1) : 1/2*(r-r*sqrt(1-4*x))

"""
La fonction dérivée theta
dtheta_r(x)
In : x :: vecteur ; r :: float
"""
function dtheta(x::Any,r::Float64)
 y=zeros(length(x))
 for i=1:size(x,1)
  y[i]= x[i]>=0 ? r/(x[i]+r)^2 : 1/r-2*x[i]/(r^2)
 end
 return y
end
dtheta(x::Float64,r::Float64)= x>=0 ? r/(x+r)^2 : 1/r-2*x/(r^2)

"""
La fonction dérivée seconde theta
ddtheta_r(x)
In : x :: vecteur ; r :: float
"""
function ddtheta(x::Any,r::Float64)
 y=zeros(length(x))
 for i=1:size(x,1)
  y[i]= x[i]>=0 ? -2*r/(x[i]+r)^3 : -2/(r^2)
 end
 return y
end
ddtheta(x::Float64,r::Float64)= x>=0 ? -2*r/(x+r)^3 : -2/(r^2)

"""
Essaye de donner la plus petite racine réelle positive du polynome ax^2+bx+c
"""
function PolyDegre2Positif(a::Float64,b::Float64,c::Float64)
prec=1e-15
 if a!=0
 disc=b^2-4*a*c
  if disc>=0
   sol=min((-b+sqrt(disc))./(2*a),(-b-sqrt(disc))./(2*a))
   if sol<prec #normal si on touche déjà la frontière
    sol=max((-b+sqrt(disc))./(2*a),(-b-sqrt(disc))./(2*a))
   end
  else
   sol=Inf #pas de racines réelles
  end
 else #a==0
  sol=-c./b
 end

 return sol
end
"""
AlphaThetaMax : résoud l'équation en alpha où psi(x,r,s,t)=s+t*theta(x,r) (n=1)
In : x float (position en x)
     dx float (déplacement en x)
     y float (position en y)
     dy float (déplacement en y)
     r float (paramètre)
     s float (paramètre)
     t float (paramètre)
"""
function AlphaThetaMax(x::Float64,dx::Float64,y::Float64,dy::Float64,r::Float64,s::Float64,t::Float64)
 #donner la formule analytique pour Theta1 :
if dx==0.0 && dy==0.0
 return [Inf,Inf]
end

alphaf1=AlphaThetaMaxOneLeaf(x,dx,y,dy,r,s,t)

alphaf2=AlphaThetaMaxOneLeaf(y,dy,x,dx,r,s,t)

 #le premier fixe yg et le deuxième yh
 return [alphaf1,alphaf2]
end
"""
Résoud l'équation en a:
(x+a*dx-s)*(y+a*dy-s+r)-t*(y+a*dy-s)=0 -> alphaf1p
et si besoin résoud l'équation en a:
r^2[x+a*dx-s-t*((y+a*dy-s)/r-(y+a*dy-s)^2/r^2)]=0 -> alphaf1n
"""
function AlphaThetaMaxOneLeaf(x::Float64,dx::Float64,y::Float64,dy::Float64,r::Float64,s::Float64,t::Float64)

 if t!=0

 #a=dx.*dy erreur de calcul
 #b=dy.*(x-s+t)+dx.*(y-s+r)
 #c=x.*y-x*s+x*r-y*s+s^2-s*r+t*y-t*s
 a=dx.*dy
 b=dy.*(x-s-t)+dx.*(y-s+r)
 c=x.*y-x*s+x*r-y*s+s^2-s*r-t*y+t*s
 alphaf1p=PolyDegre2Positif(a,b,c)

  if alphaf1p>=0 && (x+alphaf1p*dx>=s && y+alphaf1p*dy>=s)
   alphaf1=alphaf1p
  elseif alphaf1p>=0 && (x+alphaf1p*dx<=s && y+alphaf1p*dy<=s)
   #contrainte candidate mais trop bas
   #a=-dy^2 erreur de calcul
   #b=r^2*dx+r*t*dy-2*y*dy+2*s*dy
   #c=r^2*x-r^2*s+r*t*y-r*t*s-y^2-s^2+2*s*y
   a=dy^2
   b=r^2*dx-r*t*dy+t*2*y*dy-t*2*s*dy
   c=r^2*x-r^2*s-r*t*y+r*t*s+t*y^2+t*s^2-t*2*s*y
   alphaf1n=PolyDegre2Positif(a,b,c)

   #alphaf1=alphaf1n devrait être actif si on tient compte de la contrainte "intérieur"
   alphaf1=Inf
  else
   #on ne tape pas cette contrainte
   alphaf1=Inf
  end
 else #cas où t=0
  if dx!=0 && ((s-x)./dx>0 || (s==x && dx>0))
   alphaf1=(s-x)./dx
  else
   alphaf1=Inf
  end
 end #fin de condition t!=0

 return alphaf1
end

#fin du module
end


text:
module LaborAPI

using HTTP
using JSON3

const BASE_URL = "https://apiv3.dol.gov/v3/get/"
"""
    get(endpoint::String, params::Dict{String, Any})

Call the DOL API at the provided endpoint and with the provided parameters.
The `params` dictionary must contain at a minimum your API key.

# Examples

```
using LaborAPI

params = Dict(
    "X-API-KEY" => "adslfj3987234kjaf987fdas",
    "fields" => "case_id, flsa_violtn_cnt"
)

endpoint = "whd/wh_enforcement"

res = LaborAPI.get(endpoint, params)
```
"""
function get(endpoint::String, params::Dict{String, Any})
    url = BASE_URL * endpoint
    try
        res = HTTP.get(
            url,
            ["Content-Type" => "application/json"],
            query=params,
            require_ssl_verification=false
        )
        if HTTP.status(res) == 200
            json_str = String(res.body)
            return JSON3.read(json_str)
        else
            return error("ERROR: Expected status code 200 but got status code: $(HTTP.status(res))")
        end
    catch e
        return e
    end
end

end


text:
# Math Functions

# TODO: vector types
const generic_types = [Cfloat,Cdouble]
const generic_types_float = [Cfloat]
const generic_types_double = [Cdouble]

# generically typed

for gentype in generic_types
@eval begin

@device_override Base.acos(x::$gentype) = @builtin_ccall("acos", $gentype, ($gentype,), x)
@device_override Base.acosh(x::$gentype) = @builtin_ccall("acosh", $gentype, ($gentype,), x)
@device_function acospi(x::$gentype) = @builtin_ccall("acospi", $gentype, ($gentype,), x)

@device_override Base.asin(x::$gentype) = @builtin_ccall("asin", $gentype, ($gentype,), x)
@device_override Base.asinh(x::$gentype) = @builtin_ccall("asinh", $gentype, ($gentype,), x)
@device_function asinpi(x::$gentype) = @builtin_ccall("asinpi", $gentype, ($gentype,), x)

@device_override Base.atan(y_over_x::$gentype) = @builtin_ccall("atan", $gentype, ($gentype,), y_over_x)
@device_override Base.atan(y::$gentype, x::$gentype) = @builtin_ccall("atan2", $gentype, ($gentype, $gentype), y, x)
@device_override Base.atanh(x::$gentype) = @builtin_ccall("atanh", $gentype, ($gentype,), x)
@device_function atanpi(x::$gentype) = @builtin_ccall("atanpi", $gentype, ($gentype,), x)
@device_function atanpi(y::$gentype, x::$gentype) = @builtin_ccall("atan2pi", $gentype, ($gentype, $gentype), y, x)

@device_override Base.cbrt(x::$gentype) = @builtin_ccall("cbrt", $gentype, ($gentype,), x)

@device_override Base.ceil(x::$gentype) = @builtin_ccall("ceil", $gentype, ($gentype,), x)

@device_override Base.copysign(x::$gentype, y::$gentype) = @builtin_ccall("copysign", $gentype, ($gentype, $gentype), x, y)

@device_override Base.cos(x::$gentype) = @builtin_ccall("cos", $gentype, ($gentype,), x)
@device_override Base.cosh(x::$gentype) = @builtin_ccall("cosh", $gentype, ($gentype,), x)
@device_function cospi(x::$gentype) = @builtin_ccall("cospi", $gentype, ($gentype,), x)

@device_override SpecialFunctions.erfc(x::$gentype) = @builtin_ccall("erfc", $gentype, ($gentype,), x)
@device_override SpecialFunctions.erf(x::$gentype) = @builtin_ccall("erf", $gentype, ($gentype,), x)

@device_override Base.exp(x::$gentype) = @builtin_ccall("exp", $gentype, ($gentype,), x)
@device_override Base.exp2(x::$gentype) = @builtin_ccall("exp2", $gentype, ($gentype,), x)
@device_override Base.exp10(x::$gentype) = @builtin_ccall("exp10", $gentype, ($gentype,), x)
@device_override Base.expm1(x::$gentype) = @builtin_ccall("expm1", $gentype, ($gentype,), x)

@device_override Base.abs(x::$gentype) = @builtin_ccall("fabs", $gentype, ($gentype,), x)

@device_function dim(x::$gentype, y::$gentype) = @builtin_ccall("fdim", $gentype, ($gentype, $gentype), x, y)

@device_override Base.floor(x::$gentype) = @builtin_ccall("floor", $gentype, ($gentype,), x)

@device_override Base.fma(a::$gentype, b::$gentype, c::$gentype) = @builtin_ccall("fma", $gentype, ($gentype, $gentype, $gentype), a, b, c)

@device_override Base.max(x::$gentype, y::$gentype) = @builtin_ccall("fmax", $gentype, ($gentype, $gentype), x, y)

@device_override Base.min(x::$gentype, y::$gentype) = @builtin_ccall("fmin", $gentype, ($gentype, $gentype), x, y)

# NOTE: Julia's mod behaves differently than fmod
#@device_override Base.mod(x::$gentype, y::$gentype) = @builtin_ccall("fmod", $gentype, ($gentype, $gentype), x, y)
# fract(x::$gentype, $gentype *iptr) = @builtin_ccall("fract", $gentype, ($gentype, $gentype *), x, iptr)

@device_override Base.hypot(x::$gentype, y::$gentype) = @builtin_ccall("hypot", $gentype, ($gentype, $gentype), x, y)

@device_override SpecialFunctions.lgamma(x::$gentype) = @builtin_ccall("lgamma", $gentype, ($gentype,), x)

@device_override Base.log(x::$gentype) = @builtin_ccall("log", $gentype, ($gentype,), x)
@device_override Base.log2(x::$gentype) = @builtin_ccall("log2", $gentype, ($gentype,), x)
@device_override Base.log10(x::$gentype) = @builtin_ccall("log10", $gentype, ($gentype,), x)
@device_override Base.log1p(x::$gentype) = @builtin_ccall("log1p", $gentype, ($gentype,), x)
@device_function logb(x::$gentype) = @builtin_ccall("logb", $gentype, ($gentype,), x)

@device_function mad(a::$gentype, b::$gentype, c::$gentype) = @builtin_ccall("mad", $gentype, ($gentype, $gentype, $gentype), a, b, c)

@device_function maxmag(x::$gentype, y::$gentype) = @builtin_ccall("maxmag", $gentype, ($gentype, $gentype), x, y)
@device_function minmag(x::$gentype, y::$gentype) = @builtin_ccall("minmag", $gentype, ($gentype, $gentype), x, y)

# modf(x::$gentype, $gentype *iptr) = @builtin_ccall("modf", $gentype, ($gentype, $gentype *), x, iptr)

@device_function nextafter(x::$gentype, y::$gentype) = @builtin_ccall("nextafter", $gentype, ($gentype, $gentype), x, y)

@device_override Base.:(^)(x::$gentype, y::$gentype) = @builtin_ccall("pow", $gentype, ($gentype, $gentype), x, y)
@device_function powr(x::$gentype, y::$gentype) = @builtin_ccall("powr", $gentype, ($gentype, $gentype), x, y)

@device_override Base.rem(x::$gentype, y::$gentype) = @builtin_ccall("remainder", $gentype, ($gentype, $gentype), x, y)

@device_function rint(x::$gentype) = @builtin_ccall("rint", $gentype, ($gentype,), x)

@device_override Base.round(x::$gentype) = @builtin_ccall("round", $gentype, ($gentype,), x)

@device_function rsqrt(x::$gentype) = @builtin_ccall("rsqrt", $gentype, ($gentype,), x)

@device_override Base.sin(x::$gentype) = @builtin_ccall("sin", $gentype, ($gentype,), x)
# sincos(x::$gentype, $gentype *cosval) = @builtin_ccall("sincos", $gentype, ($gentype, $gentype *), x, cosval)
@device_override Base.sinh(x::$gentype) = @builtin_ccall("sinh", $gentype, ($gentype,), x)
@device_function sinpi(x::$gentype) = @builtin_ccall("sinpi", $gentype, ($gentype,), x)

@device_override Base.sqrt(x::$gentype) = @builtin_ccall("sqrt", $gentype, ($gentype,), x)

@device_override Base.tan(x::$gentype) = @builtin_ccall("tan", $gentype, ($gentype,), x)
@device_override Base.tanh(x::$gentype) = @builtin_ccall("tanh", $gentype, ($gentype,), x)
@device_function tanpi(x::$gentype) = @builtin_ccall("tanpi", $gentype, ($gentype,), x)

@device_function tgamma(x::$gentype) = @builtin_ccall("tgamma", $gentype, ($gentype,), x)

@device_override Base.trunc(x::$gentype) = @builtin_ccall("trunc", $gentype, ($gentype,), x)

end
end

# generically typed -- only floats

for gentypef in generic_types_float

if gentypef !== Cfloat
@eval begin
@device_override Base.max(x::$gentypef, y::Cfloat) = @builtin_ccall("fmax", $gentypef, ($gentypef, Cfloat), x, y)
@device_override Base.min(x::$gentypef, y::Cfloat) = @builtin_ccall("fmin", $gentypef, ($gentypef, Cfloat), x, y)
end
end

end

# generically typed -- only doubles

for gentyped in generic_types_double

if gentyped !== Cdouble
@eval begin
@device_override Base.min(x::$gentyped, y::Cdouble) = @builtin_ccall("fmin", $gentyped, ($gentyped, Cdouble), x, y)
@device_override Base.max(x::$gentyped, y::Cdouble) = @builtin_ccall("fmax", $gentyped, ($gentyped, Cdouble), x, y)
end
end

end

# specifically typed

# frexp(x::Cfloat{n}, Cint{n} *exp) = @builtin_ccall("frexp", Cfloat{n}, (Cfloat{n}, Cint{n} *), x, exp)
# frexp(x::Cfloat, Cint *exp) = @builtin_ccall("frexp", Cfloat, (Cfloat, Cint *), x, exp)
# frexp(x::Cdouble{n}, Cint{n} *exp) = @builtin_ccall("frexp", Cdouble{n}, (Cdouble{n}, Cint{n} *), x, exp)
# frexp(x::Cdouble, Cint *exp) = @builtin_ccall("frexp", Cdouble, (Cdouble, Cint *), x, exp)

# ilogb(x::Cfloat{n}) = @builtin_ccall("ilogb", Cint{n}, (Cfloat{n},), x)
@device_function ilogb(x::Cfloat) = @builtin_ccall("ilogb", Cint, (Cfloat,), x)
# ilogb(x::Cdouble{n}) = @builtin_ccall("ilogb", Cint{n}, (Cdouble{n},), x)
@device_function ilogb(x::Cdouble) = @builtin_ccall("ilogb", Cint, (Cdouble,), x)

# ldexp(x::Cfloat{n}, k::Cint{n}) = @builtin_ccall("ldexp", Cfloat{n}, (Cfloat{n}, Cint{n}), x, k)
# ldexp(x::Cfloat{n}, k::Cint) = @builtin_ccall("ldexp", Cfloat{n}, (Cfloat{n}, Cint), x, k)
@device_override Base.ldexp(x::Cfloat, k::Cint) = @builtin_ccall("ldexp", Cfloat, (Cfloat, Cint), x, k)
# ldexp(x::Cdouble{n}, k::Cint{n}) = @builtin_ccall("ldexp", Cdouble{n}, (Cdouble{n}, Cint{n}), x, k)
# ldexp(x::Cdouble{n}, k::Cint) = @builtin_ccall("ldexp", Cdouble{n}, (Cdouble{n}, Cint), x, k)
@device_override Base.ldexp(x::Cdouble, k::Cint) = @builtin_ccall("ldexp", Cdouble, (Cdouble, Cint), x, k)

# lgamma_r(x::Cfloat{n}, Cint{n} *signp) = @builtin_ccall("lgamma_r", Cfloat{n}, (Cfloat{n}, Cint{n} *), x, signp)
# lgamma_r(x::Cfloat, Cint *signp) = @builtin_ccall("lgamma_r", Cfloat, (Cfloat, Cint *), x, signp)
# lgamma_r(x::Cdouble{n}, Cint{n} *signp) = @builtin_ccall("lgamma_r", Cdouble{n}, (Cdouble{n}, Cint{n} *), x, signp)
# Cdouble lgamma_r(x::Cdouble, Cint *signp) = @builtin_ccall("lgamma_r", Cdouble, (Cdouble, Cint *), x, signp)

# nan(nancode::uintn) = @builtin_ccall("nan", Cfloat{n}, (uintn,), nancode)
@device_function nan(nancode::Cuint) = @builtin_ccall("nan", Cfloat, (Cuint,), nancode)
# nan(nancode::Culong{n}) = @builtin_ccall("nan", Cdouble{n}, (Culong{n},), nancode)
@device_function nan(nancode::Culong) = @builtin_ccall("nan", Cdouble, (Culong,), nancode)

# pown(x::Cfloat{n}, y::Cint{n}) = @builtin_ccall("pown", Cfloat{n}, (Cfloat{n}, Cint{n}), x, y)
@device_override Base.:(^)(x::Cfloat, y::Cint) = @builtin_ccall("pown", Cfloat, (Cfloat, Cint), x, y)
# pown(x::Cdouble{n}, y::Cint{n}) = @builtin_ccall("pown", Cdouble{n}, (Cdouble{n}, Cint{n}), x, y)
@device_override Base.:(^)(x::Cdouble, y::Cint) = @builtin_ccall("pown", Cdouble, (Cdouble, Cint), x, y)

# remquo(x::Cfloat{n}, y::Cfloat{n}, Cint{n} *quo) = @builtin_ccall("remquo", Cfloat{n}, (Cfloat{n}, Cfloat{n}, Cint{n} *), x, y, quo)
# remquo(x::Cfloat, y::Cfloat, Cint *quo) = @builtin_ccall("remquo", Cfloat, (Cfloat, Cfloat, Cint *), x::Cfloat, y, quo)
# remquo(x::Cdouble{n}, y::Cdouble{n}, Cint{n} *quo) = @builtin_ccall("remquo", Cdouble{n}, (Cdouble{n}, Cdouble{n}, Cint{n} *), x, y, quo)
# remquo(x::Cdouble, y::Cdouble, Cint *quo) = @builtin_ccall("remquo", Cdouble, (Cdouble, Cdouble, Cint *), x, y, quo)

# rootn(x::Cfloat{n}, y::Cint{n}) = @builtin_ccall("rootn", Cfloat{n}, (Cfloat{n}, Cint{n}), x, y)
@device_function rootn(x::Cfloat, y::Cint) = @builtin_ccall("rootn", Cfloat, (Cfloat, Cint), x, y)
# rootn(x::Cdouble{n}, y::Cint{n}) = @builtin_ccall("rootn", Cdouble{n}, (Cdouble{n}, Cint{n}), x, y)
# rootn(x::Cdouble, y::Cint) = @builtin_ccall("rootn", Cdouble{n}, (Cdouble, Cint), x, y)

# TODO: half and native


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

function testdistribute(bnd1::Function, bnd2::Function, Nx::Int, Ny::Int, M::Int)::Float64
    dbase = distribute(bnd1, bnd2, (x,y)->0, Nx, Ny, M)
    sumL2 = 0.0
    for m in 1:M, n in 1:M
        sPatch = interpolatePatch(dbase[[m,n]], 2Nx, 2Ny).value
        fPatch = Float64[bnd1(x) + bnd2(y) for x in chebgrid(2Nx, M, m), y in chebgrid(2Ny, M, n)]
        sumL2 += (L2norm(fPatch, sPatch, chebweights(2Nx)/M, chebweights(2Ny)/M))^2
    end
    return sqrt(sumL2)
end

@test testdistribute(x->sin(pi*x), y->sin(pi*y), 20, 20, 2) < 1e-8


text:
## ----------------------------------------------------------------------------
# remote bk
const _REMOTE_BK_URL_CONFIG_KEY = "remote_url"

function _remote_url()::String
    config = read_configfile()
    url = get(config, _REMOTE_BK_URL_CONFIG_KEY, "")
    return url
end

_remote_backup_dir() = joinpath(assistant_dir(), "remote")

## ----------------------------------------------------------------------------
# bk registry

_bk_registry_dir(gl::GitLink) = joinpath(repo_dir(gl), "registry")

function _get_bk_registry(gl::GitLink)
    rdir = _bk_registry_dir(gl)
    isdir(rdir) ? readdir(rdir) : String[]
end

function _addto_bk_registry!(gl::GitLink, bk::String)
    rfile = joinpath(_bk_registry_dir(gl), basename(bk))
    mkpath(dirname(rfile))
    touch(rfile)
    return rfile
end

_clear_bk_registry!(gl::GitLink) = rm(_bk_registry_dir(gl); force = true, recursive = true)

function _update_bk_registry!(gl::GitLink, kb::String)

    lbk_dir = local_backups_dir(kb)
    !isdir(lbk_dir) && error("local backup dir not found: '$(lbk_dir)'")

    for file in readdir(lbk_dir)
        !endswith(file, ".zip") && continue
        _addto_bk_registry!(gl, file)
        _addto_bk_registry!(gl, _to_ascii(file))
    end

    _get_bk_registry(gl)
end

## ----------------------------------------------------------------------------
function _copy_tagged_release_action(repodir)
    srcfile = joinpath(pkgdir(RemAssistant), ".github", "workflows", "tagged-release.yml")
    destfile = joinpath(repodir, ".github", "workflows", "tagged-release.yml")
    mkpath(dirname(destfile))
    cp(srcfile, destfile; force = true)
    return destfile
end

## ----------------------------------------------------------------------------
# TODO: add 'run_git' at GitLinks
function _tag_backup!(gl::GitLink)
    lock(gl) do
        cd(repo_dir(gl)) do
            tagname = string("v", Dates.format(today(), "yyyy.m.d")) # this might comply semver regex
            try;
                p = run(`git tag $(tagname)`; wait = false)
                wait(p)
            catch ignored end
            try;
                p = run(`git push origin $(tagname)`; wait = false)
                wait(p)
            catch ignored end
        end
    end
end

## ----------------------------------------------------------------------------
_to_ascii(str) = replace(str, !isascii => "")

function _is_registered(newest, registered)
    findfirst(
        isequal(_to_ascii(basename(newest))),
        _to_ascii.(basename.(registered))
    ) !== nothing
end

## ----------------------------------------------------------------------------
function _update_remote_backup(bkdir, remote_url;
        force = false, loglevel = :verbose
    )

    ## gitlink
    gl = GitLink(bkdir, remote_url)

    ## copy/unzip local backup
    kb = _knowledgebase()
    newest = _resolve_newest_local_bk!(kb)

    # log
    if isempty(newest)
        loglevel == :verbose && @info("No local backup detected!!!")
        return nothing
    end

    # check registered
    registered = _get_bk_registry(gl)
    is_registered = _is_registered(newest, registered)
    if !force && is_registered
        loglevel == :verbose && @info("Nothing to update!!!", newest = basename(newest), force, is_registered)
        return nothing
    end

    ## update remote backup
    instantiate(gl; verbose = false)
    ok_flag = GitLinks.writewdir(gl) do wdir
        unzip(newest, wdir)
        _copy_tagged_release_action(wdir)
        _update_bk_registry!(gl, kb)
    end

    ## Add tag
    _tag_backup!(gl)

    ## log
    if (loglevel == :verbose || loglevel == :upload)
        ok_flag ?
            @info("Backup uploaded!!!", newest = basename(newest), force, is_registered) :
            @warn("Backup updating failed!!!", newest = basename(newest), force, is_registered)
    end
end


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

mutable struct MHState{
    Q<:AbstractProposalDist,
    S<:DensitySample
} <: AbstractMCMCState
    pdist::Q
    current_sample::S
    proposed_sample::S
    proposal_accepted::Bool
    current_nreject::Int64
    nsamples::Int64
    nsteps::Int64
end

function MHState(
    pdist::AbstractProposalDist,
    current_sample::DensitySample
)
    proposed_sample = DensitySample(
        similar(current_sample.params),
        convert(typeof(current_sample.log_value), NaN),
        zero(current_sample.weight)
    )

    MHState(
        pdist,
        current_sample,
        proposed_sample,
        false,
        zero(Int64),
        zero(Int64),
        one(Int64)
    )
end

nparams(state::MHState) = nparams(state.pdist)

nsteps(state::MHState) = state.nsteps

nsamples(state::MHState) = state.nsamples

eff_acceptance_ratio(state::MHState) = nsamples(state) / nsteps(state)

function next_cycle!(state::MHState)
    state.current_sample.weight = one(state.current_sample.weight)
    state.nsamples = zero(Int64)
    state.nsteps = zero(Int64)
    state
end

density_sample_type(state::MHState{Q,S}) where {Q,S} = S

_accrej_current_sample(state::MHState) =
    ifelse(state.proposal_accepted, state.current_sample, state.proposed_sample)

function nsamples_available(chain::MCMCIterator{<:MCMCAlgorithm{MHState}}, nonzero_weights::Bool = false)
    if nonzero_weights
        (_accrej_current_sample(chain.state).weight > 0) ? 1 : 0
    else
        1
    end
end

function get_samples!(appendable, chain::MCMCIterator{<:MCMCAlgorithm{MHState}}, nonzero_weights::Bool)::typeof(appendable)
    sample = _accrej_current_sample(chain.state)
    if (!nonzero_weights || sample.weight > 0)
        push!(appendable, sample)
    end
    appendable
end

function get_sample_ids!(appendable, chain::MCMCIterator{<:MCMCAlgorithm{MHState}}, nonzero_weights::Bool)::typeof(appendable)
    state = chain.state
    sample = _accrej_current_sample(state)
    if (!nonzero_weights || sample.weight > 0)
        sampletype = state.proposal_accepted ? 2 : 1
        sampleid = MCMCSampleID(chain.id, chain.cycle, state.nsteps, sampletype)
        push!(appendable, sampleid)
    end
    appendable
end

function MCMCIterator(
    algorithm::MCMCAlgorithm{MHState},
    likelihood::AbstractDensity,
    prior::AbstractDensity,
    id::Int64,
    rng::AbstractRNG,
    initial_params::AbstractVector{P} = Vector{P}(),
    exec_context::ExecContext = ExecContext(),
) where {P<:Real}
    target = likelihood * prior

    cycle = zero(Int)
    reset_rng_counters!(rng, id, cycle, 0)

    params_vec = Vector{P}(nparams(target))
    if isempty(initial_params)
        rand_initial_params!(rng, algorithm, prior, params_vec)
    else
        params_vec .= initial_params
    end

    !(params_vec in param_bounds(target)) && throw(ArgumentError("Parameter(s) out of bounds"))

    m = length(params_vec)

    log_value = density_logval(target, params_vec, exec_context)
    L = typeof(log_value)
    W = sample_weight_type(typeof(algorithm))

    current_sample = DensitySample(
        params_vec,
        log_value,
        one(W)
    )

    proposed_sample = DensitySample(
        similar(current_sample.params),
        convert(typeof(current_sample.log_value), NaN),
        zero(W)
    )

    state = MHState(
        algorithm,
        target,
        current_sample
    )

    chain = MCMCIterator(
        algorithm,
        target,
        state,
        rng,
        id,
        cycle,
        false,
        false
    )

    chain
end

exec_capabilities(mcmc_step!, callback::AbstractMCMCCallback, chain::MCMCIterator{<:MCMCAlgorithm{MHState}}) =
    exec_capabilities(density_logval, chain.target, chain.state.proposed_sample.params)

function mcmc_step!(
    callback::AbstractMCMCCallback,
    chain::MCMCIterator{<:MCMCAlgorithm{MHState}},
    exec_context::ExecContext,
    ll::LogLevel
)
    state = chain.state
    algorithm = chain.algorithm

    if !mcmc_compatible(algorithm, chain.state.pdist, param_bounds(chain.target))
        error("Implementation of algorithm $algorithm does not support current parameter bounds with current proposal distribution")
    end

    current_sample = state.current_sample
    proposed_sample = state.proposed_sample

    state.nsteps += 1
    reset_rng_counters!(chain)

    mcmc_propose_accept_reject!(callback, chain, exec_context)

    if state.proposal_accepted
        copy!(current_sample, proposed_sample)
        state.current_nreject = 0
        state.proposal_accepted = false
    end

    chain
end

abstract type MHWeightingScheme{T<:Real} end
export MHWeightingScheme

struct MHMultiplicityWeights{T<:Real} <: MHWeightingScheme{T} end
export MHMultiplicityWeights
MHMultiplicityWeights() = MHMultiplicityWeights{Int}()

struct MHAccRejProbWeights{T<:AbstractFloat} <: MHWeightingScheme{T} end
export MHAccRejProbWeights
MHAccRejProbWeights() = MHAccRejProbWeights{Float64}()

struct MHPosteriorFractionWeights{T<:AbstractFloat} <: MHWeightingScheme{T} end
export MHPosteriorFractionWeights
MHPosteriorFractionWeights() = MHPosteriorFractionWeights{Float64}()

struct MetropolisHastings{
    Q<:ProposalDistSpec,
    W<:Real,
    WS<:MHWeightingScheme{W}
} <: MCMCAlgorithm{MHState}
    q::Q
    weighting_scheme::WS
end

export MetropolisHastings

MetropolisHastings(
    q::Q,
    weighting_scheme::WS
) where {
    Q<:ProposalDistSpec,
    W<:Real,
    WS<:MHWeightingScheme{W}
} = MetropolisHastings{Q,W,WS}(q, weighting_scheme)

MetropolisHastings(q::ProposalDistSpec = MvTDistProposalSpec()) = MetropolisHastings(q, MHMultiplicityWeights())
MetropolisHastings(weighting_scheme::MHWeightingScheme) = MetropolisHastings(MvTDistProposalSpec(), weighting_scheme)

mcmc_compatible(::MetropolisHastings, ::AbstractProposalDist, ::NoParamBounds) = true

mcmc_compatible(::MetropolisHastings, pdist::AbstractProposalDist, bounds::HyperRectBounds) =
    issymmetric(pdist) || all(x -> x == hard_bounds, bounds.bt)

sample_weight_type(::Type{MetropolisHastings{Q,W,WS}}) where {Q,W,WS} = W

function MHState(
    algorithm::MetropolisHastings,
    target::AbstractDensity,
    current_sample::DensitySample{P,T,W}
) where {P,T,W}
    MHState(
        algorithm.q(P, nparams(current_sample)),
        current_sample
    )
end

function mcmc_propose_accept_reject!(
    callback::AbstractMCMCCallback,
    chain::MCMCIterator{<:MetropolisHastings},
    exec_context::ExecContext
)
    state = chain.state
    rng = chain.rng
    target = chain.target
    pdist = state.pdist

    current_sample = state.current_sample
    proposed_sample = state.proposed_sample

    current_params = current_sample.params
    proposed_params = proposed_sample.params

    current_log_value = current_sample.log_value
    T = typeof(current_log_value)

    # Propose new parameters:
    proposed_sample.weight = 0
    proposal_rand!(rng, pdist, proposed_params, current_params)
    apply_bounds!(proposed_params, param_bounds(target), false)

    p_accept = if proposed_params in param_bounds(target)
        # Evaluate target density at new parameters:
        proposed_log_value = T(density_logval(target, proposed_params, exec_context))

        # log of ratio of forward/reverse transition probability
        log_tpr = if issymmetric(pdist)
            T(0)
        else
            log_tp_fwd = distribution_logpdf(pdist, proposed_params, current_params)
            log_tp_rev = distribution_logpdf(pdist, current_params, proposed_params)
            T(log_tp_fwd - log_tp_rev)
        end

        current_log_value = current_sample.log_value
        proposed_sample.log_value = proposed_log_value

        p_accept = if proposed_log_value > -Inf
            clamp(T(exp(proposed_log_value - current_log_value - log_tpr)), zero(T), one(T))
        else
            zero(T)
        end
    else
        p_accept = zero(T)
        proposed_sample.log_value = -Inf

        zero(T)
    end

    mh_acc_rej!(chain, p_accept)

    callback(1, chain)

    chain
end

function mh_acc_rej!(
    chain::MCMCIterator{<:MetropolisHastings{Q,W,WS}},
    p_accept::Real
) where {Q, W, WS <: MHMultiplicityWeights}
    @assert p_accept >= 0
    state = chain.state
    state.current_sample.weight += 1
    if rand(chain.rng, float(typeof(p_accept))) < p_accept
        state.proposal_accepted = true
        state.nsamples += 1
    else
        state.current_nreject += 1
    end
    chain
end

function mh_acc_rej!(
    chain::MCMCIterator{<:MetropolisHastings{Q,W,WS}},
    p_accept::Real
) where {Q, W, WS <: MHAccRejProbWeights}
    @assert p_accept >= 0
    if p_accept ≈ 1
        p_accept = one(p_accept)
    elseif p_accept ≈ 0
        p_accept = zero(p_accept)
    end

    state = chain.state
    state.current_sample.weight += (1 - p_accept)
    state.proposed_sample.weight = p_accept
    if rand(chain.rng, float(typeof(p_accept))) < p_accept
        state.proposal_accepted = true
        state.nsamples += 1
    else
        state.current_nreject += 1
    end
    chain
end

function mh_acc_rej!(
    chain::MCMCIterator{<:MetropolisHastings{Q,W,WS}},
    p_accept::Real
) where {Q, W, WS <: MHPosteriorFractionWeights}
    @assert p_accept >= 0

    state = chain.state

    r = rand(chain.rng, float(typeof(p_accept)))

    if p_accept ≈ 0
        state.current_sample.weight += 1
        state.proposed_sample.weight = 0
    else
        # Renormalize posterior values:
        logval_1 = state.current_sample.log_value
        logval_2 = state.proposed_sample.log_value
        max_lv = max(logval_1, logval_2)
        v_1 = exp(logval_1 - max_lv)
        v_2 = exp(logval_2 - max_lv)
        v_sum = v_1 + v_2

        state.current_sample.weight += v_1 / v_sum
        state.proposed_sample.weight = v_2 / v_sum
        if r < p_accept
            state.proposal_accepted = true
            state.nsamples += 1
        end
    end
    chain
end


text:
function read_xyz(ifile::String)
    """
    Reads in an xyz file of possibly multiple geometries, returning the header, atom labels,
    and coordinates as arrays of strings and Float64s for the coordinates.
    """
    @time file_contents = readlines(ifile)
    header = Array{String, 1}()
    atom_labels = Array{Array{String, 1}, 1}()
    geoms = Array{Array{Float64, 2}, 1}()
    for (i, line) in enumerate(file_contents)
        if isa(tryparse(Int, line), Int)
            # allocate the geometry for this frame
            N = parse(Int, line)
            head = string(N)
            labels = String[]
            # store the header for this frame
            head = string(line, file_contents[i+1])
            i += 1
            push!(header, head)
            # loop through the geometry storing the vectors and atom labels as you go
            geom = zeros((3, N))
            for j = 1:N
                coords = split(file_contents[i+1])
                i += 1
                push!(labels, coords[1])
                geom[:,j] = parse.(Float64, coords[2:end])
            end
            push!(geoms, geom)
            push!(atom_labels, labels)
        end
    end
    return header, atom_labels, geoms
end

number, labels, geoms = read_xyz("C:\\Users\\Nutt\\Desktop\\github\\JuliaCC.jl\\test\\xyz\\benzene.xyz")
print(number)
print(labels)
println(geoms)


text:
type TopographicPoint2D <: Singleton
  Heatmap::Array{Float64,2}
  radiusOfInterest::Float64
end

function evalPotential(topo::TopographicPoint2D, Xi::Array{Graphs.ExVertex,1}; N::Int=100)
  warn("evalPotential TopographicPoint2D is not finished yet")
  return topo.Heatmap
end


text:
function loadleds(filename::AbstractString)
    (split.(split(line, " | ")) for line in eachline(filename))
end

function part1(data)
    sum(count(length(d) in (2, 3, 4, 7) for (_, ds) in data for d in ds))
end

using Test
if basename(pwd()) == "aoc"
    cd("2021\\8")
end
@test part1(loadleds("example.txt")) == 26
println("The answer to part 1 is: ", part1(loadleds("input.txt")))

part2(filename::AbstractString) = sum(decodeline(line) for line in eachline(filename))

function decodeline(line::AbstractString)
    signals, nums = split.(split(line, " | "))

    # Canonical configuration of signals is
    #    AAAA
    #   B    C
    #   B    C
    #    DDDD
    #   E    F
    #   E    F
    #    GGGG

    options = Dict(c => Set('a':'g') for c in 'A':'G')

    canon = Dict{Union{Set{Char},Int},Union{Set{Char},Int}}((
        Set("ABCEFG") => 0,
        Set("CF") => 1,
        Set("ACDEG") => 2,
        Set("ACDFG") => 3,
        Set("BCDF") => 4,
        Set("ABDFG") => 5,
        Set("ABDEFG") => 6,
        Set("ACF") => 7,
        Set("ABCDEFG") => 8,
        Set("ABCDFG") => 9
    ))
    push!(canon, (v => k for (k, v) in pairs(canon))...)

    deduce(canonchar::Char, chars) = intersect!(options[canonchar], chars)
    function deduce(canonchars, chars)
        foreach((cc) -> deduce(cc, chars), canonchars)
        options
    end

    function cleardeduced(options)
        for (cc, chars) in pairs(options)
            if length(chars) == 1
                for c in setdiff('A':'G', cc)
                    delete!(options[c], only(chars))
                end
            end
        end
        options
    end

    for signal in signals
        if length(signal) == 2 # v == 1
            deduce(canon[1], signal)
        elseif length(signal) == 3 # v == 7
            deduce(canon[7], signal)
        elseif length(signal) == 4 # v == 4
            deduce(canon[4], signal)
        elseif length(signal) == 5 # v in [2, 3, 5]
            deduce("ADG", signal)
        elseif length(signal) == 6 # v in [6, 9, 0]
            deduce("ABFG", signal)
        end # v == 8 tells us nothing
    end

    while any(!isone ∘ length, values(options))
        cleardeduced(options)
    end

    decoder = Dict(only(v) => k for (k, v) in pairs(options))

    mapfoldl(
        (n) -> canon[Set(replace(n, decoder...))],
        (a, b) -> 10a + b,
        nums,
        init=0
    )
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module SparseCntr

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOParametricInterval

nx = 3
A = spzeros(Float64,nx,nx)
c = [3.0, 0.5, 1.0]
P = spzeros(Float64,nx,nx)
A[1,1] = P[1,1] = 0.5
A[1,2] = P[1,2] = 0.2
A[1,3] = P[1,3] = 0.1
A[2,1] = P[2,1] = 0.3
A[2,2] = P[2,2] = 0.5
A[2,3] = P[2,3] = 0.5
A[3,1] = P[3,1] = 0.4
A[3,2] = P[3,2] = 0.5
A[3,3] = P[3,3] = 0.3

lufct = lufact(P)
At = copy(A)
ct = copy(c)
Pt = copy(P)
preA = inv(lufct)*A
prec = inv(lufct)*c
SSto = SparseInSto()
SSto.Xh = zeros(Float64,nx)
SSto.Yh = zeros(Float64,nx)
SSto.Zh = zeros(Float64,nx)
SSto.Xj = spzeros(Float64,nx,nx)
SSto.Yj = spzeros(Float64,nx,nx)
SSto.Zj = spzeros(Float64,nx,nx)
SSto.nx = nx
Sparse_Precondition!(c,A,P,SSto)

@test (7.7272-10E-4 <= c[1] <= 7.7272+10E-4)
@test (-5.0-10E-4 <= c[2] <= -5.0+10E-4)
@test (1.3636-10E-4 <= c[3] <= 1.3636+10E-4)

A1 = zeros(Float64,nx,nx)
c1 = [3.0, 0.5, 1.0]
P1 = zeros(Float64,nx,nx)
A1[1,1] = P1[1,1] = 0.5
A1[1,2] = P1[1,2] = 0.2
A1[1,3] = P1[1,3] = 0.1
A1[2,1] = P1[2,1] = 0.3
A1[2,2] = P1[2,2] = 0.5
A1[2,3] = P1[2,3] = 0.5
A1[3,1] = P1[3,1] = 0.4
A1[3,2] = P1[3,2] = 0.5
A1[3,3] = P1[3,3] = 0.3
#EAGOParametricInterval.Dense_Precondition!(c1,A1,P1,nx)
#@test (7.7272-10E-4 <= c1[1] <= 7.7272+10E-4)
#@test (-5.0-10E-4 <= c1[2] <= -5.0+10E-4)
#@test (1.3636-10E-4 <= c1[3] <= 1.3636+10E-4)

end


text:
push!(LOAD_PATH, homedir()*"/Projects/QuIPS/src")

# operator library

module Operators

using LinearAlgebra, Tensor

export Operator, Gate, ControlGate, Measurement
export Gates, addgate, addgates

#
#
# Constants
#
#

const GATES = Dict(:I => [1 0;
                          0 1],

                   :X => [0 1;
                          1 0], :NOT => [0 1;
                                         1 0],
                   :Y => [0 -im;
                          im 0],

                   :Z => [1 0;
                          0 -1],

                   :H => [1 1;
                          1 -1] / sqrt(2),

                   :S => [1 0;
                          0 im],

                   :T => [1 0;
                          0 exp(im * π / 4)])

const PARAM_GATES = Dict(:RX => γ -> exp(im * γ * GATES[:X]),

                         :RY => β -> exp(im * β * GATES[:Y]),

                         :RZ => α -> exp(im * α * GATES[:Z]),

                         :PHASE => θ -> [1 0; 0 exp(im * θ)])

#
# Gates: Set of all usable gates
#

struct Gates
    Us::Dict{Symbol,Matrix}
    PUs::Dict{Symbol, Function}

    Gates() = new(GATES, PARAM_GATES)
end

#### Functions to add additional gates to the set

# Add one gate
function addgate!(gates:: Gates, U::Tuple{Symbol, Matrix}):
    merge!(gates.Us, Dict(U))
end

function addgate!(gates:: Gates, PU::Tuple{Symbol, Function}):
    merge!(gates.PUs, Dict(PU))
end

# Add a list of gates
function addgates!(gates::Gates, Us::Vector{Tuple{Symbol, Matrix}}):
    for U in Us:
        adgate!(gates, U)
    end
end

function addgates!(gates::Gates, PUs::Vector{Tuple{Symbol, Function}}):
    for PU in PUs:
        adgate!(gates, PU)
    end
end

const COMP_BASIS = ([1, 0],
                    [0, 1])

const C32 = Complex{Float32}

#
#
# Operators
#
#

abstract type Operator end

#
# Measurement Operator
#

struct Measurement <: Operator
    k::Int
    basis::NTuple{2, Vector{C32}}

    Measurement(k::Int, β=COMP_BASIS) = new(k, β)
end

function (M::Measurement)(ψ::Vector{C32}, N::Int)
    k = M.k
    β₀ = M.basis[1]
    M₀ = C32.(β₀ * β₀')
    M̃₀ = tensor(M₀, k, N)
    P₀ = abs(ψ' * M̃₀ * ψ)
    u = rand(Float32)
    if u < P₀
        p = sqrt(P₀)
        return M̃₀/p * ψ, 0
    else
        β₁ = M.basis[2]
        M₁ = C32.(β₁ * β₁')
        M̃₁ = tensor(M₁, k, N)
        p = sqrt(1 - P₀)
        return M̃₁/p * ψ, 1
    end
end

#
# Single Qubit Gate
#

struct Gate <: Operator
    name::Symbol
    p::Union{Float32, Nothing} # parameter
    k::Int                     # target qubit
    U::Matrix{C32}

    Gate(tag::Tuple{Symbol,Int}) = begin
        name, k = tag;
        new(name, nothing, k, GATES[name])
    end

    Gate(tag::Tuple{Symbol,Float32,Int}) = begin
        name, p, k = tag;
        new(name, p, k, PARAM_GATES[name](p))
    end
end

(G::Gate)(ψ::Vector{C32}, N::Int) = tensor(G.U, G.k, N) * ψ

#
# ControlGate Operator
#

struct ControlGate <: Operator
    name::Symbol
    gate::Gate
    indx::Tuple{Vararg{Int}}

    function ControlGate(tag::Tuple{Symbol,Tuple{Vararg{Int}}})
        name, indx = tag
        U = Symbol(String(filter!(l -> l != 'C', [String(name)...])))
        k = indx[end]
        gate = Gate((U, k))
        new(name, gate, indx)
    end

    function ControlGate(tag::Tuple{Symbol,Float32,Tuple{Vararg{Int}}})
        name, p, indx = tag
        U = Symbol(String(filter!(l -> l != 'C', [String(name)...])))
        k = indx[end]
        gate = Gate((U, p, k))
        new(name, gate, indx)
    end
end

function (CG::ControlGate)(ψ::Vector{C32}, N::Int)
    I = GATES[:I]
    P̂₀ = COMP_BASIS[1] * COMP_BASIS[1]'
    P̂₁ = COMP_BASIS[2] * COMP_BASIS[2]'
    Û = CG.gate.U
    CÛ = P̂₀ ⊗ I + P̂₁ ⊗ Û
    Cs = length(CG.indx) - 1
    if Cs > 1
        for i = 2:Cs
            Iᵢ = diagm(ones(C32, 2^i))
            CU = P̂₀ ⊗ Iᵢ + P̂₁ ⊗ CÛ
        end
    end
    tensor(CÛ, CG.indx, N) * ψ
end

Î = Gate((:I, 1))

Base.:^(V̂::Operator, n::Int) = n < 1 ? [Î] : [V̂ for i = 1:n]

function Base.show(V::Operator)
    if typeof(V) == ControlGate
        println(V.name, " ", ["$j " for j in V.indx]...)
    elseif typeof(V) == Gate
        if V.name != :I
            println(V.name, " ", V.k)
        end
    else
        println("MEASURE ", V.k)
    end
end

end

max_stars_repo_name: aarontrowbridge/QuIPS

max_issues_repo_name: aarontrowbridge/QuIPS

max_forks_repo_name: aarontrowbridge/QuIPS


text:
# Run advection tests

push!(LOAD_PATH, abspath(joinpath(pwd(), "..")))

using PDESolver
#using Base.Test
using Base.Test
using ODLCommonTools
using PdePumiInterface  # common mesh interface - pumi
using SummationByParts  # SBP operators
using Utils
using SimpleODEMod
using NonlinearSolvers   # non-linear solvers
using ArrayViews
import ArrayViews.view
using Input

function clean_dict(collection)
  for i in keys(collection)
    delete!(collection, i)
  end
end

#------------------------------------------------------------------------------
# define tests and tags

using TestSystem
# define tags that will be used
include("../tags.jl")

# test list
global const SimpleODETests = TestList()

function test_eq4()
  @testset "---- testing SimpleODE ----" begin
    start_dir = pwd()
    cd("./eqn4/")
    fname = "input_vals_simpleODE.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    for i = 1:length(eqn.q_vec)
      @test isapprox( eqn.q_vec[i], 4.0) atol=1e-10
    end

    cd(start_dir)
  end  # end facts block

  return nothing
end

add_func1!(SimpleODETests, test_eq4, [TAG_SHORTTEST])

#------------------------------------------------------------------------------
# run tests
@testset "----- Running SimpleODE tests -----" begin
  runTestSystem(SimpleODETests, solvePDE, ARGS)
end

#------------------------------------------------------------------------------
# cleanup