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

# EMS simulation animation

# global variables; hacky...
const Client = HTTP.WebSockets.WebSocket{HTTP.ConnectionPool.Transaction{Sockets.TCPSocket}}
global animClients = [] # store open connections
global animSimQueue = Vector{Union{Simulation,String}}() # to store sims and sim filenames between animation request and start
global animPorts = Set{Int}() # localhost ports for animation, to be set

function decodeMessage(msg)
	return String(msg)
end

# parse message from html, extract values
function parseMessage(msg::String)
	msgSplit = readdlm(IOBuffer(msg))
	# msgType = msgSplit[1]
	# msgData = msgSplit[2:end]
	return msgSplit[1], msgSplit[2:end]
end

# create dictionary for sending messages to js
function createMessageDict(message::String)
	messageDict = Dict()
	messageDict["message"] = message
	return messageDict
	# return Dict([("message", message)])
end

# change message of messageDict
function changeMessageDict!(messageDict::Dict, message::String)
	messageDict["message"] = message
end

function writeClient!(client::Client, messageDict::Dict, message::String)
	# common enough lines to warrant the use of a function, I guess
	changeMessageDict!(messageDict, message)
	write(client, json(messageDict))
end

# set icons for ambulances, hospitals, etc.
function animSetIcons(client::Client)
	messageDict = createMessageDict("set_icons")
	pngFileUrl(filename) = string("data:image/png;base64,", filename |> read |> base64encode)
	iconPath = joinpath(@__DIR__, "..", "..", "assets", "animation", "icons")
	icons = JSON.parsefile(joinpath(iconPath, "icons.json"))
	# set iconUrl for each icon
	for (name, icon) in icons
		icon["options"]["iconUrl"] = pngFileUrl(joinpath(iconPath, string(name, ".png")))
	end
	merge!(messageDict, icons)
	write(client, json(messageDict))
end

# adds nodes from fGraph
function animAddNodes(client::Client, nodes::Vector{Node})
	messageDict = createMessageDict("add_node")
	for node in nodes
		messageDict["node"] = node
		write(client, json(messageDict))
	end
end

# adds arcs from rGraph, should only be called after animAddNodes()
function animAddArcs(client::Client, net::Network)
	messageDict = createMessageDict("add_arc")
	for arc in net.rGraph.arcs
		messageDict["arc"] = arc
		messageDict["node_indices"] = net.rArcFNodes[arc.index]
		write(client, json(messageDict))
	end
end

# calculate and set speeds for arcs in rGraph, should only be called after animAddArcs()
function animSetArcSpeeds(client::Client, map::Map, net::Network)
	# shorthand:
	fNodes = net.fGraph.nodes
	rNetTravels = net.rNetTravels

	messageDict = createMessageDict("set_arc_speed")
	for arc in net.rGraph.arcs
		# calculate arc distance
		dist = 0.0
		fNodeIndices = net.rArcFNodes[arc.index]
		for i = 1:length(fNodeIndices)-1
			dist += normDist(map, fNodes[fNodeIndices[i]].location, fNodes[fNodeIndices[i+1]].location)
		end

		# select minimum travel time
		arcTime = Inf
		for i = 1:length(rNetTravels)
			if arcTime > rNetTravels[i].arcTimes[arc.index]
				arcTime = rNetTravels[i].arcTimes[arc.index]
			end
		end

		messageDict["arc"] = arc
		messageDict["speed"] = dist / (arcTime*24) # km / hr
		write(client, json(messageDict))
	end
end

function animAddBuildings(client::Client, sim::Simulation)
	messageDict = createMessageDict("add_hospital")
	for h in sim.hospitals
		messageDict["hospital"] = h
		write(client, json(messageDict))
	end
	delete!(messageDict, "hospital")

	changeMessageDict!(messageDict, "add_station")
	for s in sim.stations
		messageDict["station"] = s
		write(client, json(messageDict))
	end
	# delete!(messageDict, "station")
end

function animAddAmbs!(client::Client, sim::Simulation)
	messageDict = createMessageDict("add_ambulance")
	for amb in sim.ambulances
		copy!(amb.currentLoc, getRouteCurrentLocation!(sim.net, amb.route, sim.time))
		copy!(amb.destLoc, amb.currentLoc)
		messageDict["ambulance"] = amb
		write(client, json(messageDict))
	end
end

# write frame updates to client
function updateFrame!(client::Client, sim::Simulation, time::Float)

	# check which ambulances have moved since last frame
	# need to do this before showing call locations
	messageDict = createMessageDict("move_ambulance")
	for amb in sim.ambulances
		ambLocation = getRouteCurrentLocation!(sim.net, amb.route, time)
		if ambLocation != amb.currentLoc
			copy!(amb.currentLoc, ambLocation)
			amb.movedLoc = true
			# move ambulance
			messageDict["ambulance"] = amb
			write(client, json(messageDict))
		else
			amb.movedLoc = false
		end
	end
	delete!(messageDict, "ambulance")

	changeMessageDict!(messageDict, "update_ambulance_dest")
	for amb in sim.ambulances
		if amb.status == ambMobilising
			copy!(amb.destLoc, sim.calls[amb.callIndex].location)
		else
			copy!(amb.destLoc, amb.route.endLoc)
		end
		messageDict["ambulance"] = amb
		write(client, json(messageDict))
	end
	delete!(messageDict, "ambulance")

	# determine which calls to remove, update, and add
	# need to do this after finding new ambulance locations
	# shorthand variable names:
	previousCalls = sim.previousCalls
	currentCalls = sim.currentCalls
	removeCalls = setdiff(previousCalls, currentCalls)
	updateCalls = intersect(previousCalls, currentCalls)
	addCalls = setdiff(currentCalls, previousCalls)
	changeMessageDict!(messageDict, "remove_call")
	for call in removeCalls
		copy!(call.currentLoc, Location())
		messageDict["call"] = call
		write(client, json(messageDict))
	end
	changeMessageDict!(messageDict, "move_call")
	for call in updateCalls
		updateCallLocation!(sim, call)
		messageDict["call"] = call
		write(client, json(messageDict))
	end
	changeMessageDict!(messageDict, "add_call")
	for call in addCalls
		copy!(call.currentLoc, call.location)
		call.movedLoc = false
		updateCallLocation!(sim, call)
		messageDict["call"] = call
		write(client, json(messageDict))
	end
	sim.previousCalls = copy(sim.currentCalls) # update previousCalls
end

# update call current location
function updateCallLocation!(sim::Simulation, call::Call)
	# consider moving call if the status indicates location other than call origin location
	if call.status == callGoingToHospital || call.status == callAtHospital
		amb = sim.ambulances[call.ambIndex]
		call.movedLoc = amb.movedLoc
		if amb.movedLoc
			copy!(call.currentLoc, amb.currentLoc)
		end
	end
end

function animateClient(client::Client)
	global animClients, animSimQueue

	push!(animClients, client)
	println("Client connected")

	# get first item in animSimQueue, or get filename now
	nextAnimItem = (length(animSimQueue) > 0 ? popfirst!(animSimQueue) : selectXmlFile())
	sim = nothing # init
	if typeof(nextAnimItem) == Simulation
		sim = nextAnimItem
	else
		@assert(typeof(nextAnimItem) == String)
		configFilename = nextAnimItem
		println("Initialising simulation from config:", configFilename)
		sim = initSim(configFilename; allowResim = true)
		println("...initialised")
	end

	# check if sim can be animated
	@assert(sim.initialised, "simulation has not been initialised; see initSim function")
	@assert(!sim.animating, "simulation is already being animated")
	@assert(!sim.writeOutput, "cannot animate simulation which is writing to output files")
	if !isdefined(sim, :backup)
		if !sim.used
			backup!(sim)
		else
			@warn("will not be able to restart animation because the simulation has been partially run and has no backup")
		end
	end
	if sim.complete
		reset!(sim)
	end

	# set map
	messageDict = createMessageDict("set_map_view")
	messageDict["map"] = sim.map
	write(client, json(messageDict))

	# set sim start time
	messageDict = createMessageDict("set_start_time")
	messageDict["time"] = sim.startTime
	write(client, json(messageDict))

	animSetIcons(client) # set icons before adding items to map
	animAddBuildings(client, sim)
	animAddAmbs!(client, sim)

	if sim.time > sim.startTime
		# set animation to current sim state
		@assert(!sim.complete) # not sure what would happen otherwise
		sim.previousCalls = Set{Call}() # in case animation has been re-opened from same sim state
		messageDict = createMessageDict("jump_to_time")
		messageDict["time"] = sim.time
		write(client, json(messageDict))
	end

	sim.animating = true

	messageDict = createMessageDict("")
	while !eof(client)
		msg = readavailable(client) # waits for message from client?
		msgString = decodeMessage(msg)
		(msgType, msgData) = parseMessage(msgString)

		if msgType == "prepare_next_frame"
			simTime = Float(msgData[1])
			simulateToTime!(sim, simTime)
			sim.time = simTime # otherwise sim.time only stores time of last event
			messageDict["time"] = simTime
			writeClient!(client, messageDict, "prepared_next_frame")

		elseif msgType == "get_next_frame"
			simTime = Float(msgData[1])
			updateFrame!(client, sim, simTime) # show updated amb locations, etc
			if !sim.complete
				messageDict["time"] = simTime
				writeClient!(client, messageDict, "got_next_frame")
			else
				# no events left, finish animation
				writeClient!(client, messageDict, "got_last_frame")
			end

		elseif msgType == "pause"

		elseif msgType == "stop"
			# reset
			reset!(sim)
			animAddAmbs!(client, sim)

		elseif msgType == "update_icons"
			try
				animSetIcons(client)
			catch e
				@warn("Could not update animation icons")
				@warn(e)
			end

		elseif msgType == "get_arcs"
			animAddNodes(client, sim.net.fGraph.nodes)
			animAddArcs(client, sim.net)
			animSetArcSpeeds(client, sim.map, sim.net)

		elseif msgType == "disconnect"
			sim.animating = false
			close(client)
			deleteat!(animClients, findfirst(isequal(client), animClients))
			println("Client disconnected")
			break
		else
			sim.animating = false
			error("Unrecognised message: ", msgString)
		end
	end
end

"""
	function animate!(sim::Union{Simulation,Nothing} = nothing;
		configFilename::String = "", port::Int = 8001, openWindow::Bool = true)
Open a web browser window to animate the simulation.
Will animate for either `sim` or `configFilename`. If neither of these are given then there will be a prompt for the simulation configuration filename once the browser window has opened.

# Keyword arguments
- `configFilename` is the name of the configuration file to load the simulation from; can be used instead of `sim`
- `port` is the port number for the local host url, e.g. `port = 8001` will use localhost:8001; this can only be set once for all animation windows
- `openWindow` can be set to `false` to prevent the window from being opened automatically, which is useful if you wish to use a non-default browser
"""
function animate!(sim::Union{Simulation,Nothing} = nothing;
	configFilename::String = "", port::Int = 8001, openWindow::Bool = true)
	@assert(sim == nothing || configFilename == "", "can only set one of: sim, configFilename")
	global animSimQueue
	if runAnimServer(port)
		if sim != nothing
			push!(animSimQueue, sim)
		elseif configFilename != ""
			push!(animSimQueue, configFilename)
		end
		openWindow ? openLocalhost(port) : println("waiting for window with port $port to be opened")
	end
end

function animate(; configFilename::String = "", port::Int = 8001, openWindow::Bool = true)
	animate!(configFilename = configFilename, port = port, openWindow = openWindow)
end

# creates and runs server for given port
# returns true if server is running, false otherwise
function runAnimServer(port::Int)
	@assert(port >= 0)

	# check if port already in use
	global animPorts
	if in(port, animPorts)
		return true # port already used for animation
	end
	try
		socket = Sockets.connect(port)
		if socket.status == Base.StatusOpen
			println("port $port is already in use, try another")
			return false
		end
	catch
	end

	# create and run server
	onepage = read("$sourceDir/animation/index.html", String)
	@async HTTP.listen(Sockets.localhost, port, readtimeout = 0) do http::HTTP.Stream
		if HTTP.WebSockets.is_upgrade(http.message)
			HTTP.WebSockets.upgrade(http) do client
				animateClient(client)
			end
		else
			h = HTTP.Handlers.RequestHandlerFunction() do req::HTTP.Request
				HTTP.Response(200, onepage)
			end
			HTTP.Handlers.handle(h, http)
		end
	end
	push!(animPorts, port)

	return true
end

# opens browser window for url
function openUrl(url::String)
	if Sys.iswindows()
		run(`$(ENV["COMSPEC"]) /c start $url`)
	elseif Sys.isapple()
		run(`open $url`)
	elseif Sys.islinux() || Sys.isbsd()
		run(`xdg-open $url`)
	end
end

# opens browser window for localhost:port
function openLocalhost(port::Int)
	openUrl("http://localhost:$(port)")
end

# JSON.lower for various types, to reduce length of string returned from json function
JSON.lower(n::Node) = Dict("index" => n.index, "location" => n.location)
JSON.lower(a::Arc) = Dict("index" => a.index)
JSON.lower(a::Ambulance) = Dict("index" => a.index, "currentLoc" => a.currentLoc, "destLoc" => a.destLoc)
JSON.lower(c::Call) = Dict("index" => c.index, "currentLoc" => c.currentLoc, "priority" => c.priority)
JSON.lower(h::Hospital) = Dict("index" => h.index, "location" => h.location)
JSON.lower(s::Station) = Dict("index" => s.index, "location" => s.location)


text:
# Autogenerated wrapper script for GMP_jll for aarch64-apple-darwin
export libgmp, libgmpxx

JLLWrappers.@generate_wrapper_header("GMP")
JLLWrappers.@declare_library_product(libgmp, "@rpath/libgmp.10.dylib")
JLLWrappers.@declare_library_product(libgmpxx, "@rpath/libgmpxx.4.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgmp,
        "lib/libgmp.10.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgmpxx,
        "lib/libgmpxx.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using DifferentialEquations
using ReactionMechanismSimulator

file_path = "mech/superminimal.rms"; # Change it to your personal path
#load mechanism dictionary
phaseDict = readinput(file_path);
#mechanism dictionaries index:  phaseDict[phasename]["Species" or "Reactions"]
spcs = phaseDict["phase"]["Species"];
rxns = phaseDict["phase"]["Reactions"];

#######################    INPUT     #######################
V0 = 8.314 * 1000 / 1e5;  # [m^3] # Comparison with chemkin
P0 = 1.0e+5; #[Pa]
# T0 = 1000; #[K]
H2 = 0.67;  #[mol]
O2 = 0.33;  #[mol]
N0 = H2 + O2;
############################################################

# Define the phase (how species thermodynamic and kinetic properties calculated)
ig = IdealGas(spcs,rxns,name="phase");
# Define initial condition of the reactor
initialconds = Dict(["V"=>V0,"P"=>P0, "H2"=>H2,"O2"=>O2]);
# Define the domain (encodes how system thermodynamic properties calculated)
domain,y0 = ConstantPDomain(phase=ig,initialconds=initialconds);

#######################    INPUT     #######################
t_final = 2.1; #[s]
solver = DifferentialEquations.CVODE_BDF();
abstol = 1e-20;
reltol = 1e-12;
############################################################

react = Reactor(domain,y0,(0.0,t_final)); #Create the reactor object

#solve the ode associated with the reactor
sol = solve(react.ode,solver, abstol=reltol,reltol=reltol);
bsol = Simulation(sol,domain);

sol(2.0000)[domain.indexes[end]]

max_stars_repo_name: jiweiqi/rmsadjoint

max_issues_repo_name: jiweiqi/rmsadjoint

max_forks_repo_name: jiweiqi/rmsadjoint


text:
"""
    pystr(o::PyObject)

Return a string representation of `o` corresponding to `str(o)` in Python.
"""
pystr(args...; kwargs...) =
    pystr(AbstractString, args...; kwargs...)
pystr(T::Type{<:AbstractString}, args...; kwargs...) =
    convert(T, pystr(PyObject, args...; kwargs...))
pystr(::Type{PyObject}, o) =
    PyObject(@pycheckn ccall(@pysym(:PyObject_Str), PyPtr, (PyPtr,), PyObject(o)))
pystr(::Type{PyObject}, o::AbstractString="") =
    convertpystr(o)

"""
    pyrepr(o::PyObject)

Return a string representation of `o` corresponding to `repr(o)` in Python.
"""
pyrepr(o) =
    pystr(AbstractString, o)
pyrepr(::Type{T}, o) where {T<:AbstractString} =
    convert(T, pyrepr(PyObject, o))
pyrepr(::Type{PyObject}, o) =
    PyObject(@pycheckn ccall(@pysym(:PyObject_Repr), PyPtr, (PyPtr,), PyObject(o)))

"""
    pyisinstance(o::PyObject, t::PyObject)

True if `o` is an instance of `t`, corresponding to `isinstance(o,t)` in Python.
"""
pyisinstance(o, t) =
    pyisinstance(Bool, o, t)
pyisinstance(::Type{Bool}, o, t) =
    pyisinstance(Bool, PyObject(o), PyObject(t))
pyisinstance(::Type{Bool}, o::PyObject, t::PyObject) =
    !ispynull(t) && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t) == 1
pyisinstance(::Type{Bool}, o::PyObject, t::Union{Ptr{Cvoid},PyPtr}) =
    t != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t) == 1
pyisinstance(::Type{PyObject}, o, t) =
    convertpybool(pyisinstance(Bool, o, t))

"""
    pyistrue(o::PyObject)

True if `o` is considered to be true, corresponding to `not not o` in Python.
"""
pyistrue(o::PyObject) =
    pyistrue(Bool, o)
pyistrue(::Type{Bool}, o::PyObject) =
    (@pycheckz ccall(@pysym(:PyObject_IsTrue), Cint, (PyPtr,), o)) == 1
pyistrue(::Type{PyObject}, o) =
    convertpybool(pyistrue(Bool, o))

"""
    pynot(o::PyObject)

True if `o` is not considered to be true, corresponding to `not o` in Python.
"""
pynot(o::PyObject) =
    pynot(Bool, o)
pynot(::Type{Bool}, o::PyObject) =
    (@pycheckz ccall(@pysym(:PyObject_Not), Cint, (PyPtr,), o)) == 1
pynot(::Type{PyObject}, o) =
    convertpybool(pynot(Bool, o))

"""
    pyint(o::PyObject)
"""
pyint(args...; kwargs...) =
    pyint(PyObject, args...; kwargs...)
pyint(T::Type, args...; kwargs...) =
    pycall(pybuiltin("int"), T, args...; kwargs...)
pyint(::Type{PyObject}, o::Integer=0) =
    convertpyint(o)

"""
    pybool(o::PyObject)
"""
pybool(args...; kwargs...) =
    pybool(PyObject, args...; kwargs...)
pybool(T::Type, args...; kwargs...) =
    pycall(pybuiltin("bool"), T, args...; kwargs...)
pybool(::Type{PyObject}, o::Bool=false) =
    convertpybool(o)

"""
    pystring(o::PyObject)

Return a string representation of `o`.  Normally, this corresponds to `repr(o)`
in Python, but unlike `repr` it should never fail (falling back to `str` or to
printing the raw pointer if necessary).
"""
function pystring(o::PyObject)
    if ispynull(o)
        return "NULL"
    else
        s = ccall((@pysym :PyObject_Repr), PyPtr, (PyPtr,), o)
        if (s == C_NULL)
            pyerr_clear()
            s = ccall((@pysym :PyObject_Str), PyPtr, (PyPtr,), o)
            if (s == C_NULL)
                pyerr_clear()
                return string(PyPtr(o))
            end
        end
        return convert(AbstractString, PyObject(s))
    end
end


text:
println("\nTesting variables\n")

n, m, k = 3, 4, 5
x1t = Variable(Float32, n)
x1 = Variable(n)
x1i = Variable(randn(n))
x2 = Variable(n, m)
x2t = Variable(Float32, n, m)
x2i = Variable(randn(n, m))
x3 = Variable(n, m, k)
x3t = Variable(Float32, n, m, k)
xx = randn(n,m,k)
x3i = Variable(xx)

@test eltype(x3i) == eltype(xx)
@test size(x3i) == size(xx)
@test size(x3i,1) == size(xx,1)
@test xx == (~x3i)

@test typeof(operator(x1)) <: Eye
@test variables(x1) == x1


text:
using LazySets

function compare_to_dreal(state_vars::Array{Symbol}, control_vars::Array{Symbol}, input_set::Dict, controller_file::String, dynamics_map::Dict, dt::T where T<:Real, output_constraints::Array{Expr}, dirname::String, experiment_name::String, N_steps::T where T<:Real; jobs=1, dreal_path="/opt/dreal/4.21.06.1/bin/dreal")
    """
    # Basic overview of what this function does:
    # assert init state variables in init set
    # for t = 0:N-1
        # u_var = add controller
        # state_vars = add dynamics
        # check property at t+1
    # end
    """
    # initializations:
    # SMTLibFormula container to store smt2
    formula = SMTLibFormula()
    # load controller as network object
    controller = read_nnet(controller_file)
    # turn controller into u expression(s) just like symbolic dynamics (maybe just RHS?)
    u_expr = network2expr(controller, state_vars)

    # assert init state variables in init set
    # timestamp state vars in init set to be e.g. x_0
    timed_input_set = Dict((Meta.parse(string(k)*"_0"),v) for (k,v) in input_set)
    push!(formula.formula, define_domain(timed_input_set, formula.stats)...)

    results = []
    write_result(dirname*experiment_name, experiment_name*"\n"; specifier="w")
    tstart = time()
    t_dreal_total = 0
    for t = 0:N_steps-1
        # actual important lines of code:
        u_t = add_controller(control_vars, u_expr::Array{Expr}, state_vars, formula::SMTLibFormula, t)
        x_tp1 = add_dynamics(state_vars, control_vars, t, dynamics_map, dt, formula)
        result, t_dreal = add_output_constraints_and_check_property(formula, output_constraints, state_vars, t+1; jobs=jobs, dreal_path=dreal_path, experiment_name=experiment_name, dirname=dirname)
        # time stuff and results recording:
        t_sofar = time() - tstart
        t_dreal_total += t_dreal
        msg = "Property holds for timestep $(t+1) ? $result . Elapsed time: $(t_sofar) sec. Total time spent IN dreal: $(t_dreal_total) \n"
        print(msg)
        write_result(dirname*experiment_name, msg; specifier="a")
        push!(results, result)
        if t_sofar > 12*60*60
            write_result(dirname*experiment_name, "Breaking after $(t_sofar/60/60) hours. TIMEOUT. \n"; specifier="a")
            break;
        end
    end
    ΔT = time() - tstart
    msg = "Property holds for all timesteps? $(all(results)). Total elpased time: $ΔT sec\n"
    println(msg)
    write_result(dirname*experiment_name, msg; specifier="a")
    return ΔT
end


text:
using Mimi

include("region_parameters.jl")
include("gross_economy.jl")
include("emissions.jl")

@Mimi.defmodel model begin

    index[time] = 2015:5:2110

    # Note that the regions of your model must be specified here
    index[regions] = [:Region1, :Region2, :Region3]

    # Order matters here. If the emissions component were defined first, the model would not run.
    component(grosseconomy)
    component(emissions)

    # Set parameters for the grosseconomy component
    grosseconomy.l     = l
    grosseconomy.tfp   = tfp
    grosseconomy.s     = s
    grosseconomy.depk  = depk
    grosseconomy.k0    = k0
    grosseconomy.share = 0.3

    # Set parameters for the emissions component
    emissions.sigma = sigma

    # Connect parameters
    grosseconomy.YGROSS => emissions.YGROSS
end


text:
# Show how to do joins. First joins on two distinct columns, then joins on a range of columns.
using PyPlot

# Set max elements in iterator
MaxElem = 1000

# Pick two columns to join
col1 = "In|1,"
col2 = "In|100,"

# Find all rows with these columns
E1 = Tedge[getrow(Tedge[:,col1]),:]
E2 = Tedge[getrow(Tedge[:,col2]),:]
Ejoin = E2[:,getcol(E1)]

# Display
figure()
spy(Ejoin)
axis("auto")
xlabel("vertex")
ylabel("edge")

# Now we are finding all the out vertices with in vertices of both 111 and 222

# Set column ranges
colRange1 = StartsWith("In|111,")
colRange2 = StartsWith("In|222,")

# Set up query iterators
TedgeIt1 = getiterator(Tedge,MaxElem)
TedgeIt2 = getiterator(Tedge,MaxElem)

# Start first query iterator
E1 = logical(TedgeIt1[:,colRange1])
E1outDeg = Assoc("","","")

while nnz(E1) > 0

    global E1, E1outDeg

    # Combine
    E1 = logical(Tedge[getrow(E1),:]);
    E1outDeg = E1outDeg + sum(E1[:,StartsWith("Out|")],1)
    # Run next query iterator
    E1 = logical(TedgeIt1[])
end

# Start second query iterator
E2 = logical(TedgeIt2[:,colRange2])
E2outDeg = Assoc("","","")

while nnz(E2) > 0

    global E2, E2outDeg

    # Combine
    E2 = logical(Tedge[getrow(E2),:]);
    E2outDeg = E2outDeg + sum(E2[:,StartsWith("Out|")],1)
    # Run next query iterator
    E2 = logical(TedgeIt2[])
end

# Add together
E12outDeg = E1outDeg + E2outDeg
# Join columns
EjoinRange = E12outDeg[:,getcol(E2outDeg[:,getcol(E1outDeg)])]


text:
using SynapseClient
using DISSEQT.AlignUtils
using DISSEQT.SynapseTools

# Name of the run. (Should corrsespond to a subfolder of bamPath if in Synapse.)
runName = "H03UAAFXX"

# set to true to upload files
doUpload = true

# Should point to MyProject/Analysis/Alignment
alignmentFolder = "syn18694207" # VignuzziLabPublic/Projects/FitnessLandscapes/Analysis/Alignment

syn = SynapseClient.login()
coverageplots(syn, alignmentFolder, runName, doUpload, @__FILE__)


text:
"""
    s = buffer(input; buf_size = Inf, timespan = 1, type_stable = false)

Create a signal who buffers updates to signal `input` until maximum size of `buf_size`
or until `timespan` seconds have passed. The signal value is the last full buffer emitted
or an empty vector if the buffer have never been filled before.

Buffer type will be `Any` unless `type_stable` is set to `true`, then it will be set
to the value of the first encountered item.
"""
function buffer(input; buf_size = Inf, timespan = 1, type_stable = false)
    _buf = type_stable ? typeof(input())[] : Any[]
    sbuf = foldp(push!, _buf, input)
    cond = Signal(sbuf; state = time()) do in, state
        last_update = state.x
        (time() - last_update > timespan) || (length(in) >= buf_size)
    end
    when(cond, sbuf; v0 = _buf) do buf
        state(cond, time())
        frame = copy(buf)
        empty!(_buf)
        frame
    end

end
export buffer

"""
    debounce(f, args...; delay = 1, v0 = nothing)

Create a `Signal` whos action `f(args...)` will be called only after `delay` seconds
have passed since the last time its `args` were updated. Only works in push based paradigm.
If `v0` is not specified then the initial value is `f(args...)`.
"""
function debounce(f, args...; delay = 1, v0 = nothing)
    timer = Timer(identity, 0;interval = 0)
    f_args = PullAction(f, args)
    debounced_signal = Signal(v0 === nothing ? f_args() : v0)
    Signal(args...) do args
        finalize(timer)
        timer = Timer(t -> debounced_signal(f_args()), delay)
    end
    debounced_signal
end
export debounce

abstract type Throttle <: PullType end
"""
    throttle(f::Function, args...; maxfps = 0.03)

Create a throttled `Signal` whos action `f(args...)` will be called only
if `1/maxfps` time has passed since the last time it updated. The resulting `Signal`
will be updated maximum of `maxfps` times per second.
"""
function throttle(f::Function, args...; maxfps = 30)
    pa = PullAction(f, args, Throttle)
    sd = SignalData(pa())
    state = Ref((1/maxfps, time()))
    Signal(sd, pa, state)
end
export throttle

function (pa::PullAction{Throttle, A})(s) where A
    (dt, last_update) = s.state.x
    args = pull_args(pa)
    if !valid(s)
        if time() - last_update < dt
            validate(s)
        else
            s.state.x = (dt, time())
            pa()
        end
    end
    value(s)
end

function activate_timer(s, dt, duration)
    signalref = WeakRef(Ref(s))
    start_time = time()
    t = Timer(dt;interval = dt) do t
        time_passed = time() - start_time
        if time_passed > duration || signalref.value === nothing
            finalize(t)
        else
            push!(signalref.value.x, time(), true)
        end
        nothing
    end
    t
end

"""
    s = every(dt; duration = Inf)

A signal that updates every `dt` seconds to the current timestamp, for `duration` seconds.
"""
function every(dt; duration = Inf)
    res = Signal(time())
    activate_timer(res, dt, duration)
    res
end
export every

"""
    s = fps(freq; duration = Inf)

A signal that updates `freq` times a second to the current timestamp, for `duration` seconds.
"""
function fps(freq; duration = Inf)
    every(1/freq; duration = duration)
end
export fps

"""
    s = fpswhen(switch::Signal, freq; duration = Inf)

A signal that updates 'freq' times a second to the current timestamp, for `duration` seconds
if and only if the value of `switch` is `true`.
"""
function fpswhen(switch::Signal, freq; duration = Inf)
    res = Signal(time())
    timer = Timer(0 ; interval = 0)
    Signal(droprepeats(switch)) do sw
        if sw === true
            timer = activate_timer(res, 1/freq, duration)
        else
            finalize(timer)
        end
    end
    res
end
export fpswhen

"""
     s = for_signal(f::Function, range, args...; fps = 1)

Create a `Signal` that updates to `f(i,args....) for i in range` every `1/fps` seconds.
`range` and `args` can be of type `Signal` or any other type. The loop starts whenever
one of the arguments updates. If the previous for loop did not
complete it gets cancelled.

# Example

    rng = Signal(1:5)
    A = Signal(2)
    B = for_signal(rng,A;fps = 30) do i,a
        println(a^i)
    end;
"""
function for_signal(f::Function, rng, args...; fps = 1)
    i_sig = Signal(first(pull!(rng)))
    timer = Timer(0; interval = 0)
    Signal(args...; v0 = timer) do args
        current_range = pull!(rng)
        finalize(timer)
        iteration = iterate(current_range)
        (iteration === nothing) && return
        (item, state) = iteration
        i_sig(item)
        timer = Timer(1/fps; interval = 1/fps) do t
            iteration = iterate(current_range,state)
            if iteration === nothing
                finalize(t)
            else
                (item, state) = iteration
                i_sig(item)
            end
            nothing
        end
    end
    Signal(i_sig; v0 = nothing) do i
        f(i, pull_args(args)...)
    end
end
export for_signal

nothing


text:
module FdeSolver

using SpecialFunctions
using LinearAlgebra
using FFTW

include("main.jl")
include("main_Jacob.jl")
include("SupFuns.jl")
include("SupFuns_Jacob.jl")
include("dummy_fun.jl")
include("structs.jl")

export(FDEsolver)

end


text:
# Compatibility with older Julia versions

# TODO: Remove when we drop support for v1.0
@static if VERSION < v"1.2"
    hasproperty(x, p) = p in propertynames(x)
end


text:
ECHART_FUNCTION_ATTRS=["formatter"]

function echart_json(v,f_mode)
    if f_mode
        return v
    else
        return JSON.json(v)
    end
end

echart_json(a::Array,f_mode)="[$(join(echart_json.(a,f_mode),","))]"

function echart_json(d::Dict,f_mode::Bool=false)
    els=[]
    for (k,v) in d
        if k in ECHART_FUNCTION_ATTRS
            j="\"$k\": $(echart_json(v,true))"
        else
            j="\"$k\":"*echart_json(v,f_mode==false ? false : true)
        end
        push!(els,j)
    end
    return "{$(join(els,","))}"
end


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d
include("../LF_common/sparse_LF3d.jl")
include("build_nabla2_z.jl")

const ⊗ = kron

function test_main( NN::Array{Int64} )

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  @printf("Using loop")
  @time Lz_v1 = build_nabla2_z(LF)
  #println(Lz_v1)

  #@printf("Using Kronecker product directly")
  #@time Lz_v2 = speye(Nx) ⊗ speye(Ny) ⊗ LF.LFz.D2jl

  # Check whether there are any differences between the two matrices
  #println(Lz_v1 - Lz_v2)

end

test_main([80,80,80])


text:
# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# Description

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function simplify_kmer_graph(kmer_graph)
#     @info "simplifying kmer graph"
#     @info "resolving untigs..."
#     @time untigs = resolve_untigs(kmer_graph)
#     @info "determining untig orientations..."
#     oriented_untigs = determine_oriented_untigs(kmer_graph, untigs)
#     simplified_graph = MetaGraphs.MetaDiGraph(length(oriented_untigs))
#     MetaGraphs.set_prop!(simplified_graph, :k, kmer_graph.gprops[:k])
#     @info "initializing graph node metadata"
#     for (vertex, untig) in enumerate(oriented_untigs)
#         MetaGraphs.set_prop!(simplified_graph, vertex, :sequence, untig.sequence)
#         MetaGraphs.set_prop!(simplified_graph, vertex, :path, untig.path)
#         MetaGraphs.set_prop!(simplified_graph, vertex, :orientations, untig.orientations)
#         MetaGraphs.set_prop!(simplified_graph, vertex, :evidence, untig.evidence)
#     end

#     # determine oriented edges of simplified graph
#     simplified_untigs = []
#     @info "creating simplified untigs to resolve connections"
#     for vertex in Graphs.vertices(simplified_graph)
#         in_kmer = simplified_graph.vprops[vertex][:path][1] => simplified_graph.vprops[vertex][:orientations][1]
#         out_kmer = simplified_graph.vprops[vertex][:path][end] => simplified_graph.vprops[vertex][:orientations][end]
#     #     @show vertex, in_kmer, out_kmer
#         push!(simplified_untigs, in_kmer => out_kmer)
#     end

#     @info "resolving connections"
#     ProgressMeter.@showprogress for (ui, u) in enumerate(simplified_untigs)
#         for (vi, v) in enumerate(simplified_untigs)
#     #         + => +
#             source_kmer_index, source_orientation = last(u)
#             destination_kmer_index, destination_orientation = first(v)
#             edge = Graphs.Edge(source_kmer_index, destination_kmer_index)
#             if Graphs.has_edge(kmer_graph, edge)
#                 source_orientation_matches = (kmer_graph.eprops[edge][:orientations].source_orientation == source_orientation)
#                 destination_orientation_matches = (kmer_graph.eprops[edge][:orientations].destination_orientation == destination_orientation)
#                 if source_orientation_matches && destination_orientation_matches
# #                     @show "right orientation!! + +"

#                     simplified_graph_edge = Graphs.Edge(ui, vi)

#                     Graphs.add_edge!(simplified_graph, simplified_graph_edge)
#                     edge_orientations = (
#                         source_orientation = source_orientation,
#                         destination_orientation = destination_orientation
#                     )
#                     set_metadata!(simplified_graph, simplified_graph_edge, :orientations, edge_orientations)
#                 end
#             end
#     #         + => -
#             source_kmer_index, source_orientation = last(u)
#             destination_kmer_index, destination_orientation = last(v)
#             destination_orientation = !destination_orientation

#             edge = Graphs.Edge(source_kmer_index, destination_kmer_index)
#             if Graphs.has_edge(kmer_graph, edge)
#                 source_orientation_matches = (kmer_graph.eprops[edge][:orientations].source_orientation == source_orientation)
#                 destination_orientation_matches = (kmer_graph.eprops[edge][:orientations].destination_orientation == destination_orientation)
#                 if source_orientation_matches && destination_orientation_matches
# #                     @show "right orientation!! + -"
#                     simplified_graph_edge = Graphs.Edge(ui, vi)

#                     Graphs.add_edge!(simplified_graph, simplified_graph_edge)
#                     edge_orientations = (
#                         source_orientation = source_orientation,
#                         destination_orientation = destination_orientation
#                     )
#                     set_metadata!(simplified_graph, simplified_graph_edge, :orientations, edge_orientations)
#                 end
#             end
#     #         - => +
#             source_kmer_index, source_orientation = first(u)
#             source_orientation = !source_orientation
#             destination_kmer_index, destination_orientation = first(v)

#             edge = Graphs.Edge(source_kmer_index, destination_kmer_index)
#             if Graphs.has_edge(kmer_graph, edge)
#                 source_orientation_matches = (kmer_graph.eprops[edge][:orientations].source_orientation == source_orientation)
#                 destination_orientation_matches = (kmer_graph.eprops[edge][:orientations].destination_orientation == destination_orientation)
#                 if source_orientation_matches && destination_orientation_matches
# #                     @show "right orientation!! - +"
#                     simplified_graph_edge = Graphs.Edge(ui, vi)

#                     Graphs.add_edge!(simplified_graph, simplified_graph_edge)
#                     edge_orientations = (
#                         source_orientation = source_orientation,
#                         destination_orientation = destination_orientation
#                     )
#                     set_metadata!(simplified_graph, simplified_graph_edge, :orientations, edge_orientations)
#                 end
#             end
#     #         - => -
#             source_kmer_index, source_orientation = first(u)
#             source_orientation = !source_orientation
#             destination_kmer_index, destination_orientation = last(v)
#             destination_orientation = !destination_orientation

#             edge = Graphs.Edge(source_kmer_index, destination_kmer_index)
#             if Graphs.has_edge(kmer_graph, edge)
#                 source_orientation_matches = (kmer_graph.eprops[edge][:orientations].source_orientation == source_orientation)
#                 destination_orientation_matches = (kmer_graph.eprops[edge][:orientations].destination_orientation == destination_orientation)
#                 if source_orientation_matches && destination_orientation_matches
# #                     @show "right orientation!! - -"
#                     simplified_graph_edge = Graphs.Edge(ui, vi)

#                     Graphs.add_edge!(simplified_graph, simplified_graph_edge)
#                     edge_orientations = (
#                         source_orientation = source_orientation,
#                         destination_orientation = destination_orientation
#                     )
#                     set_metadata!(simplified_graph, simplified_graph_edge, :orientations, edge_orientations)
#                 end
#             end
#         end
#     end
#     return simplified_graph
# end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# Description

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
# function simplify_kmer_graph(kmer_graph)
#     @info "simplifying kmer graph"
#     @info "resolving untigs..."
#     @time untigs = Mycelia.resolve_untigs(kmer_graph)

#     @info "determining untig orientations..."
#     @time oriented_untigs = Mycelia.determine_oriented_untigs(kmer_graph, untigs)

#     simplified_graph = MetaGraphs.MetaDiGraph(length(oriented_untigs))
#     MetaGraphs.set_prop!(simplified_graph, :k, kmer_graph.gprops[:k])
#     @info "initializing graph node metadata"
#     ProgressMeter.@showprogress for (vertex, untig) in enumerate(oriented_untigs)
#         MetaGraphs.set_prop!(simplified_graph, vertex, :sequence, untig.sequence)
#         MetaGraphs.set_prop!(simplified_graph, vertex, :path, untig.path)
#         MetaGraphs.set_prop!(simplified_graph, vertex, :orientations, untig.orientations)
#         MetaGraphs.set_prop!(simplified_graph, vertex, :weight, untig.weight)
#     end

#     # determine oriented edges of simplified graph
#     simplified_untigs = Vector{Pair{Pair{Int64,Bool},Pair{Int64,Bool}}}(undef, length(Graphs.vertices(simplified_graph)))
#     @info "creating simplified unitgs to help resolve connections"
#     # use a pre-allocated array here to speed up
#     ProgressMeter.@showprogress for vertex in Graphs.vertices(simplified_graph)
#         in_kmer = simplified_graph.vprops[vertex][:path][1] => simplified_graph.vprops[vertex][:orientations][1]
#         out_kmer = simplified_graph.vprops[vertex][:path][end] => simplified_graph.vprops[vertex][:orientations][end]
#     #     @show vertex, in_kmer, out_kmer
#         simplified_untigs[vertex] = in_kmer => out_kmer
#     #     push!(simplified_untigs, )
#     end

#     # make a dictionary mapping endcap to oriented_untig index

#     end_mer_map = Dict()
#     ProgressMeter.@showprogress for (i, oriented_untig) in enumerate(oriented_untigs)
#         end_mer_map[first(oriented_untig.path)] = i
#         end_mer_map[last(oriented_untig.path)] = i
#     end

#     ProgressMeter.@showprogress for (untig_index, oriented_untig) in enumerate(oriented_untigs)
#     #     @show untig_index
#         true_in_overlap = oriented_untig.sequence[1:simplified_graph.gprops[:k]-1]

#         non_backtracking_neighbors = Graphs.neighbors(kmer_graph, oriented_untig.path[1])
#         if length(oriented_untig.path) > 1
#             non_backtracking_neighbors = setdiff(non_backtracking_neighbors, oriented_untig.path[2])
#         end
#         for non_backtracking_neighbor in non_backtracking_neighbors
#             neighboring_untig_index = end_mer_map[non_backtracking_neighbor]
#             neighboring_untig = oriented_untigs[neighboring_untig_index]

#             neighbor_true_out_overlap = neighboring_untig.sequence[end-simplified_graph.gprops[:k]+2:end]
#             if neighbor_true_out_overlap == true_in_overlap
#                 e = Graphs.Edge(neighboring_untig_index, untig_index)
#     #             o = true => true
#                 o = (source_orientation = true, destination_orientation = true)
#                 Graphs.add_edge!(simplified_graph, e)
#                 Mycelia.set_metadata!(simplified_graph, e, :orientations, o)
#             end

#             neighbor_false_out_overlap = BioSequences.reverse_complement(neighboring_untig.sequence)[end-simplified_graph.gprops[:k]+2:end]
#             if neighbor_false_out_overlap == true_in_overlap
#                 e = Graphs.Edge(neighboring_untig_index, untig_index)
#     #             o = false => true
#                 o = (source_orientation = false, destination_orientation = true)
#                 Graphs.add_edge!(simplified_graph, e)
#                 Mycelia.set_metadata!(simplified_graph, e, :orientations, o)
#             end
#         end

#         true_out_overlap = oriented_untig.sequence[end-simplified_graph.gprops[:k]+2:end]
#         non_backtracking_neighbors = Graphs.neighbors(kmer_graph, oriented_untig.path[end])
#         if length(oriented_untig.path) > 1
#             non_backtracking_neighbors = setdiff(non_backtracking_neighbors, oriented_untig.path[end-1])
#         end
#         for non_backtracking_neighbor in non_backtracking_neighbors
#             neighboring_untig_index = end_mer_map[non_backtracking_neighbor]
#             neighboring_untig = oriented_untigs[neighboring_untig_index]

#             neighbor_true_in_overlap = neighboring_untig.sequence[1:simplified_graph.gprops[:k]-1]
#             if true_out_overlap == neighbor_true_in_overlap
#                 e = Graphs.Edge(untig_index, neighboring_untig_index)
#     #             o = true => true
#                 o = (source_orientation = true, destination_orientation = true)
#                 Graphs.add_edge!(simplified_graph, e)
#                 Mycelia.set_metadata!(simplified_graph, e, :orientations, o)
#             end

#             neighbor_false_in_overlap = BioSequences.reverse_complement(neighboring_untig.sequence)[1:simplified_graph.gprops[:k]-1]
#             if true_out_overlap == neighbor_false_in_overlap
#                 e = Graphs.Edge(untig_index, neighboring_untig_index)
#     #             o = true => false
#                 o = (source_orientation = true, destination_orientation = false)
#                 Graphs.add_edge!(simplified_graph, e)
#                 Mycelia.set_metadata!(simplified_graph, e, :orientations, o)
#             end
#         end

#         false_in_overlap = BioSequences.reverse_complement(oriented_untig.sequence)[1:simplified_graph.gprops[:k]-1]

#         non_backtracking_neighbors = Graphs.neighbors(kmer_graph, oriented_untig.path[end])
#         if length(oriented_untig.path) > 1
#             non_backtracking_neighbors = setdiff(non_backtracking_neighbors, oriented_untig.path[end-1])
#         end
#         for non_backtracking_neighbor in non_backtracking_neighbors
#             neighboring_untig_index = end_mer_map[non_backtracking_neighbor]
#             neighboring_untig = oriented_untigs[neighboring_untig_index]

#             neighbor_true_out_overlap = neighboring_untig.sequence[end-simplified_graph.gprops[:k]+2:end]
#             if neighbor_true_out_overlap == false_in_overlap
#                 e = Graphs.Edge(neighboring_untig_index, untig_index)
#     #             o = true => false
#                 o = (source_orientation = true, destination_orientation = false)
#                 Graphs.add_edge!(simplified_graph, e)
#                 Mycelia.set_metadata!(simplified_graph, e, :orientations, o)
#             end

#             neighbor_false_out_overlap = BioSequences.reverse_complement(neighboring_untig.sequence)[end-simplified_graph.gprops[:k]+2:end]
#             if neighbor_false_out_overlap == false_in_overlap
#                 e = Graphs.Edge(neighboring_untig_index, untig_index)
#     #             o = false => false
#                 o = (source_orientation = false, destination_orientation = false)
#                 Graphs.add_edge!(simplified_graph, e)
#                 Mycelia.set_metadata!(simplified_graph, e, :orientations, o)
#             end
#         end

#         false_out_overlap = BioSequences.reverse_complement(oriented_untig.sequence)[end-simplified_graph.gprops[:k]+2:end]

#         non_backtracking_neighbors = Graphs.neighbors(kmer_graph, oriented_untig.path[1])
#         if length(oriented_untig.path) > 1
#             non_backtracking_neighbors = setdiff(non_backtracking_neighbors, oriented_untig.path[2])
#         end

#         for non_backtracking_neighbor in non_backtracking_neighbors
#             neighboring_untig_index = end_mer_map[non_backtracking_neighbor]
#             neighboring_untig = oriented_untigs[neighboring_untig_index]

#             neighbor_true_in_overlap = neighboring_untig.sequence[1:simplified_graph.gprops[:k]-1]
#             if false_out_overlap == neighbor_true_in_overlap
#                 e = Graphs.Edge(untig_index, neighboring_untig_index)
#     #             o = false => true
#                 o = (source_orientation = false, destination_orientation = true)
#                 Graphs.add_edge!(simplified_graph, e)
#                 Mycelia.set_metadata!(simplified_graph, e, :orientations, o)
#             end

#             neighbor_false_in_overlap = BioSequences.reverse_complement(neighboring_untig.sequence)[1:simplified_graph.gprops[:k]-1]
#             if false_out_overlap == neighbor_false_in_overlap
#                 e = Graphs.Edge(untig_index, neighboring_untig_index)
#     #             o = false => false
#                 o = (source_orientation = false, destination_orientation = false)
#                 Graphs.add_edge!(simplified_graph, e)
#                 Mycelia.set_metadata!(simplified_graph, e, :orientations, o)
#             end
#         end
#     end
#     return simplified_graph
# end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# Description

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
function determine_oriented_untigs(kmer_graph, untigs)
    oriented_untigs = []
    for path in untigs
        sequence = BioSequences.LongDNASeq(kmer_graph.vprops[first(path)][:kmer])
        if length(path) == 1
            orientations = [true]
        elseif length(path) > 1
            initial_edge = Graphs.Edge(path[1], path[2])
            initial_orientation = first(kmer_graph.eprops[initial_edge][:orientations]).source_orientation
            orientations = [initial_orientation]
            if !initial_orientation
                sequence = BioSequences.reverse_complement(sequence)
            end

            for (src, dst) in zip(path[1:end-1], path[2:end])
                edge = Graphs.Edge(src, dst)
                destination = BioSequences.LongDNASeq(kmer_graph.vprops[edge.dst][:kmer])
                destination_orientation = first(kmer_graph.eprops[edge][:orientations]).destination_orientation
                push!(orientations, destination_orientation)
                if !destination_orientation
                    destination = BioSequences.reverse_complement(destination)
                end
                sequence_suffix = sequence[end-length(destination)+2:end]
                destination_prefix = destination[1:end-1]
                @assert sequence_suffix == destination_prefix
                push!(sequence, destination[end])
            end
        end

        oriented_untig =
        (
            sequence = BioSequences.canonical(sequence),
            path = BioSequences.iscanonical(sequence) ? path : reverse(path),
            orientations = BioSequences.iscanonical(sequence) ? orientations : reverse(.!orientations),
            weight = Statistics.median([kmer_graph.vprops[v][:weight] for v in path])
        )

        push!(oriented_untigs, oriented_untig)
    end
    return oriented_untigs
end

# """
# $(DocStringExtensions.TYPEDSIGNATURES)

# Description

# ```jldoctest
# julia> 1 + 1
# 2
# ```
# """
function resolve_untigs(kmer_graph)
    untigs = Vector{Int}[]
    visited = falses(Graphs.nv(kmer_graph))
    first_unvisited = findfirst(!, visited)
    while first_unvisited != nothing
        forward_walk = oriented_unbranching_walk(kmer_graph, first_unvisited, true)
        reverse_walk = oriented_unbranching_walk(kmer_graph, first_unvisited, false)
        inverted_reverse_walk = [Graphs.Edge(e.dst, e.src) for e in reverse(reverse_walk)]
        edges = vcat(inverted_reverse_walk, forward_walk)
        if isempty(edges)
            untig = [first_unvisited]
        else
            untig = vcat([first(edges).src], [edge.dst for edge in edges])
        end
        push!(untigs, untig)
        for vertex in untig
            visited[vertex] = true
        end
        first_unvisited = findfirst(!, visited)
    end
    return untigs
end

function apply_kmedoids_treshold(graph)
    kmer_counts = [MetaGraphs.get_prop(graph, v, :count) for v in Graphs.vertices(graph)]

    kmer_counts_histogram = sort(collect(StatsBase.countmap(values(kmer_counts))), by=x->x[1])

#     scale = 250
#     p = Mycelia.plot_kmer_frequency_spectra(values(kmer_counts), size=(2scale,scale), log_scale=log2, title="kmer frequencies")
#     display(p)

#     p = StatsPlots.scatter(log2.(first.(kmer_counts_histogram)))
#     display(p)

    kmer_depth_of_coverage_bins = log2.(first.(kmer_counts_histogram))

    distance_matrix = zeros((length(kmer_depth_of_coverage_bins), length(kmer_depth_of_coverage_bins)))
    for (row, depth_of_coverage_bin_1) in enumerate(kmer_depth_of_coverage_bins)
        for (col, depth_of_coverage_bin_2) in enumerate(kmer_depth_of_coverage_bins)
            distance = abs(depth_of_coverage_bin_1 - depth_of_coverage_bin_2)
            distance_matrix[row, col] = distance
        end
    end
    distance_matrix

    # max out k at the same max k we use for DNAMers
    max_k = min(length(kmer_depth_of_coverage_bins), 63)
    ks = Primes.primes(2, max_k)
    ys = map(k ->
                Statistics.mean(Statistics.mean(Clustering.silhouettes(Clustering.kmedoids(distance_matrix, k), distance_matrix)) for i in 1:100),
                ks)

    p = StatsPlots.plot(ks, ys, label="silhouette score", ylabel = "silhouette score", xlabel = "number of clusters")
    display(p)

    ymax, ymax_index = findmax(ys)
    optimal_k = ks[ymax_index]
    clusterings = [Clustering.kmedoids(distance_matrix, optimal_k) for i in 1:10]
    max_value, max_value_index = findmax(clustering -> Statistics.mean(Clustering.silhouettes(clustering, distance_matrix)), clusterings)
    optimal_clustering = clusterings[max_value_index]
    # optimal_clustering.assignments
    min_medoid_value, min_medoid_index = findmin(optimal_clustering.medoids)
    indices_to_include = map(assignment -> assignment .!= min_medoid_index, optimal_clustering.assignments)
    # kmer_depth_of_coverage_bins
    threshold = Int(ceil(2^maximum(kmer_depth_of_coverage_bins[.!indices_to_include]))) + 1

    scale = 250
    p = Mycelia.plot_kmer_frequency_spectra(values(kmer_counts), log_scale = log2, size=(2scale,scale), title="kmer frequencies")
    StatsPlots.vline!(p, log2.([threshold]))
    display(p)

    # find all vertices with count > threshold
    vertices_to_keep = [v for v in Graphs.vertices(graph) if (MetaGraphs.get_prop(graph, v, :count) > threshold)]
    # induce subgraph
    induced_subgraph, vertex_map = Graphs.induced_subgraph(graph, vertices_to_keep)

    # set kmer as indexing prop
    MetaGraphs.set_indexing_prop!(induced_subgraph, :kmer)
    return induced_subgraph
end

max_stars_repo_name: cjprybol/Eisenia

max_issues_repo_name: cjprybol/Eisenia

max_forks_repo_name: cjprybol/Eisenia


text:
using MuZero

#=
@boardgame Tac (3,3)
@actions 1:9
@pieces [
    (name=:X, icon=(color="black")),
    (name=:O, icon=(color="red"))
]
=#

struct TacState <: AbstractBoardState
    board_state::BoardState{9}
    to_play::Integer
end

TacState() = TacState(BoardState{9}(), PLAYER_1)
TacState(encoded::String) = decode(TacState, encoded)
TacState(::Vector{Nothing}, to_play::Integer) = TacState()
TacState(pieces::Vector{Piece}, to_play::Integer) = init_state(TacState, pieces, to_play)
TacState(pieces::Vector{Union{Nothing, Piece}}, to_play::Integer) = init_state(TacState, pieces, to_play)

#@todo change this to player "A" and "B"
const PLAYER_1_PIECES = [Piece(:X, "X", icon_symbol="X", icon_color="black", icon_size=60)]
const PLAYER_2_PIECES = [Piece(:O, "O", icon_symbol="O", icon_color="red", icon_size=60)]
const ALL_PIECES = [PLAYER_1_PIECES..., PLAYER_2_PIECES...]

MuZero.get_pieces(::Type{TacState}) = ALL_PIECES
MuZero.num_rows(::Type{TacState}) = 3
MuZero.num_cols(::Type{TacState}) = 3
MuZero.image_scale(::Type{TacState}) = 750
MuZero.image_title_offset(::Type{TacState}) = (75,45)
MuZero.actions(::Type{TacState}) = collect(1:9)
MuZero.num_players(::Type{TacState}) = 2

struct TacStateSVG
    s::TacState
end

struct Tac <: AbstractBoardGame
    history::Vector{<:Integer}
    states::Vector{TacState}
    rewards::Vector{Float32}
    function Tac(history, states, rewards)
        @assert length(history) == length(rewards)
        @assert length(states) == (length(history) + 1)
        new(history, states, rewards)
    end
end

function Tac()
    Tac(
        Vector{Int8}(),
        [TacState()],
        Vector{Float32}(),
    )
end

Tac(actions::Vector{<:Integer}) = init_game(Tac, actions)

MuZero.getstatetype(::Type{Tac}) = TacState

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

# @todo change to get_coord(::TacState, ::AbstractAction)
function MuZero.get_coord(s::TacState, action::Integer)
    action_to_coord = (
        a7=(1,1), a8=(1,2), a9=(1,3),
        a4=(2,1), a5=(2,2), a6=(2,3),
        a1=(3,1), a2=(3,2), a3=(3,3))
    action_to_coord[Symbol("a$action")]
end

# @todo simplify this to get_marker(::Tac, ::AbstractAction)
# @todo break out into get_piece(::Tac, ::AbstractAction)
function MuZero.get_marker(s::TacState, player::Integer, action::Integer)
    row_idx, col_idx = get_coord(s, action)
    name = player == PLAYER_1 ? :X : :O
    piece = get_piece(s, name)
    Marker(row_idx, col_idx, piece)
end

function MuZero.get_winner(s::TacState)
    for player in [PLAYER_1, PLAYER_2]
        target = player * ones(UInt8, 3)
        for row_idx in 1:num_rows(s)
            if s[row_idx, 1:3] == target
                return player
            end
        end
        for col_idx in 1:num_cols(s)
            if s[1:3, col_idx] == target
                return player
            end
        end
        if s[1, 1] == s[2, 2] == s[3, 3] == player
            return player
        end
        if s[3, 1] == s[2, 2] == s[1, 3] == player
            return player
        end
    end
    nothing
end

function MuZero.legal_actions(s::TacState)
    if game_is_won(s)
        return []
    end
    filter(a->begin
        coord = get_coord(s, a)
        iszero(s[coord...])
    end, actions(s))
end

function MuZero.getconfig(::Type{Tac})
    config = make_board_game_config(
        Tac,
        network_type=Incremental,
        max_moves=9,
        root_dirichlet_alpha=0.3,
        action_space_size=9, # Number of all possible actions in the game
        support_size = 10, # Value and reward are scaled (with almost sqrt) and encoded on a vector with a range of -support_size to support_size. Choose it so that support_size <= sqrt(max(abs(discounted reward)))
        observation_shape = (3, 3, 3),
        stacked_observations = 0, # Number of previous observations and previous actions to add to the current observation
        hidden_state_size = 32, # The 1-dimensional size of the encoded game state
        checkpoint_interval = 100, #  Number of training steps before using the model for self-playing
        num_actors = 3,
        weight_decay = 1e-3,
        lr_init=0.002f0,
        save_interval = 10000,
        num_unroll_steps=1,
    )
    config.batch_size = 10
    config.num_simulations = 10
    config.training_steps = 12000
    config
end


text:
elementTypeID_NodeCount = [ 2,3,4,4,8,6,5,3,6,9,10,27,18,14,
                            1,8,20,15,13,9,10,12,15,25,21,
                            4,5,6,20,35,56 ]

elementType_TypeID = Dict([ ("L2", 1), ("Tri3", 2), ("Quad4", 3), ("Tet4", 4),
                            ("Hex8", 5), ("L3", 8), ("Tri6", 9), ("Quad9", 10),
                            ("Tet10", 11), ("Hex27", 12)])

gmsh_tag_type= Dict([("Point", 0), ("Curve", 1), ("Surface", 2), ("Volume", 3)])

struct GMSHFileInfo{T1<:Real,T2<:Int}
    version_number::T1;
    file_type::T2;
    data_size::T2;
end

struct PhysicalEntity
    type::Int
    number::Int
    name::String
end

struct PhysicalEntities
    count::Int
    entities::Union{Nothing,Array{PhysicalEntity,1}}
end

struct Point{T1<:Int, T2<:Real}
    pointTag::T1
    x::T2
    y::T2
    z::T2
    numPhysicalTags::T1
    physicalTag::Array{T1,1}
end

struct Curve{T1<:Int, T2<:Real}
    curveTag::T1
    xmin::T2
    ymin::T2
    zmin::T2
    xmax::T2
    ymax::T2
    zmax::T2
    numPhysicalTags::T1
    physicalTag::Array{T1,1}
    numBoundingPoints::T1
    pointTag::Array{T1,1}
end

struct Surface{T1<:Int, T2<:Real}
    surfaceTag::T1
    xmin::T2
    ymin::T2
    zmin::T2
    xmax::T2
    ymax::T2
    zmax::T2
    numPhysicalTags::T1
    physicalTag::Array{T1,1}
    numBoundingCurves::T1
    curveTag::Array{T1,1}
end

struct Volume{T1<:Int, T2<:Real}
    volumeTag::T1
    xmin::T2
    ymin::T2
    zmin::T2
    xmax::T2
    ymax::T2
    zmax::T2
    numPhysicalTags::T1
    physicalTag::Array{T1,1}
    numBoundingSurfaces::T1
    surfaceTag::Array{T1,1}
end

struct MeshEntities
    points::Array{Point, 1}
    curves::Array{Curve, 1}
    surfaces::Array{Surface, 1}
    volumes::Array{Volume}
end

mutable struct Nodes{T1<:Int, T2<:Real}
    num::Vector{T1}
    x::Vector{T2}
    y::Vector{T2}
    z::Vector{T2}
    noEntities::Vector{T1}
    entityTypeID::Vector{T1}
    entityTag::Vector{T1}
end

struct Element
    number::Int
    nodes::Array{Int,1}
end

struct ElementGroup
    entityTag::Int
    entityTypeID::Int
    elementTypes::Int
    count::Int
    elements::Array{Element,1}
end

struct GMSH
    info::GMSHFileInfo
    pyhysical_entities::PhysicalEntities
    mesh_entities::MeshEntities
    nodes::Nodes
    element_groups::Array{ElementGroup,1}
end

max_stars_repo_name: c-schubert/linearMeshIO

max_issues_repo_name: c-schubert/linearMeshIO

max_forks_repo_name: c-schubert/linearMeshIO


text:
using Metalhead, Test

# Standardized testing for the models of tomorrow
@info "Starting Basic Models Tests..."
@testset "Basic Model Tests" begin
    for (T, MODEL) in [
            (Float32, VGG19),
            (Float32, SqueezeNet),
            (Float64, DenseNet),
            (Float64, GoogleNet),
        ]
	@info "Testing $(MODEL)..."
        model = MODEL()

        x_test = rand(T, 224, 224, 3, 1)
        y_test = model(x_test)

        # Test that types and shapes work out as we expect
        @test y_test isa AbstractArray
        @test length(y_test) == 1000

        # Test that the models can be indexed
        @test length(model.layers[1:4].layers) == 4
    end
end

# Test proper download and functioning of CIFAR10
@info "Starting CIFAR10 tests..."
@testset "CIFAR dataset tests" begin
    x1 = trainimgs(CIFAR10)[1]
    x2 = valimgs(CIFAR10)[1]

    # Test that the input is roughly what we expect
    @test size(x1.img) == (32, 32)
    @test size(x2.img) == (32, 32)
end

# Test printing of prediction
@info "Testing Prediction on CIFAR10..."
@testset "Prediction table display" begin
    x = valimgs(CIFAR10)[1]
    m = VGG19()
    predict(m, x)
end

# Just run the prediction code end-to-end
# TODO: Set up travis to actually run these
if length(datasets()) == 2
    for dataset in (ImageNet, CIFAR10)
        val1 = valimgs(dataset)[1]
        predict(vgg19, val1)
        classify(vgg19, val1)
    end
    predict(vgg19, testimgs(dataset(ImageNet))[1])
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const OneHotEncoder_HandleUnknown = (;[
    Symbol("ErrorOnUnknown") => Int32(0),
    Symbol("IgnoreUnknown") => Int32(1),
]...)

mutable struct OneHotEncoder <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function OneHotEncoder(; kwargs...)
        obj = new(meta(OneHotEncoder), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct OneHotEncoder
const __meta_OneHotEncoder = Ref{ProtoMeta}()
function meta(::Type{OneHotEncoder})
    ProtoBuf.metalock() do
        if !isassigned(__meta_OneHotEncoder)
            __meta_OneHotEncoder[] = target = ProtoMeta(OneHotEncoder)
            fnum = Int[1,2,10,11]
            allflds = Pair{Symbol,Union{Type,String}}[:stringCategories => StringVector, :int64Categories => Int64Vector, :outputSparse => Bool, :handleUnknown => Int32]
            oneofs = Int[1,1,0,0]
            oneof_names = Symbol[Symbol("CategoryType")]
            meta(target, OneHotEncoder, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_OneHotEncoder[]
    end
end
function Base.getproperty(obj::OneHotEncoder, name::Symbol)
    if name === :stringCategories
        return (obj.__protobuf_jl_internal_values[name])::StringVector
    elseif name === :int64Categories
        return (obj.__protobuf_jl_internal_values[name])::Int64Vector
    elseif name === :outputSparse
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :handleUnknown
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

export OneHotEncoder_HandleUnknown, OneHotEncoder


text:
# GPU memory management and pooling

using Printf
using TimerOutputs

using Base: @lock
using Base.Threads: SpinLock

# global lock for shared object dicts (allocated, requested).
# stats are not covered by this and cannot be assumed to be exact.
# each allocator needs to lock its own resources separately too.
const memory_lock = SpinLock()

const MEMDEBUG = ccall(:jl_is_memdebug, Bool, ())

## allocation statistics

mutable struct AllocStats
  # pool allocation requests
  pool_nalloc::Int
  pool_nfree::Int
  ## in bytes
  pool_alloc::Int

  # actual CUDA allocations
  actual_nalloc::Int
  actual_nfree::Int
  ## in bytes
  actual_alloc::Int
  actual_free::Int

  pool_time::Float64
  actual_time::Float64
end

const alloc_stats = AllocStats(0, 0, 0, 0, 0, 0, 0, 0, 0)

Base.copy(alloc_stats::AllocStats) =
  AllocStats((getfield(alloc_stats, field) for field in fieldnames(AllocStats))...)

## CUDA allocator

const alloc_to = TimerOutput()

"""
    alloc_timings()

Show the timings of the CUDA allocator. Assumes [`CuArrays.enable_timings()`](@ref) has been
called.
"""
alloc_timings() = (show(alloc_to; allocations=false, sortby=:name); println())

const usage = Threads.Atomic{Int}(0)
const usage_limit = Ref{Union{Nothing,Int}}(nothing)

const allocated = Dict{CuPtr{Nothing},Mem.DeviceBuffer}()

function actual_alloc(bytes)
  # check the memory allocation limit
  if usage_limit[] !== nothing
    if usage[] + bytes > usage_limit[]
      return nothing
    end
  end

  # try the actual allocation
  time, buf = try
    time = Base.@elapsed begin
      @timeit_debug alloc_to "alloc" buf = Mem.alloc(Mem.Device, bytes)
    end
    Threads.atomic_add!(usage, bytes)
    time, buf
  catch err
    (isa(err, CuError) && err.code == CUDAdrv.ERROR_OUT_OF_MEMORY) || rethrow()
    return nothing
  end
  @assert sizeof(buf) == bytes
  ptr = convert(CuPtr{Nothing}, buf)

  # record the buffer
  @lock memory_lock begin
    @assert !haskey(allocated, ptr)
    allocated[ptr] = buf
  end

  alloc_stats.actual_time += time
  alloc_stats.actual_nalloc += 1
  alloc_stats.actual_alloc += bytes

  return ptr
end

function actual_free(ptr::CuPtr{Nothing})
  # look up the buffer
  buf = @lock memory_lock begin
    buf = allocated[ptr]
    delete!(allocated, ptr)
    buf
  end
  bytes = sizeof(buf)

  # free the memory
  @timeit_debug alloc_to "free" begin
    time = Base.@elapsed Mem.free(buf)
    Threads.atomic_sub!(usage, bytes)
  end

  alloc_stats.actual_time += time
  alloc_stats.actual_nfree += 1
  alloc_stats.actual_free += bytes

  return
end

## memory pools

const pool_to = TimerOutput()

macro pool_timeit(args...)
    TimerOutputs.timer_expr(CuArrays, true, :($CuArrays.pool_to), args...)
end

"""
    pool_timings()

Show the timings of the currently active memory pool. Assumes
[`CuArrays.enable_timings()`](@ref) has been called.
"""
pool_timings() = (show(pool_to; allocations=false, sortby=:name); println())

# pool API:
# - init()
# - alloc(sz)::CuPtr{Nothing}
# - free(::CuPtr{Nothing})
# - reclaim(nb::Int=typemax(Int))::Int
# - used_memory()
# - cached_memory()

include("memory/binned.jl")
include("memory/simple.jl")
include("memory/split.jl")
include("memory/dummy.jl")

const pool = Ref{Module}(BinnedPool)

## interface

export OutOfGPUMemoryError

const requested = Dict{CuPtr{Nothing},Vector}()

"""
    OutOfGPUMemoryError()

An operation allocated too much GPU memory for either the system or the memory pool to
handle properly.
"""
struct OutOfGPUMemoryError <: Exception
  sz::Int
end

function Base.showerror(io::IO, err::OutOfGPUMemoryError)
    println(io, "Out of GPU memory trying to allocate $(Base.format_bytes(err.sz))")
    memory_status(io)
end

"""
    alloc(sz)

Allocate a number of bytes `sz` from the memory pool. Returns a `CuPtr{Nothing}`; may throw
a [`OutOfGPUMemoryError`](@ref) if the allocation request cannot be satisfied.
"""
@inline function alloc(sz)
  # 0-byte allocations shouldn't hit the pool
  sz == 0 && return CU_NULL

  time = Base.@elapsed begin
    @pool_timeit "pooled alloc" ptr = pool[].alloc(sz)::Union{Nothing,CuPtr{Nothing}}
  end
  ptr === nothing && throw(OutOfGPUMemoryError(sz))

  # record the allocation
  if Base.JLOptions().debug_level >= 2
    bt = backtrace()
    @lock memory_lock begin
      @assert !haskey(requested, ptr)
      requested[ptr] = bt
    end
  end

  alloc_stats.pool_time += time
  alloc_stats.pool_nalloc += 1
  alloc_stats.pool_alloc += sz

  if MEMDEBUG && ptr == CuPtr{Cvoid}(0xbbbbbbbbbbbbbbbb)
    error("Allocated a scrubbed pointer")
  end

  return ptr
end

"""
    free(sz)

Releases a buffer pointed to by `ptr` to the memory pool.
"""
@inline function free(ptr::CuPtr{Nothing})
  # 0-byte allocations shouldn't hit the pool
  ptr == CU_NULL && return

  if MEMDEBUG && ptr == CuPtr{Cvoid}(0xbbbbbbbbbbbbbbbb)
    Core.println("Freeing a scrubbed pointer!")
  end

  # this function is typically called from a finalizer, where we can't switch tasks,
  # so perform our own error handling.
  try
    # record the allocation
    if Base.JLOptions().debug_level >= 2
      @lock memory_lock begin
        @assert haskey(requested, ptr)
        delete!(requested, ptr)
      end
    end

    time = Base.@elapsed begin
      @pool_timeit "pooled free" pool[].free(ptr)
    end

    alloc_stats.pool_time += time
    alloc_stats.pool_nfree += 1
  catch ex
    Base.showerror_nostdio(ex, "WARNING: Error while freeing $ptr")
    Base.show_backtrace(Core.stdout, catch_backtrace())
    Core.println()
  end

  return
end

"""
    reclaim([sz=typemax(Int)])

Reclaims `sz` bytes of cached memory. Use this to free GPU memory before calling into
functionality that does not use the CuArrays memory pool. Returns the number of bytes
actually reclaimed.
"""
reclaim(sz::Int=typemax(Int)) = pool[].reclaim(sz)

"""
    @retry_reclaim fail ex

Run a block of code `ex` repeatedly until it successfully allocates the memory it needs;
Failure to do so indicated by returning `fail`. At each try, more and more memory is freed
from the CuArrays memory pool. When that is not possible anymore, `fail` will be returned.

This macro is intended for use with CUDA APIs, which sometimes allocate (outside of the
CuArrays memory pool) and return a specific error code when failing to.
"""
macro retry_reclaim(fail, ex)
  quote
    ret = nothing
    for phase in 1:3
      ret = $(esc(ex))
      ret == $(esc(fail)) || break

      # incrementally more costly reclaim of cached memory
      if phase == 1
        reclaim()
      elseif phase == 2
        GC.gc(false)
        reclaim()
      elseif phase == 3
        GC.gc(true)
        reclaim()
      end
    end
    ret
  end
end

## utilities

"""
    @allocated

A macro to evaluate an expression, discarding the resulting value, instead returning the
total number of bytes allocated during evaluation of the expression.
"""
macro allocated(ex)
    quote
        let
            local f
            function f()
                b0 = alloc_stats.pool_alloc
                $(esc(ex))
                alloc_stats.pool_alloc - b0
            end
            f()
        end
    end
end

"""
    @time ex

Run expression `ex` and report on execution time and GPU/CPU memory behavior. The GPU is
synchronized right before and after executing `ex` to exclude any external effects.

"""
macro time(ex)
    quote
        # @time might surround an application, so be sure to initialize CUDA before that
        CUDAnative.prepare_cuda_call()

        # coarse synchronization to exclude effects from previously-executed code
        CUDAdrv.synchronize()

        local gpu_mem_stats0 = copy(alloc_stats)
        local cpu_mem_stats0 = Base.gc_num()
        local cpu_time0 = time_ns()

        # fine-grained synchronization of the code under analysis
        local val = @sync $(esc(ex))

        local cpu_time1 = time_ns()
        local cpu_mem_stats1 = Base.gc_num()
        local gpu_mem_stats1 = copy(alloc_stats)

        local cpu_time = (cpu_time1 - cpu_time0) / 1e9
        local gpu_gc_time = gpu_mem_stats1.pool_time - gpu_mem_stats0.pool_time
        local gpu_alloc_count = gpu_mem_stats1.pool_nalloc - gpu_mem_stats0.pool_nalloc
        local gpu_lib_time = gpu_mem_stats1.actual_time - gpu_mem_stats0.actual_time
        local gpu_alloc_size = gpu_mem_stats1.pool_alloc - gpu_mem_stats0.pool_alloc
        local cpu_mem_stats = Base.GC_Diff(cpu_mem_stats1, cpu_mem_stats0)
        local cpu_gc_time = cpu_mem_stats.total_time / 1e9
        local cpu_alloc_count = Base.gc_alloc_count(cpu_mem_stats)
        local cpu_alloc_size = cpu_mem_stats.allocd

        Printf.@printf("%10.6f seconds", cpu_time)
        for (typ, gctime, libtime, bytes, allocs) in
            (("CPU", cpu_gc_time, 0, cpu_alloc_size, cpu_alloc_count),
             ("GPU", gpu_gc_time, gpu_lib_time, gpu_alloc_size, gpu_alloc_count))
          if bytes != 0 || allocs != 0
              allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units), Int64(1000))
              if ma == 1
                  Printf.@printf(" (%d%s %s allocation%s: ", allocs, Base._cnt_units[ma], typ, allocs==1 ? "" : "s")
              else
                  Printf.@printf(" (%.2f%s %s allocations: ", allocs, Base._cnt_units[ma], typ)
              end
              print(Base.format_bytes(bytes))
              if gctime > 0
                  Printf.@printf(", %.2f%% gc time", 100*gctime/cpu_time)
                if libtime > 0
                    Printf.@printf(" of which %.2f%% spent allocating", 100*libtime/gctime)
                end
              end
              print(")")
          elseif gctime > 0
              Printf.@printf(", %.2f%% %s gc time", 100*gctime/cpu_time, typ)
          end
        end
        println()

        val
    end
end

"""
    memory_status([io=stdout])

Report to `io` on the memory status of the current GPU and the active memory pool.
"""
function memory_status(io::IO=stdout)
  free_bytes, total_bytes = CUDAdrv.Mem.info()
  used_bytes = total_bytes - free_bytes
  used_ratio = used_bytes / total_bytes
  @printf(io, "Effective GPU memory usage: %.2f%% (%s/%s)\n",
              100*used_ratio, Base.format_bytes(used_bytes),
              Base.format_bytes(total_bytes))

  @printf(io, "CuArrays allocator usage: %s", Base.format_bytes(usage[]))
  if usage_limit[] !== nothing
    @printf(io, " (capped at %s)", Base.format_bytes(usage_limit[]))
  end
  println(io)

  alloc_used_bytes = pool[].used_memory()
  alloc_cached_bytes = pool[].cached_memory()
  alloc_total_bytes = alloc_used_bytes + alloc_cached_bytes
  @printf(io, "%s usage: %s (%s allocated, %s cached)\n", nameof(pool[]),
              Base.format_bytes(alloc_total_bytes), Base.format_bytes(alloc_used_bytes),
              Base.format_bytes(alloc_cached_bytes))

  # check if the memory usage as counted by the CUDA allocator wrapper
  # matches what is reported by the pool implementation
  discrepancy = abs(usage[] - alloc_total_bytes)
  if discrepancy != 0
    println(io, "Discrepancy of $(Base.format_bytes(discrepancy)) between memory pool and allocator!")
  end

  if Base.JLOptions().debug_level >= 2
    requested′, allocated′ = @lock memory_lock begin
      copy(requested), copy(allocated)
    end
    for (ptr, bt) in requested′
      buf = allocated′[ptr]
      @printf(io, "\nOutstanding memory allocation of %s at %p",
              Base.format_bytes(sizeof(buf)), Int(ptr))
      stack = stacktrace(bt, false)
      StackTraces.remove_frames!(stack, :alloc)
      Base.show_backtrace(io, stack)
      println(io)
    end
  end
end

## init

"""
    enable_timings()

Enable the recording of debug timings.
"""
enable_timings() = (TimerOutputs.enable_debug_timings(CuArrays); return)
disable_timings() = (TimerOutputs.disable_debug_timings(CuArrays); return)

function __init_memory__()
  # memory limit configuration
  memory_limit_str = if haskey(ENV, "JULIA_CUDA_MEMORY_LIMIT")
    ENV["JULIA_CUDA_MEMORY_LIMIT"]
  elseif haskey(ENV, "CUARRAYS_MEMORY_LIMIT")
    Base.depwarn("The CUARRAYS_MEMORY_LIMIT environment flag is deprecated, please use JULIA_CUDA_MEMORY_LIMIT instead.", :__init_memory__)
    ENV["CUARRAYS_MEMORY_LIMIT"]
  else
    nothing
  end
  if memory_limit_str !== nothing
    usage_limit[] = parse(Int, memory_limit_str)
  end

  # memory pool configuration
  memory_pool_str = if haskey(ENV, "JULIA_CUDA_MEMORY_POOL")
    ENV["JULIA_CUDA_MEMORY_POOL"]
  elseif haskey(ENV, "CUARRAYS_MEMORY_POOL")
    Base.depwarn("The CUARRAYS_MEMORY_POOL environment flag is deprecated, please use JULIA_CUDA_MEMORY_POOL instead.", :__init_memory__)
    ENV["CUARRAYS_MEMORY_POOL"]
  else
    nothing
  end
  if memory_pool_str !== nothing
    pool[] =
      if memory_pool_str == "binned"
        BinnedPool
      elseif memory_pool_str == "simple"
        SimplePool
      elseif memory_pool_str == "split"
        SplittingPool
      elseif memory_pool_str == "none"
        DummyPool
      else
        error("Invalid allocator selected")
      end

    # the user hand-picked an allocator, so be a little verbose
    atexit(()->begin
      Core.println("""
        CuArrays.jl $(nameof(pool[])) statistics:
         - $(alloc_stats.pool_nalloc) pool allocations: $(Base.format_bytes(alloc_stats.pool_alloc)) in $(round(alloc_stats.pool_time; digits=2))s
         - $(alloc_stats.actual_nalloc) CUDA allocations: $(Base.format_bytes(alloc_stats.actual_alloc)) in $(round(alloc_stats.actual_time; digits=2))s""")
    end)
  end

  # initialization
  pool[].init()
  reset_timers!()
end

"""
    reset_timers!()

Reset all debug timers. This is automatically called at initialization time,
"""
function reset_timers!()
  TimerOutputs.reset_timer!(alloc_to)
  TimerOutputs.reset_timer!(pool_to)
end


text:
module BenchDivideAtMeanWithMap
using Random
using BenchmarkTools
using TapirBenchmarks

Random.seed!(1234)

const SUITE = BenchmarkGroup()

@inline square(x) = x^2

for f in [sin, square]
    n = 2^20
    SUITE[:f=>Symbol(f)] = s0 = BenchmarkGroup()
    xs = rand(n)
    s0[:impl=>:seq] = @benchmarkable divide_at_mean_seq($f, $xs)
    s0[:impl=>:threads] = @benchmarkable divide_at_mean_threads($f, $xs)
    s0[:impl=>:tapir] = @benchmarkable divide_at_mean_tapir($f, $xs)
end

end
BenchDivideAtMeanWithMap.SUITE


text:
#=
mkdir /dev/cpuset
mount -t cpuset cpuset /dev/cpuset
cd /dev/cpuset
mkdir default
mkdir rr
cd default
echo 0-6 > cpuset.cpus
while read i; do /bin/echo $i; done < ../tasks > tasks
cd ..
cd rr
echo 7 > cpuset.cpus
=#
using RR; using Cxx
using RR: current_session
session = icxx"""rr::ReplaySession::create("");"""
timeline = icxx"""new rr::ReplayTimeline{std::move($session),rr::ReplaySession::Flags{}};""";
session = nothing
icxx"$timeline->maybe_add_reverse_exec_checkpoint(rr::ReplayTimeline::LOW_OVERHEAD);"
open("/dev/cpuset/rr/tasks","w") do f
    print(f,icxx"$(current_task(current_session(timeline)))->real_tgid();")
end

function check_consistency()
    session = current_session(timeline)
    task = current_task(session)
    ticks = icxx"$task->tick_count();"
    traceticks = icxx"$session->current_trace_frame().ticks();"
    @show (ticks, traceticks)
end

cxx"""
#include <cxxabi.h>
#include "llvm/MC/SubtargetFeature.h"
#include "llvm/MC/MCObjectFileInfo.h"
#include "llvm/MC/MCContext.h"
#include "llvm/MC/MCDisassembler.h"
#include "llvm/MC/MCInst.h"
#include "llvm/MC/MCInstPrinter.h"
#include "llvm/MC/MCAsmInfo.h"
#include "llvm/Support/TargetSelect.h"
#include "llvm/Support/TargetRegistry.h"
#include "llvm/MC/MCInstrInfo.h"
using namespace llvm;
"""

include(Pkg.dir("RR","scratch","disassembler.jl"))

const ctx = DisAsmContext()

for i = 1:100
  RR.single_step!(current_session(timeline))
end

function do_and_check_single!()
    task = current_task(current_session(timeline))
    regs = icxx"$task->regs();"
    ip = Gallium.ip(regs)
    insts = RR.load(current_task(current_session(timeline)),
      Gallium.RemotePtr{UInt8}(UInt(ip)), 20)
    @show insts
    @show UInt(ip)
    (Inst, InstSize) = getInstruction(insts, 0; ctx = ctx)
    @show (Inst, InstSize)
    @show icxx"$Inst->getOpcode();"
    ticks = icxx"$task->tick_count();"
    iscondbranch = isConditionalBranch(Inst, ctx)
    free(Inst)
    RR.single_step!(current_session(timeline))
    @assert icxx"$task->tick_count();" == ticks + iscondbranch
end


text:
@testset "Blue Style" begin
    @testset "ops" begin
        # `//` and `^` are binary ops without whitespace around them
        @test bluefmt("1 // 2 + 3 ^ 4") == "1//2 + 3^4"
        @test bluefmt("a.//10") == "a .// 10"
    end

    @testset "nest to one line" begin
        str_ = """
        var = [arg1, #com
        arg2]
        """
        str = """
        var = [
            arg1, #com
            arg2,
        ]
        """
        @test bluefmt(str_) == str

        str_ = """
        var = (arg1,
        arg2)
        """
        str = """
        var = (arg1, arg2)
        """
        @test bluefmt(str_, 4, 18) == str

        str_ = """
        var = {arg1,
        arg2}
        """
        str = """
        var = {
            arg1, arg2
        }
        """
        @test bluefmt(str_, 4, 17) == str
        @test bluefmt(str_, 4, 14) == str

        str_ = """
        var = call(arg1,
        arg2)
        """
        str = """
        var = call(
            arg1, arg2
        )
        """
        @test bluefmt(str_, 4, 14) == str

        str = """
        var = call(
            arg1,
            arg2,
        )
        """
        @test bluefmt(str_, 4, 13) == str

        str_ = """
        var = ref[arg1,
        arg2]
        """
        str = """
        var = ref[
            arg1, arg2
        ]
        """
        @test bluefmt(str_, 4, 14) == str

        str = """
        var = ref[
            arg1,
            arg2,
        ]
        """
        @test bluefmt(str_, 4, 13) == str

        str_ = """
        var = ABC{arg1,
        arg2}
        """
        str = """
        var = ABC{
            arg1,arg2
        }
        """
        @test bluefmt(str_, 4, 13) == str

        str = """
        var = ABC{
            arg1,
            arg2,
        }
        """
        @test bluefmt(str_, 4, 12) == str

        str_ = """
        var = @call(arg1,
        arg2)
        """
        str = """
        var = @call(
            arg1, arg2
        )
        """
        @test bluefmt(str_, 4, 14) == str

        str = """
        var = @call(
            arg1,
            arg2
        )
        """
        @test bluefmt(str_, 4, 13) == str

        str_ = """
        function long_name_of_function_because_i_am_writing_an_example(
            arg1, arg2, arg3, arg4, arg5, arg6
        )
            # code
        end
        """
        @test bluefmt(str_, 4, 38) == str_

        str = """
        function long_name_of_function_because_i_am_writing_an_example(
            arg1,
            arg2,
            arg3,
            arg4,
            arg5,
            arg6,
        )
            # code
        end
        """
        @test bluefmt(str_, 4, 37) == str

        str = """
        Dict(
            "options" => join(
                imap(Iterators.filter(keep_option, connection_options)) do (k, v)
                    "-c k=(show_option(v))"
                end,
                " ",
            ),
        )
        """
        @test bluefmt(str, 4, 92) == str

        str_ = """
        var = foo(
            map(arr) do x
                x * 10
            end, "")
        """
        str = """
        var = foo(
            map(arr) do x
                x * 10
            end,
            "",
        )
        """
        @test bluefmt(str_, 4, 92) == str
    end

    @testset "do not prepend return in `do` blocks" begin
        str = """
        map(arg1, arg2) do x, y
            expr1
            expr2
        end"""
        @test bluefmt(str, always_use_return = true) == str
    end

    @testset "separate kw args with semicolon" begin
        str_ = "xy = f(x, y=3)"
        str = "xy = f(x; y = 3)"
        @test bluefmt(str_) == str

        str_ = "xy = f(x=1, y=2)"
        str = "xy = f(; x = 1, y = 2)"
        @test bluefmt1(str_) == str
        @test bluefmt(str_) == str
        @test bluefmt(str) == str

        str_ = "xy = f(x=1; y=2)"
        @test bluefmt1(str_) == str
        @test bluefmt(str_) == str

        str = """
        function g(x, y = 1)
            return x + y
        end
        macro h(x, y = 1)
            nothing
        end
        shortdef1(MatrixT, VectorT = nothing) = nothing
        shortdef2(MatrixT, VectorT = nothing) where {T} = nothing
        """
        @test bluefmt1(str) == str
        @test bluefmt(str) == str
    end

    @testset "no chained ternary allowed !!!" begin
        str = """
        E1 ? A : B
        """
        @test bluefmt(str) == str

        str_ = """
        E1 ? A : E2 ? B : C
        """
        str = """
        if E1
            A
        elseif E2
            B
        else
            C
        end
        """
        @test bluefmt(str_) == str
    end

    @testset "use `return nothing` instead of `return`" begin
        str_ = """
        function foo()
            return
        end
        """
        str = """
        function foo()
            return nothing
        end
        """
        @test bluefmt(str_) == str

        str_ = "a || return"
        str = "a || return nothing"
        @test bluefmt(str_) == str
    end
end


text:
module AutoAccessorMethods

export @auto_accessor_methods

using MacroTools
using MacroTools: postwalk

"""
    @auto_accessor_methods typ [ex=true]

Automatically define accessor methods for fields of a type `typ`. Disable accessor methods
from being exported by setting `ex = false`.

`@auto_accessor_methods` transforms the following code:

```julia
@auto_accessor_methods struct T
    x
    y
end
```

into:

```julia
export x, y

struct T
    x
    y
end

x(obj::T) = obj.x
y(obj::T) = obj.y
```

Examples:

```jldoctest
julia> using AutoAccessorMethods

julia> @auto_accessor_methods struct MyType
           first_field
           second_field
       end

julia> first_field
first_field (generic function with 1 method)

julia> second_field
second_field (generic function with 1 method)

julia> mt = MyType("Hello", "Goodbye");

julia> first_field(mt)
"Hello"

julia> second_field(mt)
"Goodbye"

```

```jldoctest
julia> using AutoAccessorMethods

julia> @auto_accessor_methods struct MyType{T}
           first_field::T
           second_field::Tuple{T,T}
       end

julia> mt = MyType(1, (2, 3));

julia> first_field(mt)
1

julia> second_field(mt)
(2, 3)

```
"""
macro auto_accessor_methods(expr, ex=true)
    postwalk(expr) do expr
        # Capture the type definition
        @capture(expr, (struct T_{__} fields__ end) | (struct T_ fields__ end)) || return expr

        # Fields
        length(fields) > 0 ||
            throw(ArgumentError("`$T` has no fields"))

        field_names = map(fields) do field
            @capture(field, var_::__) || @capture(field, var_)
            var
        end

        # Accessor methods
        accessors = Expr[]
        for fn in field_names
            push!(
                accessors,
                :($(esc(fn))(x::$(esc(T))) = getfield(x, Symbol($(esc(fn))))),
            )

            if ex
                push!(
                    accessors,
                    :(export $(esc(fn)))
                )
            end
        end

        return quote
            $(esc(expr))
            $(accessors...)
        end
    end
end

# Trival macro to test nested macro invocations
# TODO get this to work when defined in `test/runtests.jl`
macro noop(expr)
    esc(expr)
end

end # module