text:
# based on go implementation in src/net/http/cookie.go

# Copyright (c) 2009 The Go Authors. All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#    * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#    * Redistributions in binary form must reproduce the above
# copyright notice, this list of conditions and the following disclaimer
# in the documentation and/or other materials provided with the
# distribution.
#    * Neither the name of Google Inc. nor the names of its
# contributors may be used to endorse or promote products derived from
# this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

module Cookies

export Cookie

import Base.==
import HTTP.isurlchar

"""
    Cookie()
    Cookie(; kwargs...)
    Cookie(name, value; kwargs...)

A Cookie represents an HTTP cookie as sent in the Set-Cookie header of an
HTTP response or the Cookie header of an HTTP request. Supported fields
(which can be set using keyword arguments) include:

  * `name`: name of the cookie
  * `value`: value of the cookie
  * `path`: applicable path for the cookie
  * `domain`: applicable domain for the cookie
  * `expires`: a `DateTime` representing when the cookie should expire
  * `maxage`: `maxage == 0` means no max age, `maxage < 0` means delete cookie now, `max age > 0` means the # of seconds until expiration
  * `secure::Bool`: secure cookie attribute
  * `httponly::Bool`: httponly cookie attribute
  * `hostonly::Bool`: hostonly cookie attribute

See http:#tools.ietf.org/html/rfc6265 for details.
"""
mutable struct Cookie
    name::String
    value::String

    path::String      # optional
    domain::String    # optional
    expires::DateTime # optional

    # MaxAge=0 means no 'Max-Age' attribute specified.
    # MaxAge<0 means delete cookie now, equivalently 'Max-Age: 0'
    # MaxAge>0 means Max-Age attribute present and given in seconds
    maxage::Int
    secure::Bool
    httponly::Bool
    hostonly::Bool
    unparsed::Vector{String} # Raw text of unparsed attribute-value pairs
end

function Cookie(cookie::Cookie; kwargs...)
    for (k, v) in kwargs
        setfield!(cookie, k, convert(fieldtype(Cookie, k), v))
    end
    return cookie
end
Cookie(; kwargs...) = Cookie(Cookie("", ""); kwargs...)

Cookie(name, value; args...) = Cookie(Cookie(name, value, "", "", DateTime(), 0, false, false, false, String[]); args...)

Base.isequal(a::Cookie, b::Cookie) = a.name == b.name && a.path == b.path && a.domain == b.domain
Base.hash(x::Cookie, h::UInt) = hash(x.name, hash(x.path, hash(x.domain, h)))

==(a::Cookie,b::Cookie) = (a.name     == b.name)    &&
                          (a.value    == b.value)   &&
                          (a.path     == b.path)    &&
                          (a.domain   == b.domain)  &&
                          (a.expires  == b.expires) &&
                          (a.maxage   == b.maxage)  &&
                          (a.secure   == b.secure)  &&
                          (a.httponly == b.httponly)

# request cookie stringify-ing
function Base.String(c::Cookie, isrequest::Bool=true)
    io = IOBuffer()
    nm = strip(c.name)
    !iscookienamevalid(nm) && return ""
    write(io, sanitizeCookieName(nm), '=', sanitizeCookieValue(c.value))
    if !isrequest
        length(c.path) > 0 && write(io, "; Path=", sanitizeCookiePath(c.path))
        length(c.domain) > 0 && validCookieDomain(c.domain) && write(io, "; Domain=", c.domain[1] == '.' ? c.domain[2:end] : c.domain)
        validCookieExpires(c.expires) && write(io, "; Expires=", Dates.format(c.expires, Dates.RFC1123Format), " GMT")
        c.maxage > 0 && write(io, "; Max-Age=", string(c.maxage))
        c.maxage < 0 && write(io, "; Max-Age=0")
        c.httponly && write(io, "; HttpOnly")
        c.secure && write(io, "; Secure")
    end
    return String(take!(io))
end

function Base.string(cookiestring::String, cookies::Vector{Cookie}, isrequest::Bool=true)
    io = IOBuffer()
    !isempty(cookiestring) && write(io, cookiestring, cookiestring[end] == ';' ? "" : ";")
    len = length(cookies)
    for (i, cookie) in enumerate(cookies)
        write(io, String(cookie, isrequest), ifelse(i == len, "", "; "))
    end
    return String(take!(io))
end

validcookiepathbyte(b) = (' ' <= b < '\x7f') && b != ';'
validcookievaluebyte(b) = (' ' <= b < '\x7f') && b != '"' && b != ';' && b != '\\'

function parsecookievalue(raw, allowdoublequote::Bool)
    if allowdoublequote && length(raw) > 1 && raw[1] == '"' && raw[end] == '"'
        raw = raw[2:end-1]
    end
    for i = 1:length(raw)
        !validcookievaluebyte(raw[i]) && return "", false
    end
    return raw, true
end

iscookienamevalid(raw) = raw == "" ? false : any(isurlchar, raw)

const AlternateRFC1123Format = Dates.DateFormat("e, dd-uuu-yyyy HH:MM:SS G\\MT")

# readSetCookies parses all "Set-Cookie" values from
# the header h and returns the successfully parsed Cookies.
readsetcookies(host, cookies) = Cookie[readsetcookie(host, c) for c in cookies]

function readsetcookie(host, cookie)
    parts = split(strip(cookie), ';')
    length(parts) == 1 && parts[1] == "" && return Cookie()
    parts[1] = strip(parts[1])
    j = findfirst(parts[1], '=')
    j < 1 && return Cookie()
    name, value = parts[1][1:j-1], parts[1][j+1:end]
    iscookienamevalid(name) || return Cookie()
    value, ok = parsecookievalue(value, true)
    ok || return Cookie()
    c = Cookie(name, value)
    for x = 2:length(parts)
        parts[x] = strip(parts[x])
        length(parts[x]) == 0 && continue
        attr, val = parts[x], ""
        j = findfirst(parts[x], '=')
        if j > 0
            attr, val = attr[1:j-1], attr[j+1:end]
        end
        lowerattr = lowercase(attr)
        val, ok = parsecookievalue(val, false)
        if !ok
            push!(c.unparsed, parts[x])
            continue
        end
        if lowerattr == "secure"
            c.secure = true
        elseif lowerattr == "httponly"
            c.httponly = true
        elseif lowerattr == "domain"
            c.domain = val
        elseif lowerattr == "max-age"
            secs = tryparse(Int, val)
            (isnull(secs) || val[1] == '0') && continue
            c.maxage = max(Base.get(secs), -1)
        elseif lowerattr == "expires"
            try
                c.expires = DateTime(val, Dates.RFC1123Format)
            catch
                try
                    c.expires = DateTime(val, AlternateRFC1123Format)
                catch
                    continue
                end
            end
        elseif lowerattr == "path"
            c.path = val
        else
            push!(c.unparsed, parts[x])
        end
    end
    c.domain, c.hostonly = domainandtype(host == "" ? c.domain : host, c.domain)
    return c
end

# shouldsend determines whether e's cookie qualifies to be included in a
# request to host/path. It is the caller's responsibility to check if the
# cookie is expired.
function shouldsend(cookie::Cookie, https::Bool, host, path)
    return domainmatch(cookie, host) && pathmatch(cookie, path) && (https || !cookie.secure)
end

# domainMatch implements "domain-match" of RFC 6265 section 5.1.3.
function domainmatch(cookie::Cookie, host)
    cookie.domain == host && return true
    return !cookie.hostonly && hasdotsuffix(host, cookie.domain)
end

# hasdotsuffix reports whether s ends in "."+suffix.
function hasdotsuffix(s, suffix)
    return length(s) > length(suffix) && s[length(s)-length(suffix)] == '.' && s[length(s)-length(suffix)+1:end] == suffix
end

# pathMatch implements "path-match" according to RFC 6265 section 5.1.4.
function pathmatch(cookie::Cookie, requestpath)
    requestpath == cookie.path && return true
    if startswith(requestpath, cookie.path)
        if cookie.path[end] == '/'
            return true # The "/any/" matches "/any/path" case.
        elseif length(requestpath) >= length(cookie.path) + 1 && requestpath[length(cookie.path)+1] == '/'
            return true # The "/any" matches "/any/path" case.
        end
    end
    return false
end

function isIP(host)
    try
        Base.parse(IPAddr, host)
        return true
    catch e
        isa(e, ArgumentError) && return false
        rethrow(e)
    end
end

# domainAndType determines the cookie's domain and hostOnly attribute.
function domainandtype(host, domain)
    if domain == ""
        # No domain attribute in the SetCookie header indicates a
        # host cookie.
        return host, true
    end

    if isIP(host)
        # According to RFC 6265 domain-matching includes not being
        # an IP address.
        # TODO: This might be relaxed as in common browsers.
        return "", false
    end

    # From here on: If the cookie is valid, it is a domain cookie (with
    # the one exception of a public suffix below).
    # See RFC 6265 section 5.2.3.
    if domain[1] == '.'
        domain = domain[2:end]
    end

    if length(domain) == 0 || domain[1] == '.'
        # Received either "Domain=." or "Domain=..some.thing",
        # both are illegal.
        return "", false
    end
    domain = lowercase(domain)

    if domain[end] == '.'
        # We received stuff like "Domain=www.example.com.".
        # Browsers do handle such stuff (actually differently) but
        # RFC 6265 seems to be clear here (e.g. section 4.1.2.3) in
        # requiring a reject.  4.1.2.3 is not normative, but
        # "Domain Matching" (5.1.3) and "Canonicalized Host Names"
        # (5.1.2) are.
        return "", false
    end

    #TODO:
    # See RFC 6265 section 5.3 #5.
    # if j.psList != nil
    #     if ps := j.psList.PublicSuffix(domain); ps != "" && !hasDotSuffix(domain, ps)
    #         if host == domain
    #             # This is the one exception in which a cookie
    #             # with a domain attribute is a host cookie.
    #             return host, true, nil
    #         end
    #         return "", false
    #     end
    # end

    # The domain must domain-match host: www.mycompany.com cannot
    # set cookies for .ourcompetitors.com.
    if host != domain && !hasdotsuffix(host, domain)
        return "", false
    end

    return domain, false
end

# readCookies parses all "Cookie" values from the header h and
# returns the successfully parsed Cookies.
# if filter isn't empty, only cookies of that name are returned
function readcookies(h::Dict{String,String}, filter::String)
    if !haskey(h, "Cookie") && !haskey(h, "cookie")
        return Cookie[]
    end
    lines = Base.get(h, "Cookie", Base.get(h, "cookie", ""))

    cookies = Cookie[]
    for part in split(lines, ';')
        part = strip(part)
        length(part) <= 1 && continue
        j = findfirst(part, '=')
        if j >= 0
            name, val = part[1:j-1], part[j+1:end]
        else
            name, val = part, ""
        end
        !iscookienamevalid(name) && continue
        filter != "" && filter != name && continue
        val, ok = parsecookievalue(val, true)
        !ok && continue
        push!(cookies, Cookie(name, val))
    end
    return cookies
end

# validCookieExpires returns whether v is a valid cookie expires-value.
function validCookieExpires(dt)
    # IETF RFC 6265 Section 5.1.1.5, the year must not be less than 1601
    return Dates.year(dt) >= 1601
end

# validCookieDomain returns whether v is a valid cookie domain-value.
function validCookieDomain(v::String)
    isCookieDomainName(v) && return true
    isIP(v) && !contains(v, ":") && return true
    return false
end

# isCookieDomainName returns whether s is a valid domain name or a valid
# domain name with a leading dot '.'.  It is almost a direct copy of
# package net's isDomainName.
function isCookieDomainName(s::String)
    length(s) == 0 && return false
    length(s) > 255 && return false
    s = s[1] == '.' ? s[2:end] : s
    last = '.'
    ok = false
    partlen = 0
    for c in s
        if 'a' <= c <= 'z' || 'A' <= c <= 'Z'
            ok = true
            partlen += 1
        elseif '0' <= c <= '9'
            partlen += 1
        elseif c == '-'
            last == '.' && return false
            partlen += 1
        elseif c == '.'
            (last == '.' || last == '-') && return false
            (partlen > 63 || partlen == 0) && return false
            partlen = 0
        else
            return false
        end
        last = c
    end
    (last == '-' || partlen > 63) && return false
    return ok
end

sanitizeCookieName(n::String) = replace(replace(n, '\n', '-'), '\r', '-')
sanitizeCookieName(n) = sanitizeCookieName(String(n))

# http:#tools.ietf.org/html/rfc6265#section-4.1.1
# cookie-value      = *cookie-octet / ( DQUOTE *cookie-octet DQUOTE )
# cookie-octet      = %x21 / %x23-2B / %x2D-3A / %x3C-5B / %x5D-7E
#           ; US-ASCII characters excluding CTLs,
#           ; whitespace DQUOTE, comma, semicolon,
#           ; and backslash
# We loosen this as spaces and commas are common in cookie values
# but we produce a quoted cookie-value in when value starts or ends
# with a comma or space.
# See https:#golang.org/issue/7243 for the discussion.
function sanitizeCookieValue(v::String)
    v = String(filter(validcookievaluebyte, [Char(b) for b in Vector{UInt8}(v)]))
    length(v) == 0 && return v
    if v[1] == ' ' || v[1] == ',' || v[end] == ' ' || v[end] == ','
        return string('"', v, '"')
    end
    return v
end

sanitizeCookiePath(v) = filter(validcookiepathbyte, v)

end # module


text:
module DashCarbonComponents
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.1.4"

include("button.jl")
include("card.jl")
include("checkbox.jl")
include("column.jl")
include("datepicker.jl")
include("datepickerrange.jl")
include("dropdown.jl")
include("grid.jl")
include("link.jl")
include("multiselect.jl")
include("numberinput.jl")
include("radiobuttongroup.jl")
include("row.jl")
include("slider.jl")
include("tab.jl")
include("tabs.jl")
include("uishell.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_carbon_components",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "dash_carbon_components.min.js",
    external_url = nothing,
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_carbon_components.min.js.map",
    external_url = nothing,
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

ext: jl

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
# Set plot globals
ENV["PLOTS_TEST"] = "true"
ENV["GKSwstype"] = "100"

using Plots
gr()
default(show=false)

# This function show mirror that in ControlExamplePlots.jl/genplots.jl
# to make sure that the plots in these tests can be tested for accuracy
"""funcs, names = getexamples()
Get the example functions and names
"""
function getexamples()
    tf1 = tf([1],[1,1])
    tf2 = tf([1/5,2],[1,1,1])
    sys = [tf1 tf2]
    sysss = ss([-1 2; 0 1], [1 0; 1 1], [1 0; 0 1], [0.1 0; 0 -0.2])

    ws = 10.0 .^range(-2,stop=2,length=200)
    ts = 0:0.01:5
    bodegen() = begin
      setPlotScale("dB")
      bodeplot(sys,ws)
    end
    nyquistgen() = nyquistplot(sysss,ws)
    sigmagen() = sigmaplot(sysss,ws)
    #Only siso for now
    nicholsgen() = nicholsplot(tf1,ws)

    stepgen() = stepplot(sys, ts[end], l=(:dash, 4))
    impulsegen() = impulseplot(sys, ts[end], l=:blue)
    L = lqr(sysss.A, sysss.B, [1 0; 0 1], [1 0; 0 1])
    lsimgen() = lsimplot(sysss, (x,i)->-L*x, ts, [1;2])

    margingen() = marginplot([tf1, tf2], ws)
    gangoffourgen() = begin
      setPlotScale("log10");
      gangoffourplot(tf1, [tf(1), tf(5)])
    end
    pzmapgen() = pzmap(tf2, xlims=(-15,5))
    rlocusgen() = rlocusplot(tf2)

    refs = ["bode.png", "nyquist.png", "sigma.png", "nichols.png", "step.png",
            "impulse.png", "lsim.png", "margin.png", "gangoffour.png", "pzmap.png", "rlocus.png"]
    funcs = [bodegen, nyquistgen, sigmagen, nicholsgen, stepgen,
             impulsegen, lsimgen, margingen, gangoffourgen, pzmapgen, rlocusgen]

    funcs, refs
end

@testset "test_plots" begin
  funcs, names = getexamples()

  for i in eachindex(funcs)
    println("run_plot_$(i): $(names[i])")
    @test funcs[i]() isa Plots.Plot
  end

end


text:
import Distributions.MvNormal
import Distributions.Categorical
import Distributions.Multinomial

mvn = MvNormal([1.2,2.3], [1.4 1.0; 1.0 1.2])

@test isapprox(pdf(mvn, [1.0, 1.0]),  0.04794537749882221)

println("There should NOT be a warning between here")
r = rand(MersenneTwister(12), mvn)
println("and here.")
# @test isa(r, sampletype(mvn))

println("There should be a warning below.")
mn = Multinomial(4, 8)
r = rand(MersenneTwister(12), mn)
# @test isa(r, sampletype(mn))

println("There should NOT be a warning between here")
cat = Categorical(3)
r = rand(MersenneTwister(12), cat)
println("and here.")
# @test isa(r, sampletype(cat))


text:
# using Distributions, PyPlot
push!(LOAD_PATH, "../../VICAS/src/")
push!(LOAD_PATH, "../../DeepCorrection")
push!(LOAD_PATH, "../../DeepCorrection/src")
push!(LOAD_PATH, "../TrajectoryViz")

include("../TrajectoryViz/trajViz.jl")

th = linspace(0, 2*pi, 50)

conflictRange = SENSING_RANGE * 1.5

function check_NMAC(i, acArr, finished)
	ac = acArr[i]
	for j in 1:length(acArr)
		if j != i && finished[i] == false && finished[j] == false
			d = sqrt((ac.x - acArr[j].x)^2 + (ac.y - acArr[j].y)^2)
			if d <= NMAC_RANGE
				acArr[i].dmin = d
				acArr[j].dmin = d
				return true, d
			end
		end
	end
	return false, Inf
end

function check_colli(ac, acArr)
	for existing_ac in acArr
		if sqrt((ac.x - existing_ac.x)^2 + (ac.y - existing_ac.y)^2) <= NMAC_RANGE
			return true
		end
	end
end

function create_ac(spawn_range, CASType, policy)
	# Annulus type, centralized encounter
	spawn_range = rand(2000:4000)
	start_angle = rand(th)
	x0 = spawn_range * cos(start_angle)
	y0 = spawn_range * sin(start_angle)
	# conflictPoint_x = conflictRange * (- 2 * rand() + 1)
	# conflictPoint_y  = sqrt(conflictRange^2 - conflictPoint_x^2) * (- 2 * rand() + 1)
	# dest_angle = atan2(conflictPoint_y - y0, conflictPoint_x - x0)
	# dest_dist = (4000 - 2 * spawn_range) * rand() + 2 * spawn_range
	dest_dist = rand(2000:4000)
	# x_dest = x0 + cos(dest_angle) * dest_dist
	# y_dest = y0 + sin(dest_angle) * dest_dist
	x_dest = dest_dist * cos(start_angle + pi)
	y_dest = dest_dist * sin(start_angle + pi)

	# Circular type encounter
	# x0 = spawn_range * (- 2 * rand() + 1)
	# y0max  = sqrt(spawn_range^2 - x0^2)
	# y0     = y0max * (- 2 * rand() + 1)
	# dest_angle = rand(th)
	# dest_dist = 2*spawn_range
	# x_dest = cos(dest_angle) * dest_dist
	# y_dest = sin(dest_angle) * dest_dist

	# dest_angle = rand(th)
	# x_dest = spawn_range * cos(dest_angle)
	# y_dest = spawn_range * sin(dest_angle)
	ac = initOwnship(x0, y0, x_dest, y_dest, CASType, policy)
	return ac
end

function resRatio(;
				CASType::Symbol=:correctedSector,
				policyFile::AbstractString="",
				numAC::Int64=2,
				numSim::Int64=10000,
				timeout::Int64=200,
				correction_weight::Float64=0.09)

	problem_file = "../../DeepCorrection/logs/" * policyFile * "/final_problem.jld"
	weights_file = "../../DeepCorrection/logs/" * policyFile * "/final_weights.jld"
	policy = restore(problem_file=problem_file, weights_file=weights_file)
	policy.env.correction_weight = correction_weight
	spawn_range = SENSING_RANGE * 1.5
	resolvedStats = 0
	unresolved_NMAC_Stats = Int64[]
	unresolved_Timeout_Stats = 0
	resTimeArr = Vector{Float64}()
	avgResTime = 0
	criticality = Float64[]
	v_closing_v_cruise = Float64[]

	for sim in 1 : numSim
		# Initiate
		t = 1
		NMACFlag = false
		fail = false
		finished = falses(numAC)
		acArr = Ownship[]
		for i = 1 : numAC
			ac = create_ac(spawn_range, CASType, policy)
			while check_colli(ac, acArr) == true
				ac = create_ac(spawn_range, CASType, policy)
			end
			push!(acArr, ac)
		end

		push!(unresolved_NMAC_Stats, 0)

		while !all(finished) && NMACFlag == false
			for i in 1 : numAC
				if !finished[i]
					ownship = acArr[i]
					# Get intruder coords for ownship
					intCoords = Vector{Float64}()
					for j in 1 : length(acArr)
						if j != i && !finished[j] # not ownship and not finished
							intCoords = vcat(intCoords, [acArr[j].x, acArr[j].y], acArr[j].heading)
						end
					end
					numInt = Int(length(intCoords) / 3)

					if numInt > 0 && CASType != :NOCAS
						# Get state
						numObservedInt = getState!(ownship, numInt, intCoords)
						if numObservedInt > 0
							# Get action and advisory
							getAdvisory!(ownship)
						else
							ownship.advisory = :COC
						end
					else
						ownship.advisory = :COC
					end
					# print(i, " : ", ownship.state, " : ", ownship.advisory, "\n")
				end # if !finished[i]
			end # for

			for i in 1 : numAC
				if !finished[i]
					dest_angle = norm_angle(atan2(acArr[i].y_dest - acArr[i].y, acArr[i].x_dest - acArr[i].x))
					heading_error = norm_angle(dest_angle - acArr[i].heading)
					push!(v_closing_v_cruise, cos(heading_error))
					updateOwnshipCoords!(acArr[i])
				end

				if sqrt((acArr[i].x - acArr[i].x_dest)^2 + (acArr[i].y - acArr[i].y_dest)^2) < 80
					finished[i] = true
				end
			end # for

			for i in 1 : numAC
				nmac, d = check_NMAC(i, acArr, finished)
				if finished[i] == false && nmac == true # FAIL: NMAC
					push!(criticality, max(0.0, 1 - acArr[i].dmin / NMAC_RANGE))
					fail = true
					NMACFlag = true
					unresolved_NMAC_Stats[sim] += 1
				end
			end

			if all(finished) && NMACFlag == false # SUCCESS
				resolvedStats += 1
			end

			t += 1
			if t > timeout && NMACFlag == false # FAIL: timeout
				fail = true
				unresolved_Timeout_Stats += 1
				break # while
			end

		end #while

		if t <= timeout && NMACFlag == false # only count SUCCESS cases
			avgResTime += t
		end
	end # for sim

	return resolvedStats, unresolved_NMAC_Stats ./ 2, unresolved_Timeout_Stats, avgResTime / resolvedStats, criticality, v_closing_v_cruise
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
"""
	$(SIGNATURES)

Object that holds result from traversing an object and applying a function to each property.

Contains the name of the object (typically the field name in the parent object), its `DataType`, the returned function value. Also contains a `Vector{ApplyFctResult}` that holds the same information for each child object.
"""
mutable struct ApplyFctResult
    name :: Symbol
    type :: DataType
    fVal
    children :: Union{Nothing, Vector{ApplyFctResult}}
end

"""
	$(SIGNATURES)

Name of object.
"""
obj_name(x :: ApplyFctResult) = x.name;

"""
	$(SIGNATURES)

Type of object.
"""
obj_type(x :: ApplyFctResult) = x.type;

"""
	$(SIGNATURES)

Function value or `nothing`.
"""
fct_value(x :: ApplyFctResult) = x.fVal;

"""
	$(SIGNATURES)

Child objects, excluding those that contain no function values (themselves, or their children).
"""
children(x :: ApplyFctResult) = x.children;

function child_names(x :: ApplyFctResult)
    if isnothing(children(x))
        return nothing;
    else
        return [obj_name(child)  for child in children(x)]
    end
end

Base.show(io :: IO, a :: ApplyFctResult) =
    print(io, "ApplyFctResult $(obj_name(a)) with children $(child_names(a))");

"""
	$(SIGNATURES)

Function returns false when parsing object tree should stop at this object. Predefined as `false` for common "built-in" types, such as `Number`, `Symbol`, `Array`.
"""
continue_parsing(x) = true;
continue_parsing(::Number) = false;
continue_parsing(::Symbol) = false;
continue_parsing(::AbstractDict) = false;
continue_parsing(::AbstractString) = false;
continue_parsing(::AbstractArray) = false;
continue_parsing(::DataType) = false;

"""
	$(SIGNATURES)

Apply function `f` to all properties of an object and its children.
Stop parsing child objects if `continue_parsing(x) == false`.
Returns an `ApplyFctResult` object.

Beware: Some functions, such as `sum`, have fallback methods that apply to `Any`. That leads to crashes.

Argument `maxLevels` limits the number of recursions and avoids stack overflow if something goes wrong.
"""
function apply_fct_to_object(x, xName :: Symbol, f :: Function;
    maxLevels = 5)
    if continue_parsing(x)
        childV = apply_fct_to_children(x, f; maxLevels = (maxLevels - 1));
    else
        childV = nothing;
    end
    if hasmethod(f, (typeof(x), ))
        fValue = f(x);
    else
        fValue = nothing;
    end
    return ApplyFctResult(xName, typeof(x), fValue, childV)
end

# Returns a vector of entries, skipping those with `f(x) == nothing` and no children with entries.
# Returns `nothing` for an object without useful info.
function apply_fct_to_children(x, f; maxLevels = 5)
    if maxLevels < 1
        return nothing
    end
    nameV = propertynames(x);
    if isempty(nameV)
        return nothing
    end
    resultV = Vector{ApplyFctResult}();
    for cName in nameV
        child = getproperty(x, cName);
        if continue_parsing(child)
            grandChildV = apply_fct_to_children(child, f;
                maxLevels = (maxLevels - 1));
        else
            grandChildV = nothing;
        end
        if hasmethod(f, (typeof(child), ))
            fValue = f(child);
        else
            fValue = nothing;
        end
        if !isnothing(fValue) || !isnothing(grandChildV)
            push!(resultV, ApplyFctResult(cName, typeof(child), fValue, grandChildV));
        end
    end
    if isempty(resultV)
        return nothing
    else
        return resultV
    end
end

## ------------  Display structure of an object

"""
	$(SIGNATURES)

Returns a description for an object. `Nothing` is the default.

Define this for user defined types to return a `Vector{String}` or `Matrix{String}` that can be formatted as a table. [`desribe_object`](@ref) also shows descriptions of child objects.
"""
describe(x) = nothing;

"""
	$(SIGNATURES)

Write object description, including child objects, to IO.
"""
function show_description(io :: IO, x)
    [println(io, line)  for line in describe_object(x)];
    return nothing
end

show_description(x) = show_description(stdout, x);

"""
	$(SIGNATURES)

Describe an object and its children. Recursively calls `describe` and returns a vector of strings.

The default return value for `describe` is nothing. Then the object is not displayed.
Otherwise, the return is passed through `format_describe` for formatting. The result is an indented list which visualizes the structure of parent-child relationships as well.
"""
function describe_object(x)
    xResult = apply_fct_to_object(x, :Top, describe);
    return describe_result(xResult);
end

function describe_result(x :: ApplyFctResult; level = 0)
    outV = Vector{String}();
    indentStr = repeat("  ", level);
    push!(outV, "$indentStr$(obj_name(x))");

    # Own description or function values
    fVals = fct_value(x);
    if !isnothing(fVals)
        fmtStr = format_describe(fVals, level + 1);
        if isa(fmtStr, AbstractVector)
            for fVal in fmtStr
                push!(outV, fVal);
            end
        else
            push!(outV, fmtStr);
        end
    end
    if !isnothing(children(x))
        for child in children(x)
            append!(outV, describe_result(child; level = level + 1));
        end
    end
    return outV
end

indent_string(level) = repeat("  ", level);

format_describe(x, level :: Integer) = indent_string(level) * "$x";
format_describe(x :: AbstractString, level :: Integer) = indent_string(level) * x;
format_describe(x :: AbstractVector{S1}, level :: Integer) where S1 <: AbstractString =
    [format_describe(s, level)  for s in x];

format_describe(x :: AbstractMatrix{S1}, level :: Integer) where S1 <: AbstractString =
    [format_describe(x[ir,1] * ":  " * x[ir, 2], level) for ir = 1 : size(x, 1)];

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


text:
using HDF5
using Test

function learn_open_hdf5()
    fname = tempname(pwd())
    @test !isfile(fname)

    fid = h5open(fname, "cw")
    @test isfile(fname)

    close(fid)
    @test isfile(fname)

    rm(fname)
    @test !isfile(fname)

    # The "do" syntax
    h5open(fname, "cw") do fid
        @test isfile(fname)
    end
    rm(fname)
    @test !isfile(fname)
end

function learn_writing_a_group_or_dataset()
    # Test const data
    data = [13.0, 17.0, 19.0]

    ## Create aHDF5 file
    # Must close at the end
    # One can use the do syntax
    #   h5open(fname, "cw") do fid
    #   end
    fname = tempname(pwd())
    fid = h5open(fname, "cw")

    ## Create group: /mygroup/
    # Group is like a folder
    g = create_group(fid, "mygroup")       # Create
    # Access non-existing group
    @test_throws KeyError fid["thegroup"]
    # Access a group
    x = fid["mygroup"]                     # Get
    @test g    != x
    @test g.id != x.id
    @test g.file == x.file

    ## Create dataset 1
    # From the source code
    # Create datasets and attributes with "native" types, but don't write the data.
    # The return syntax is: dset, dtype = create_dataset(parent, name, data; properties...)
    name = "mydata1"
    xd, dtype = create_dataset(fid, name, data)
    @test xd isa HDF5.Dataset
    @test dtype isa HDF5.Datatype #@show dtype # HDF5.Datatype: H5T_IEEE_F64LE
    x = read(fid, name)             # Read
    @test x isa Vector{Float64}
    @test x == [0.0, 0.0, 0.0]
    @test x != [13.0, 17.0, 19.0]   # x is not [13.0, 17.0, 19.0]
    x[1] = 7.0                      # Modify x does not modify the file
    x = read(fid, name)             # Read
    @test x ==  [0.0, 0.0, 0.0]     #
    xd[1] = 7.0                     # Modify dset modifies the file
    x = read(fid, name)             # Read
    @test x ==  [7.0, 0.0, 0.0]     #

    ## Create dataset 2
    # Create by specifing data type and size.
    name = "mydata2"
    xd = create_dataset(fid, name, datatype(eltype(data)), dataspace(size(data)))
    @test xd isa HDF5.Dataset
    x = read(fid, name)
    @test x isa Vector{Float64}
    @test x == [0.0, 0.0, 0.0]      # Initialized to zeros
    x[1] = 7.0                      # Modify x does not modify the file
    x = read(fid, name)             # Read
    @test x ==  [0.0, 0.0, 0.0]
    xd[1] = 7.0                     # Modify dset modifies the file
    x = read(fid, name)
    @test x ==  [7.0, 0.0, 0.0]

    ## Create dataset 3
    # Create by [] operator
    name = "mydata3"
    fid[name] = data               # Create
    x = read(fid, name)            # Read
    @test x isa Vector{Float64}
    @test x == [13.0, 17.0, 19.0]  # x is set correctly
    x[1] = 7.0                     # This does not change the data
    x = read(fid, name)
    @test x == [13.0, 17.0, 19.0]
    d = open_dataset(fid, name)    # Get
    d[1] = 7.0                     # Write
    x = read(fid, name)            # Read
    @test x == [7.0, 17.0, 19.0]   # Data is set correctly

    ## Create dataset 4
    # Create by write()
    name = "mydata4"
    write(fid, name, data)
    x = read(fid, name)
    @test x == [13.0, 17.0, 19.0]  # x is set correctly
    d = open_dataset(fid, name)    # Get
    d[1] = 7.0                     # Write
    x = read(fid, name)            # Read
    @test x == [7.0, 17.0, 19.0]   # Data is set correctly

    close(fid)

    rm(fname)
end

@testset "Basic" begin
    learn_open_hdf5()
    learn_writing_a_group_or_dataset()
end

nothing


text:
using Makie, Statistics

# Point3f0 is a 3-tuple of 32-bit floats for 3-dimensional space, and all Points are 3D.
Point = Point3f0

# a Face is defined by the points that are its vertices, in order.
Face = Vector{Point}

# an Edge is a line segment where the points are sorted
struct Edge
    p1::Point
    p2::Point
    Edge(a, b) = new(min(a, b), max(a, b))
end

edgemidpoint(edge) = (edge.p1  + edge.p2) / 2.0
facesforpoint(p, faces) = [f for f in faces if p in f]
facesforedge(e, faces) = [f for f in faces if (e.p1 in f) && (e.p2 in f)]
nexttohole(edge, faces) = length(facesforedge(edge, faces)) < 2

function newedgepoint(edge, faces)
    f = facesforedge(edge, faces)
    p1, p2, len = edge.p1, edge.p2, length(f)
    if len == 2
        return (p1 + p2 + mean(f[1]) + mean(f[2])) / 4.0
    elseif len == 1
        return (p1 + p2 + mean(f[1])) / 3.0
    end
    return (p1 + p2) / 2.0
end

function edgesforface(face)
    ret, indices = Vector{Edge}(), collect(1:length(face))
    for i in 1:length(face)-1
        push!(ret, Edge(face[indices[1]], face[indices[2]]))
        indices .= circshift(indices, 1)
    end
    ret
end

function edgesforpoint(p, faces)
    f = filter(x -> p in x, faces)
    return filter(e -> p == e.p1 || p == e.p2, mapreduce(edgesforface, vcat, f))
end

function adjacentpoints(point, face)
    a = indexin([point], face)
    if a[1] != nothing
        adjacent = (a[1] == 1) ? [face[end], face[2]] :
            a[1] == length(face) ? [face[end-1], face[1]] :
            [face[a[1] - 1], face[a[1] + 1]]
        return sort(adjacent)
    else
        throw("point $point not in face $face")
    end
end

adjacentedges(point, face) = [Edge(point, x) for x in adjacentpoints(point, face)]
facewrapped(face) = begin f = deepcopy(face); push!(f, f[1]); f end
drawface(face, colr) = lines(facewrapped(face); color=colr)
drawface!(face, colr) = lines!(facewrapped(face); color=colr)
drawfaces!(faces, colr) = for f in faces drawface!(f, colr) end
const colors = [:red, :green, :blue, :gold]

function drawfaces(faces, colr)
    scene = drawface(faces[1], colr)
    if length(faces) > 1
        for f in faces[2:end]
            drawface!(f, colr)
        end
    end
    scene
end

function catmullclarkstep(faces)
    d, E, dprime = Set(reduce(vcat, faces)), Dict{Vector, Point}(), Dict{Point, Point}()
    for face in faces, (i, p) in enumerate(face)
        edge = (p == face[end]) ? Edge(p, face[1]) : Edge(p, face[i + 1])
        E[[edge, face]] = newedgepoint(edge, faces)
    end
    for p in d
        F = mean([mean(face) for face in facesforpoint(p, faces)])
        pe = edgesforpoint(p, faces)
        R = mean(map(edgemidpoint, pe))
        n = length(pe)
        dprime[p] = (F + 2 * R + p * (n - 3)) / n
    end
    newfaces = Vector{Face}()
    for face in faces
        v = mean(face)
        for point in face
            fp1, fp2 = map(x -> E[[x, face]], adjacentedges(point, face))
            push!(newfaces, [fp1, dprime[point], fp2, v])
        end
    end
    return newfaces
end

"""
    catmullclark(faces, iters, scene)

Perform a multistep Catmull-Clark subdivision of a surface. See Wikipedia or page 53
of http://graphics.stanford.edu/courses/cs468-10-fall/LectureSlides/10_Subdivision.pdf
Plots the iterations, with colors for each iteration as set in the colors array.
Uses a Makie Scene of scene to plot the iters iterations.
"""
function catmullclark(faces, iters, scene)
    nextfaces = deepcopy(faces)
    for i in 1:iters
        nextfaces = catmullclarkstep(nextfaces)
        drawfaces!(nextfaces, colors[i])
        display(scene)
        sleep(1)
    end
end

const inputpoints = [
    [-1.0, -1.0, -1.0],
    [-1.0, -1.0, 1.0],
    [-1.0, 1.0, -1.0],
    [-1.0, 1.0, 1.0],
    [1.0, -1.0, -1.0],
    [1.0, -1.0, 1.0],
    [1.0, 1.0, -1.0],
    [1.0, 1.0, 1.0]
]

const inputfaces = [
    [0, 4, 5, 1],
    [4, 6, 7, 5],
    [6, 2, 3, 7],
    [2, 0, 1, 3],
    [1, 5, 7, 3],
    [0, 2, 6, 4]
]

const faces = [map(x -> Point3f0(inputpoints[x]), p .+ 1) for p in inputfaces]

scene = drawfaces(faces, :black)
display(scene)
sleep(1)

catmullclark(faces, 4, scene)

println("Press Enter to continue", readline())

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Keymap used to encode information regarding keyboard layout, and country and language codes.

A string representation can be obtained from a `Keymap` by using `String(keymap)`.
"""
mutable struct Keymap <: Handle
    h::Ptr{xkb_keymap}
    ctx::Ptr{xkb_context}
    state::Ptr{xkb_state}
    conn::Connection
    device_id
    function Keymap(h, ctx, state, conn, id)
        km = new(h, ctx, state, conn, id)

        # get notified of state changes
        @check xcb_xkb_select_events_aux(conn, id, XCB_XKB_EVENT_TYPE_STATE_NOTIFY, true, true, false, 0, Ref(event_details_xkb(Dict("state" => true))))

        finalizer(km) do x
            xkb_state_unref(x.state)
            xkb_keymap_unref(x.h)
            xkb_context_unref(x.ctx)
        end
    end
end

function event_details_xkb(fields::Dict{String, Bool})
    names = fieldnames(xcb_xkb_select_events_details_t)
    flags = zeros(length(names))
    fields_affect = Dict("affect" * uppercasefirst(k) => v for (k, v) ∈ fields)
    for i ∈ 1:2:length(names)
        access_field = string(names[i])
        if access_field ∈ keys(fields_affect)
            flags[i] = 1 # affect
            flags[i + 1] = fields_affect[access_field] # details
        end
    end
    xcb_xkb_select_events_details_t(flags...)
end

"""
Construct a `Keymap` using `conn` as the connection to the X Server.

It uses XKB, the X Keyboard extension, which must be initialized for each connection, typically when creating a keymap (it is unlikely to be used before then). If this is not your first call to this constructor with this connection, you should set `setup_xkb` to `false`.
"""
function Keymap(conn::Connection; setup_xkb=true)
    if setup_xkb
        ret_code = xkb_x11_setup_xkb_extension(conn, 1, 0, XKB_X11_SETUP_XKB_EXTENSION_NO_FLAGS, C_NULL, C_NULL, C_NULL, C_NULL)
        ret_code == 0 && error("XKB extension setup failed")
    end
    ctx = xkb_context_new(XKB_CONTEXT_NO_DEFAULT_INCLUDES)
    ctx == C_NULL && error("Context creation failed")
    device_id = xkb_x11_get_core_keyboard_device_id(conn)
    device_id == -1 && error("No keyboard device ID could be retrieved")
    keymap = xkb_x11_keymap_new_from_device(ctx, conn, device_id, XKB_KEYMAP_COMPILE_NO_FLAGS)
    state = xkb_x11_state_new_from_device(keymap, conn, device_id)
    state == C_NULL && error("State creation failed")
    Keymap(keymap, ctx, state, conn, device_id)
end

function Base.String(km::Keymap)
    km_name_ptr = xkb_keymap_get_as_string(km, XKB_KEYMAP_FORMAT_TEXT_V1)
    km_name_ptr == C_NULL && error("Failed to obtain a string from the keymap $km")
    unsafe_string(km_name_ptr)
end

"""
Return a `String` representation of an integer XCB keysym.
"""
function keysym_string(keysym::Integer; max_chars=50)
    char_tmp = zeros(UInt8, max_chars)
    char_tmp_ptr = pointer(char_tmp)
    GC.@preserve char_tmp begin
        val = xkb_keysym_get_name(keysym, pointer(char_tmp), sizeof(eltype(char_tmp)) * length(char_tmp))
        val == -1 && @error("Failed to obtain a keysym string for $(hex(keysym))")
        str = unsafe_string(char_tmp_ptr)
    end
    str
end

"""
Obtain the name of a physical key identified by its `keycode` using a keymap.
"""
function key_name(km::Keymap, keycode::Integer)
    ptr = xkb_keymap_key_get_name(km, keycode)
    ptr == C_NULL && error("Failed to obtain a name from the keycode $keycode")
    Symbol(unsafe_string(ptr))
end

"""
Obtain a keycode from a key `name` using a keymap.
"""
function keycode(km::Keymap, name::Symbol)
    keycode = xkb_keymap_key_by_name(km, string(name))
    keycode == typemax(UInt32) && error("Failed to obtain a keycode from the key name $name")
    keycode
end

function keystroke_info(km::Keymap, event_details::EventDetails)
    @unpack key_name, key, input = event_details.data
    keycode = xkb_keymap_key_by_name(km, string(key_name))
    "Key \e[31m$key_name\e[m (code \e[33m$keycode\e[m): input \"\e[36m$input\e[m\" from symbol \e[36m$key\e[m"
end

Base.Char(km::Keymap, keycode::Integer) = Char(xkb_state_key_get_utf32(km.state, keycode))

"""
Generate a `KeySymbol` from a keycode.
"""
KeySymbol(km::Keymap, keycode::Integer) = key_symbol_from_keysym_string(keysym_string(xkb_state_key_get_one_sym(km.state, keycode)))

"""
Generate a `KeySymbol` from a key name.
"""
KeySymbol(km::Keymap, key_name::Symbol) = KeySymbol(km, xkb_keymap_key_by_name(km, string(key_name)))

"""
Generate a `KeySymbol` from a XCB keysym string.
"""
function key_symbol_from_keysym_string(keysym::AbstractString)
    keysym = Symbol(lowercase(keysym))
    KeySymbol(standardize_key_symbol(keysym))
end

"""
Transform a key symbol obtained via XCB to a standard symbol defined in `WindowAbstractions`.
"""
standardize_key_symbol(key_symbol::Symbol) = key_symbol ∈ keys(keysym_names_translation) ? keysym_names_translation[key_symbol] : key_symbol

const keysym_names_translation = Dict(
    :shift_r          => :shift_right,
    :shift_l          => :shift_left,
    :control_r        => :control_right,
    :control_l        => :control_left,
    :iso_level3_shift => :alt_gr,
    :alt_l            => :alt_left,
    :return           => :enter,
)

"""
Produce a key event based on a key name, a modifier state and an action using a keymap.
"""
function key_event_from_name(km::Keymap, key_name::Symbol, modifiers::KeyModifierState, action::KeyAction)
    KeyEvent(key_name, KeySymbol(km, key_name), Char(km, keycode(km, key_name)), modifiers, action)
end


text:
"""
Implementation of dominated Coupling From The Past (dCFTP) developed by [KeMo99](@cite) and [KeMo00](@cite) for [`PRS.AbstractSpatialPointProcess`](@ref).

**See also**

- [Hub16](@cite)
- [Kendall's notes on perfect simulation](https://warwick.ac.uk/fac/sci/statistics/staff/academic-research/kendall/personal/ppt/428.pdf)
"""

@doc raw"""
    papangelou_conditional_intensity(pp::AbstractSpatialPointProcess, x, X)

Compute the [Papangelou conditional intensity](https://en.wikipedia.org/wiki/Point_process#Papangelou_intensity_function) of `pp`, as the ratio of densities ``\frac{f(X \cup x)}{f(X)}``.

**See also**

- Section 6.1.1 [MoWa04](@cite)
"""
function papangelou_conditional_intensity(pp::AbstractSpatialPointProcess, x, X) end

"""
    upper_bound_papangelou_conditional_intensity(pp::AbstractSpatialPointProcess)

Compute an upper bound on the [`PRS.papangelou_conditional_intensity`](@ref) of `pp`

**See also**

- Equation 2.1 [KeMo99](@cite)
"""
function upper_bound_papangelou_conditional_intensity(pp::AbstractSpatialPointProcess) end

"""
    isrepulsive(pp::AbstractSpatialPointProcess)

Property of a [`PRS.AbstractSpatialPointProcess`](@ref).
"""
function isrepulsive(pp::AbstractSpatialPointProcess) end

"""
    isattractive(pp::AbstractSpatialPointProcess)

Property of a [`PRS.AbstractSpatialPointProcess`](@ref).
"""
function isattractive(pp::AbstractSpatialPointProcess) end

"""
    generate_sample_dcftp(
        [rng::Random.AbstractRNG,]
        pp::AbstractSpatialPointProcess{T};
        win::Union{Nothing,AbstractSpatialWindow}=nothing,
        n₀::Int=1
    )::Vector{T} where {T}

Generate an exact sample from a spatial point process `pp` on window `win` using dominated coupling from the past.

- Default window (`win=nothing`) is `window(pp)=pp.window`
- Initial coalescence check performed after `n₀` steps.
- Seed or random number generator is passed via `rng`.

**See also**

- [KeMo99](@cite), [KeMo00](@cite)
- Section 11.2.6 [MoWa04](@cite)
"""
function generate_sample_dcftp(
    rng::Random.AbstractRNG,
    pp::AbstractSpatialPointProcess{T},
    win::Union{Nothing,AbstractSpatialWindow}=nothing,
    n₀::Int=1
)::Vector{T} where {T}
    @assert n₀ >= 1

    window_ = isnothing(win) ? window(pp) : win
    β = upper_bound_papangelou_conditional_intensity(pp)
    birth_rate = β * volume(window_)

    # Dominating process
    hppp = HomogeneousPoissonPointProcess(β, window_)
    D = Set{T}(eachcol(generate_sample(rng, hppp)))

    M = Float64[]  # Marking process
    R = T[]        # Recording process

    steps = -1:-1:-n₀
    while true
        backward_extend!(rng, D, M, R, steps, birth_rate, window_)
        coupling, L = forward_coupling(D, M, R, pp, β)
        coupling && return collect(L)
        steps = (steps.stop-1):-1:(2*steps.stop)
    end
end

function generate_sample_dcftp(pp::AbstractPointProcess, args...)
    return generate_sample_dcftp(Random.default_rng(), pp, args...)
end

"""
    backward_extend!(
        rng::Random.AbstractRNG,
        D::Set{T},          # Dominating process
        M::Vector{Float64}, # Marking process
        R::Vector{T},       # Recording process
        steps::StepRange,   # Number of backward steps
        birth_rate::Real,
        win::AbstractSpatialWindow
    ) where {T}

Sample from the dominating birth-death process backwards in time according to `steps`.
The marks `M` and the points `R` which were added (uniform mark) / deleted (mark=0) along the run are recorded (`pushfirst!`).
The final state of the dominating process is updated to `D`.

**See also**

- `SBDevolve()` [KeMo99](@cite)
"""
function backward_extend!(
    rng::Random.AbstractRNG,
    D::Set{T},          # Dominating process
    M::Vector{Float64}, # Marking process
    R::Vector{T},       # Recording process
    steps::StepRange,   # Number of backward steps
    birth_rate::Real,
    win::AbstractSpatialWindow
) where {T}
    for _ in steps
        card_D = length(D)
        if rand(rng) < card_D / (birth_rate + card_D)
            # forward death (delete) ≡ backward birth (pushfirst)
            x = rand(rng, D)
            delete!(D, x)
            pushfirst!(R, x)
            pushfirst!(M, rand(rng))
        else
            # forward birth (push) ≡ backward death (pushfirst)
            x = rand(rng, win)
            push!(D, x)
            pushfirst!(R, x)
            pushfirst!(M, 0.0)
        end
    end
end

"""
    forward_coupling(
        D::Set{T},           # Dominating process
        M::Vector{Float64},  # Marking process
        R::Vector{T},        # Recording process
        pp::AbstractSpatialPointProcess{T},
        β::Real              # Upper bound on papangelou conditional intensity
    ) where {T}

Check if coalescence occured between the lower and upper bounding processes and return the state of the lower bounding process at time 0.

**See also**

- `SBDadd()` [KeMo99](@cite)
"""
function forward_coupling(
    D::Set{T},           # Dominating process
    M::Vector{Float64},  # Marking process
    R::Vector{T},        # Recording process
    pp::AbstractSpatialPointProcess{T},
    β::Real              # Upper bound on papangelou conditional intensity
) where {T}
    # L ⊆ X ⊆ U ⊆ D, where X is the target process
    L, U = empty(D), copy(D)
    for (m, x) in zip(M, R)
        if m > 0  # if birth occured in D
            if isrepulsive(pp)
                if m < papangelou_conditional_intensity(pp, x, U) / β
                    push!(L, x)
                    push!(U, x)
                elseif m < papangelou_conditional_intensity(pp, x, L) / β
                    push!(U, x)
                end
            elseif isattractive(pp)
                if m < papangelou_conditional_intensity(pp, x, L) / β
                    push!(L, x)
                    push!(U, x)
                elseif m < papangelou_conditional_intensity(pp, x, U) / β
                    push!(U, x)
                end
            else
                error("The current implementation of Dominated Coupling From The Past "
                       * "requires the target point process to be attractive or repulsive.")
            end
        else  # if death occured in D
            delete!(L, x)
            delete!(U, x)
        end
    end
    # Check coalescence L = U (Since L ⊆ U: L = U ⟺ |L| = |U|)
    return length(L) == length(U), L
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineProperties <: SwaggerModel
    hardwareProfile::Any # spec type: Union{ Nothing, HardwareProfile } # spec name: hardwareProfile
    storageProfile::Any # spec type: Union{ Nothing, StorageProfile } # spec name: storageProfile
    osProfile::Any # spec type: Union{ Nothing, OSProfile } # spec name: osProfile
    networkProfile::Any # spec type: Union{ Nothing, NetworkProfile } # spec name: networkProfile
    diagnosticsProfile::Any # spec type: Union{ Nothing, DiagnosticsProfile } # spec name: diagnosticsProfile
    availabilitySet::Any # spec type: Union{ Nothing, SubResource } # spec name: availabilitySet
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    instanceView::Any # spec type: Union{ Nothing, VirtualMachineInstanceView } # spec name: instanceView
    licenseType::Any # spec type: Union{ Nothing, String } # spec name: licenseType
    vmId::Any # spec type: Union{ Nothing, String } # spec name: vmId

    function VirtualMachineProperties(;hardwareProfile=nothing, storageProfile=nothing, osProfile=nothing, networkProfile=nothing, diagnosticsProfile=nothing, availabilitySet=nothing, provisioningState=nothing, instanceView=nothing, licenseType=nothing, vmId=nothing)
        o = new()
        validate_property(VirtualMachineProperties, Symbol("hardwareProfile"), hardwareProfile)
        setfield!(o, Symbol("hardwareProfile"), hardwareProfile)
        validate_property(VirtualMachineProperties, Symbol("storageProfile"), storageProfile)
        setfield!(o, Symbol("storageProfile"), storageProfile)
        validate_property(VirtualMachineProperties, Symbol("osProfile"), osProfile)
        setfield!(o, Symbol("osProfile"), osProfile)
        validate_property(VirtualMachineProperties, Symbol("networkProfile"), networkProfile)
        setfield!(o, Symbol("networkProfile"), networkProfile)
        validate_property(VirtualMachineProperties, Symbol("diagnosticsProfile"), diagnosticsProfile)
        setfield!(o, Symbol("diagnosticsProfile"), diagnosticsProfile)
        validate_property(VirtualMachineProperties, Symbol("availabilitySet"), availabilitySet)
        setfield!(o, Symbol("availabilitySet"), availabilitySet)
        validate_property(VirtualMachineProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(VirtualMachineProperties, Symbol("instanceView"), instanceView)
        setfield!(o, Symbol("instanceView"), instanceView)
        validate_property(VirtualMachineProperties, Symbol("licenseType"), licenseType)
        setfield!(o, Symbol("licenseType"), licenseType)
        validate_property(VirtualMachineProperties, Symbol("vmId"), vmId)
        setfield!(o, Symbol("vmId"), vmId)
        o
    end
end # type VirtualMachineProperties

const _property_map_VirtualMachineProperties = Dict{Symbol,Symbol}(Symbol("hardwareProfile")=>Symbol("hardwareProfile"), Symbol("storageProfile")=>Symbol("storageProfile"), Symbol("osProfile")=>Symbol("osProfile"), Symbol("networkProfile")=>Symbol("networkProfile"), Symbol("diagnosticsProfile")=>Symbol("diagnosticsProfile"), Symbol("availabilitySet")=>Symbol("availabilitySet"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("instanceView")=>Symbol("instanceView"), Symbol("licenseType")=>Symbol("licenseType"), Symbol("vmId")=>Symbol("vmId"))
const _property_types_VirtualMachineProperties = Dict{Symbol,String}(Symbol("hardwareProfile")=>"HardwareProfile", Symbol("storageProfile")=>"StorageProfile", Symbol("osProfile")=>"OSProfile", Symbol("networkProfile")=>"NetworkProfile", Symbol("diagnosticsProfile")=>"DiagnosticsProfile", Symbol("availabilitySet")=>"SubResource", Symbol("provisioningState")=>"String", Symbol("instanceView")=>"VirtualMachineInstanceView", Symbol("licenseType")=>"String", Symbol("vmId")=>"String")
Base.propertynames(::Type{ VirtualMachineProperties }) = collect(keys(_property_map_VirtualMachineProperties))
Swagger.property_type(::Type{ VirtualMachineProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineProperties }, property_name::Symbol) =  _property_map_VirtualMachineProperties[property_name]

function check_required(o::VirtualMachineProperties)
    true
end

function validate_property(::Type{ VirtualMachineProperties }, name::Symbol, val)
end


text:
# Autogenerated wrapper script for Qt6Base_jll for x86_64-apple-darwin
export libqt6concurrent, libqt6core, libqt6dbus, libqt6gui, libqt6network, libqt6opengl, libqt6printsupport, libqt6sql, libqt6test, libqt6widgets, libqt6xml

using Xorg_libXext_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using Xorg_libXrender_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
using Zlib_jll
using CompilerSupportLibraries_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Qt6Base")
JLLWrappers.@declare_library_product(libqt6concurrent, "@rpath/QtConcurrent")
JLLWrappers.@declare_library_product(libqt6core, "@rpath/QtCore")
JLLWrappers.@declare_library_product(libqt6dbus, "@rpath/QtDBus")
JLLWrappers.@declare_library_product(libqt6gui, "@rpath/QtGui")
JLLWrappers.@declare_library_product(libqt6network, "@rpath/QtNetwork")
JLLWrappers.@declare_library_product(libqt6opengl, "@rpath/QtOpenGL")
JLLWrappers.@declare_library_product(libqt6printsupport, "@rpath/QtPrintSupport")
JLLWrappers.@declare_library_product(libqt6sql, "@rpath/QtSql")
JLLWrappers.@declare_library_product(libqt6test, "@rpath/QtTest")
JLLWrappers.@declare_library_product(libqt6widgets, "@rpath/QtWidgets")
JLLWrappers.@declare_library_product(libqt6xml, "@rpath/QtXml")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libXext_jll, Xorg_libxcb_jll, Xorg_xcb_util_wm_jll, Xorg_xcb_util_image_jll, Xorg_xcb_util_keysyms_jll, Xorg_xcb_util_renderutil_jll, Xorg_libXrender_jll, xkbcommon_jll, Libglvnd_jll, Fontconfig_jll, Glib_jll, Zlib_jll, CompilerSupportLibraries_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libqt6concurrent,
        "lib/QtConcurrent.framework/Versions/A/QtConcurrent",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6core,
        "lib/QtCore.framework/Versions/A/QtCore",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6dbus,
        "lib/QtDBus.framework/Versions/A/QtDBus",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6gui,
        "lib/QtGui.framework/Versions/A/QtGui",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6network,
        "lib/QtNetwork.framework/Versions/A/QtNetwork",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6opengl,
        "lib/QtOpenGL.framework/Versions/A/QtOpenGL",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6printsupport,
        "lib/QtPrintSupport.framework/Versions/A/QtPrintSupport",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6sql,
        "lib/QtSql.framework/Versions/A/QtSql",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6test,
        "lib/QtTest.framework/Versions/A/QtTest",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6widgets,
        "lib/QtWidgets.framework/Versions/A/QtWidgets",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6xml,
        "lib/QtXml.framework/Versions/A/QtXml",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
$(TYPEDEF)

Holds the three parameters of a Lorentzian peak, defined as A / sqrt((ω-ω0)² + Γ²).

# Fields
$(TYPEDFIELDS)

"""
struct Lorentzian
    ω0::Float64
    Γ::Float64
    A::Float64
    Lorentzian(;ω0::Float64, Γ::Float64) = new(ω0, Γ, 1) # default peak height is 1
    Lorentzian(ω0, Γ, A) = new(ω0, Γ, A)
 end

 """
 $(TYPEDEF)

 Holds a set of `Lorentzian` objects belonging to a variable.

 # Fields
 $(TYPEDFIELDS)

 # Constructor
 ```julia
 JacobianSpectrum(res::Result; index::Int, branch::Int)
 ```

 """
 mutable struct JacobianSpectrum
    peaks::Vector{Lorentzian}
 end

"""
$(TYPEDEF)

Holds the compiled response matrix of a system.

# Fields
$(TYPEDFIELDS)

"""
struct ResponseMatrix
    """The response matrix (compiled)."""
    matrix::Matrix{Function}
    """Any symbolic variables in `matrix` to be substituted at evaluation."""
    symbols::Vector{Num}
    """The frequencies of the harmonic variables underlying `matrix`. These are needed to transform
    the harmonic variables to the non-rotating frame."""
    ωs::Vector{Num}

    ResponseMatrix(matrix, symbols, ωs) = new(matrix, symbols, ωs)
end


text:
# mostly unchanging functions that help configure TMI
"""
    function config_from_nc(TMIversion)
    Configure TMI environment from NetCDF input file.
# Arguments
- `TMIversion`: TMI version for water-mass/circulation model
# Output
- `A`: TMI steady-state water-mass matrix
- `Alu`: LU decomposition of A
- `γ`: TMI grid properties
- `TMIfile`: TMI file name
"""
function config_from_nc(TMIversion)

    #make datdir() if it doesn't exist
    !isdir(pkgdatadir()) && mkpath(pkgdatadir())
    TMIfile = pkgdatadir("TMI_"*TMIversion*".nc")

    #if TMIfile doesn't exist, get GDrive url and download
    if !isfile(TMIfile)

        # add a workaround for large files
        if TMIversion == "modern_180x90x33_GH11_GH12"
            println("workaround for 2° x 2°")
            shellscript = pkgsrcdir("read_nc_modern_180x90x33_GH11_GH12.sh")
            run(`sh $shellscript`)
            mv(joinpath(pwd(),"TMI_"*TMIversion*".nc"),TMIfile)
        else
            println("read via GoogleDrive.jl")
            #- `url`: Google Drive URL for data
            url = ncurl(TMIversion)
            google_download(url,pkgdatadir())
        end

    end

    # get properties of grid
    lon,lat,depth = gridprops(TMIfile)

    # read Cartesian Index from file.
    I = cartesianindex(TMIfile)

    # make a mask
    # first choice: smaller but inconsistent with input grid
    #wet = falses(maximum(I)[1],maximum(I)[2],maximum(I)[3])
    wet = falses(length(lon),length(lat),length(depth))
    wet[I] .= 1

    R = linearindex(wet)

    println("A")
    @time A = watermassmatrix(TMIfile)

    # LU factorization for efficient matrix inversions
    println("Alu")
    @time Alu = lu(A)

    γ = Grid(lon,lat,depth,I,R,wet)

    # would be good to make this optional
    println("L=")
    @time L = circulationmatrix(TMIfile,A,γ)

    println("B=")
    @time B = boundarymatrix(TMIfile,γ)

    return  A, Alu, γ, TMIfile, L, B

end

"""
Configure TMI environment from original MATLAB output
"""
function config_from_mat(TMIversion)

    #- `url`: Google Drive URL for data
    url = maturl(TMIversion)
    TMIfile = pkgdatadir("TMI_"*TMIversion*".mat")
    TMIfilegz = TMIfile*".gz"
    println(TMIfile)
    !isdir(pkgdatadir()) && mkpath(pkgdatadir())
    #    !isfile(TMIfilegz) & !isfile(TMIfile) ? google_download(url,pkgdatadir()) : nothing

    if TMIversion == "modern_180x90x33_GH11_GH12"
        println("workaround for 2° x 2°")
        shellscript = pkgsrcdir("read_mat_modern_180x90x33_GH11_GH12.sh")
        run(`sh $shellscript`)
        mv(joinpath(pwd(),"TMI_"*TMIversion*".mat.gz"),TMIfilegz,force=true)
    else
        !isfile(TMIfilegz) & !isfile(TMIfile) && google_download(url,pkgdatadir())
    end

    # cloak mat file in gz to get Google Drive spam filter to shut down
    isfile(TMIfilegz) & !isfile(TMIfile) && run(`gunzip $TMIfilegz`)

    # # make a sample field from zyx cartesian indices
    Izyx = cartesianindex(TMIfile)

    # # make a mask
    wet = falses(maximum(Izyx)[1],maximum(Izyx)[2],maximum(Izyx)[3])
    wet[Izyx] .= 1

    I = cartesianindex(wet)

    R = linearindex(wet)

    Azyx = watermassmatrix(TMIfile)
    A = matrix_zyx2xyz(TMIfile,Azyx,R)

    # # LU factorization for efficient matrix inversions
    Alu = lu(A)

    # get properties of grid
    lat,lon,depth = gridprops(TMIfile)

    γ = Grid(lon,lat,depth,I,R,wet)

    # need to make this optional
    L = circulationmatrix(TMIfile,γ)

    B = boundarymatrix(TMIfile,γ)

    # consider re-ordering this.
    # some output should be optional
    # return Izyx or I or neither?
    #return  A, Alu, γ, TMIfile, I, L, B
    return  A, Alu, γ, TMIfile, L, B
end

"""
    function cartesianindex(file)
    Read and assemble the grid coordinates
    according to the legacy MATLAB code (z,y,x order).
# Arguments
- `file`: TMI NetCDF file name
# Output
- `I`: TMI Cartesian index for wet points
"""
function cartesianindex(file::String)
    # make the Cartesian tracer grid
    if file[end-1:end] == "nc"

        it = convert(Vector{Int},ncread(file,"i"))
        jt = convert(Vector{Int},ncread(file,"j"))
        kt = convert(Vector{Int},ncread(file,"k"))

        I = CartesianIndex.(it,jt,kt)

    elseif file[end-2:end] == "mat"

        matobj = matopen(file)
        haskey(matobj,"it") ? it=convert(Vector{Integer},vec(read(matobj,"it"))) : it = convert(Vector{Integer},vec(read(matobj,"i")))
        haskey(matobj,"jt") ? jt=convert(Vector{Integer},vec(read(matobj,"jt"))) : jt=convert(Vector{Integer},vec(read(matobj,"j")))
        haskey(matobj,"kt") ? kt=convert(Vector{Integer},vec(read(matobj,"kt"))) : kt=convert(Vector{Integer},vec(read(matobj,"k")))
        close(matobj)
        I = CartesianIndex.(it,jt,kt)
    end
    return I
end

"""
    function gridprops(file)
    Read and assemble the grid properties.
# Arguments
- `file`: TMI NetCDF file name
# Output
- `grid`: TMI grid coordinates
"""
function gridprops(file)
    if file[end-1:end] == "nc"

        lon = convert(Vector{Float64},ncread(file,"lon"))
        lat = convert(Vector{Float64},ncread(file,"lat"))
        depth = convert(Vector{Float64},ncread(file,"depth"))

    elseif file[end-2:end] == "mat"

        matobj = matopen(file)
        lon=convert(Vector{Float64},vec(read(matobj,"LON")))
        lat=convert(Vector{Float64},vec(read(matobj,"LAT")))
        depth=convert(Vector{Float64},vec(read(matobj,"DEPTH")))
        close(matobj)

    end

    return lon,lat,depth
end

"""
    function watermassmatrix(file)
    Read and assemble the water-mass matrix.
# Arguments
- `file`: TMI NetCDF or MATLAB file name
# Output
- `A`: water-mass matrix
"""
function watermassmatrix(file)

    # consider adding a catch if A doesn't exist in file.

    if file[end-1:end] == "nc"
        # Int or Integer?
        i = convert(Vector{Int},ncread(file,"Arow"))
        j = convert(Vector{Int},ncread(file,"Acol"))
        m = ncread(file,"m")
        A = sparse(i,j,m)
    elseif file[end-2:end] == "mat"
        matobj = matopen(file)
        if haskey(matobj,"A")
            A=read(matobj,"A")
            close(matobj)
        else
            close(matobj)
            return nothing
        end

        # But MATLAB had zyx format and we need xyz format.
        # linearindices R not available so will do conversion in higher scope

    end
    return A
end

"""
        function matrix_zyx2xyz(TMIfile,Azyx,γ)

    Transfer zyx format water-mass matrix A to xyz format
# Arguments
- `Azyx`: water-mass matrix in zyx format
- `γ`: TMI grid
# Output
- `Axyz`: water-mass matrix in xyz format
"""
function matrix_zyx2xyz(file,Azyx,R)

    izyx, jzyx, mzyx = findnz(Azyx)
    Izyx = cartesianindex(file)

    # Julia accounting x,y,z
    ixyz = updatelinearindex(izyx,Izyx,R)
    jxyz = updatelinearindex(jzyx,Izyx,R)

    # use grid indices to switch i,j values
    Axyz = sparse(ixyz,jxyz,mzyx)
    return Axyz
end

"""
    function circulationmatrix(file,γ)
    Read and assemble the circulation matrix from MATLAB.
    Transfer to updated x,y,z version
# Arguments
- `file`: TMI MATLAB file name
- `γ`: TMI grid
# Output
- `L`: circulation matrix in xyz format
"""
function circulationmatrix(file,γ)

    if file[end-2:end] == "mat"

        matobj = matopen(file)
        if haskey(matobj,"L")
            # Matlab output in zyx format
            Lzyx=read(matobj,"L")
            close(matobj)

            Izyx = cartesianindex(file)
            izyx, jzyx, Fzyx = findnz(Lzyx)
            # Julia accounting x,y,z
            ixyz = updatelinearindex(izyx,Izyx,γ.R)
            jxyz = updatelinearindex(jzyx,Izyx,γ.R)
            L = sparse(ixyz,jxyz,Fzyx)

        else
            close(matobj)
            return nothing
        end

    elseif file[end-1:end] == "nc"

        # based on function arguments, read from inefficient storage of L matrix.
        if haskey(NCDataset(file),"F")

            i = convert(Vector{Int},ncread(file,"Lrow"))
            j = convert(Vector{Int},ncread(file,"Lcol"))
            F = ncread(file,"F")
            L = sparse(i,j,F)
        else
            return nothing
        end
    end

    return L
end

"""
    function circulationmatrix(file,A,γ)
    Read and assemble the circulation matrix from the efficient storage of A and F₀ variables.
# Arguments
- `file`: TMI MATLAB file name
- `A`: TMI water-mass matrix
- `γ`: TMI grid
# Output
- `L`: circulation matrix in xyz format
"""
function circulationmatrix(file,A,γ)

    file[end-1:end] !== "nc" && error("not a NetCDF file")

    # based on function arguments, read F₀ to efficiently reproduce L matrix.

    if haskey(NCDataset(file),"F₀")
        F₀ = ncread(file,"F₀")
        F₀vec = F₀[γ.wet]

        # For each row of A, multiply by F₀
        i, j, F = findnz(A)

        # careful, this loop can be really slow
        for nn in eachindex(i)
            F[nn] *= F₀vec[i[nn]]
        end

        L = sparse(i,j,F)
        return L
    else
        return nothing
    end

end

"""
        function boundarymatrix(file,γ)
    Read and assemble the boundary matrix from MATLAB.
    Transfer to updated x,y,z version
# Arguments
- `file`: TMI MATLAB file name
- `γ`: TMI grid
# Output
- `B`: boundary condition matrix
"""
function boundarymatrix(file,γ)

    if file[end-2:end] == "mat"

        matobj = matopen(file)
        if haskey(matobj,"B")
            Bzyx=read(matobj,"B")
            close(matobj)

            # matlab in zyx format.
            # consider using Azyx2xyz here.
            Izyx = cartesianindex(file)
            izyx, jzyx, Fzyx = findnz(Bzyx)
            # for B, rows are 3D grid space, columns are for the surface index.
            # Julia accounting x,y,z
            Isfc = surfaceindex(Izyx)
            ixyz = updatelinearindex(izyx,Izyx,γ.R)
            jxyz = updatelinearindex(Isfc[jzyx],Izyx,γ.R)

            # assume surface at k = 1 (revisit for LGM problem)
            # give the full dimension of sparse matrix
            B = sparse(ixyz,jxyz,Fzyx,sum(γ.wet),sum(γ.wet[:,:,1]))
        else
            close(matobj)
            return nothing
        end

    elseif file[end-1:end] == "nc"

        # based on function arguments, read from inefficient storage of L matrix.
        if haskey(NCDataset(file),"b")

            i = convert(Vector{Int},ncread(file,"Brow"))
            j = convert(Vector{Int},ncread(file,"Bcol"))
            b = ncread(file,"b")
            B = sparse(i,j,b,sum(γ.wet),sum(γ.wet[:,:,1]))

        else
            return nothing
        end
    end
    return B
end

"""
    function updatelinearindex(izyx,Izyx,R)
    Linear index translated from z,y,x to x,y,z accounting
# Arguments
- `izyx`: index of interest in z,y,x accounting
- `Izyx`: wet Cartesian Index for z,y,x
- `R`: Linear indices for x,y,z
# Output
- `ixyz`: index of interest in x,y,z accounting
"""
function updatelinearindex(izyx,Izyx,R)
    # get Izyx Cartesian index stored from legacy MATLAB code
    ixyz = R[Izyx[izyx]]
    return ixyz
end

"""
     function surfaceregion(TMIversion::String,region::String,γ::Grid)::BoundaryCondition

    Read an oceanographically-relevant surface region from NetCDF file. (Also could be read from mat file.)
    Return a BoundaryCondition
"""
function surfaceregion(TMIversion::String,region::String,γ::Grid)::BoundaryCondition

    file = pkgdatadir("TMI_"*TMIversion*".nc")
    tracername = "d_"*region

    # Didn't use readfiled because dsfc is 2d.
    dsfc = ncread(file,tracername)
    b = BoundaryCondition(dsfc,γ.lon,γ.lat,γ.depth[1],3,1,γ.wet[:,:,1])
    return b
end

"""
    function ncurl(TMIversion)
    placeholder function to give location (URL) of NetCDF Google Drive input
    in the future, consider a struct or Dict that describes all TMI versions.
# Arguments
- `TMIversion`: version of TMI water-mass/circulation model
# Output
- `url`: location (URL) for download
"""
function ncurl(TMIname)
    if TMIname == "modern_90x45x33_GH10_GH12"
        url = "https://docs.google.com/uc?export=download&id=1Fn_cY-90_RDbBGh6kV0kpXmsvwdjp1Cd"
    elseif TMIname == "modern_180x90x33_GH11_GH12"
        url = "https://docs.google.com/uc?export=download&id=1-YEkB_YeQGqPRH6kauhBb2bi_BjVGt9b"
    elseif TMIname == "modern_90x45x33_unpub12"
        url = "https://docs.google.com/uc?export=download&id=1Kw_Mr7fiKqan0nx0dKvGHnSInP0hQ7AV"
    elseif TMIname == "modern_90x45x33_G14"
        url = "https://docs.google.com/uc?export=download&id=1aeE7EXA-vy3Cm_drt4qCFw4AlpYrdudk"
    elseif TMIname == "modern_90x45x33_G14_v2"
        url = "https://docs.google.com/uc?export=download&id=1Mwhv70soBX6-pYijU0ElNl0TZw0vSbXN"
    elseif TMIname == "LGM_90x45x33_G14"
        url = "https://docs.google.com/uc?export=download&id=1yoDi7_foBt3TVULCstlWnNLHFc2G47Fz"
    elseif TMIname == "LGM_90x45x33_G14A"
        url = "https://docs.google.com/uc?export=download&id=1ADkDI3Fc3z4Vm75K5u6hx0Yu1P0iVnW1"
    elseif TMIname == "LGM_90x45x33_GPLS1"
        url = "https://docs.google.com/uc?export=download&id=1VOrZGUsO7lp21qw6Yw0dBqGlIy_ImdRW"
    elseif TMIname == "LGM_90x45x33_GPLS2"
        url = "https://docs.google.com/uc?export=download&id=1cOCrty9kvA2s3NoD1QZjnNehlVbP0rHP"
    elseif TMIname == "LGM_90x45x33_OG18"
        url = "https://docs.google.com/uc?export=download&id=19zccG1BSdspD9rti2OttsF2Dm4P2OLjt"
    else
        url = nothing
    end
end

"""
    function maturl(TMIversion)
    Find *mat file here.
    placeholder function to give location (URL) of Google Drive input
    in the future, consider a struct or Dict that describes all TMI versions.
# Arguments
- `TMIversion`: version of TMI water-mass/circulation model
# Output
- `url`: location (URL) for download
"""
function maturl(TMIname)
    if TMIname == "modern_90x45x33_GH10_GH12"
        url = "https://docs.google.com/uc?export=download&id=1Z2knDctAmZHO2lcWTBCdR8zjkmbcyCGg"
    elseif TMIname == "modern_180x90x33_GH11_GH12"
        url = "https://docs.google.com/uc?export=download&id=11zD1nOfT6V7G0qIHdjK2pDGHFk-ExXwU"
    elseif TMIname == "modern_90x45x33_unpub12"
        url = "https://docs.google.com/uc?export=download&id=1sqkjFCPxZT_2Bm9rsp0acyxxkBri9YAT"
    elseif TMIname == "modern_90x45x33_G14"
        url = "https://docs.google.com/uc?export=download&id=1dCrDe5VXrsXiOf04mbuHID7xc5Ymm8-z"
    elseif TMIname == "modern_90x45x33_G14_v2"
        url = "https://docs.google.com/uc?export=download&id=1Axaqn88HZv3i4rYn0g5dUztVXlx_TmGt"
    elseif TMIname == "LGM_90x45x33_G14"
                url = "https://docs.google.com/uc?export=download&id=1qWnL9SrcjFRt4TUQ0k7v0jEFrKtUYpev"
    elseif TMIname == "LGM_90x45x33_G14A"
                url = "https://docs.google.com/uc?export=download&id=1E7qYWdAnoz4YXFvJ0AkVcqpHDq8aFP4R"
    elseif TMIname == "LGM_90x45x33_GPLS1"
                url = "https://docs.google.com/uc?export=download&id=1nF5KjWE3n--NkOTgIymIhP6pkafJ7jl3"
    elseif TMIname == "LGM_90x45x33_GPLS2"
                url = "https://docs.google.com/uc?export=download&id=1CHiR2J60HPRrxXTvJupyWEHmggMkgnxV"
    elseif TMIname == "LGM_90x45x33_OG18"
                url = "https://docs.google.com/uc?export=download&id=1WKAWxDYDFls6C5YzO1sMLDlUhC8_tKV2"
    else
        url = nothing
    end
    return url
end

"""
Save TMI configuration to NetCDF format for non-proprietary access
"""
function config2nc(TMIversion,A,γ,L,B)

    # make new netcdf file.
    filenetcdf = pkgdatadir("TMI_"*TMIversion*".nc")
    isfile(filenetcdf) && rm(filenetcdf)

    grid2nc(TMIversion,γ)

    matfields2nc(TMIversion,γ)

    !isnothing(A) && watermassmatrix2nc(TMIversion,A)

    !isnothing(L) && circulationmatrix2nc(TMIversion,L,γ)

    !isnothing(B) && boundarymatrix2nc(TMIversion,B)

    #= is this part of the config? Or should it go to
     a separate output? It is similar to the output fields above. Probably should be considered part of the config. =#
    regions2nc(TMIversion,γ)

    optim2nc(TMIversion)

end

"""
Save grid dictionaries of attributes for writing to NetCDF file
"""
function griddicts(γ)
    # update names and types in dictionary

    TMIgrids = Dict("lon" => γ.lon,
                    "lat" => γ.lat,
                    "depth" => γ.depth)

    TMIgridsatts = Dict("lon" => Dict("longname" => "Longitude", "units" => "°E"),
                        "lat" => Dict("longname" => "Latitude", "units" => "°N"),
                        "depth" => Dict("longname" => "depth", "units" => "m"))

    return TMIgrids, TMIgridsatts

end

"""
Read 3D fields from mat file and save to NetCDF file.
"""
function matfields2nc(TMIversion,γ)

    filenetcdf = pkgdatadir("TMI_"*TMIversion*".nc")
    filemat = pkgdatadir("TMI_"*TMIversion*".mat")
    vars = matread(filemat)

    TMIgrids, TMIgridsatts = griddicts(γ)

    T = eltype(γ.lon) # does the eltype of longitude have to equal the tracer eltype?
    #T =  Float64

    varlist = Dict("dP"=>"qPO₄","q"=>"qPO₄",
                   "Tobs"=>"θ","Tmod"=>"θ","Tlgm"=>"θ",
                   "Terr"=>"σθ",
                   "Sobs"=>"Sp","Smod"=>"Sp","Slgm"=>"Sp",
                   "Serr"=>"σSp",
                   "O18obs"=>"δ¹⁸Ow","O18mod"=>"δ¹⁸Ow","O18lgm"=>"δ¹⁸Ow",
                   "O18err"=>"σδ¹⁸Ow",
                   "Pobs"=>"PO₄","Pmod"=>"PO₄","Plgm"=>"PO₄",
                   "Perr" => "σPO₄",
                   "Nobs"=>"NO₃","Nmod"=>"NO₃","Nlgm"=>"NO₃",
                   "Nerr" => "σNO₃",
                   "Oobs"=>"O₂","Omod"=>"O₂","Olgm"=>"O₂",
                   "Oerr"=>"σO₂",
                   "C13obs"=>"δ¹³C","C13mod"=>"δ¹³C","C13lgm"=>"δ¹³C",
                   "C13err" =>  "σδ¹³C")

    # iterate over all possible variables listed above
    Izyx = cartesianindex(filemat)
    TMIfields = Dict{String,Array{T,3}}()
    for (kk,vv) in varlist
        haskey(vars,kk) ? push!(TMIfields, vv => tracerinit(vars[kk], Izyx, γ.wet)) : nothing
    end

    # also save fields that are stored in the x struct, if they exist
    if haskey(vars,"x")
        for (kk,vv) in varlist
            haskey(vars["x"],kk) ? push!(TMIfields, vv => tracerinit(vars["x"][kk], Izyx, γ.wet)) : nothing
        end
    end

    TMIfieldsatts = fieldsatts()

    # iterate in TMIgrids Dictionary to write to NetCDF.
    for (varname,varvals) in TMIfields

        nccreate(filenetcdf,varname,"lon",γ.lon,TMIgridsatts["lon"],"lat",γ.lat,TMIgridsatts["lat"],"depth",γ.depth,TMIgridsatts["depth"],atts=TMIfieldsatts[varname])
        println("write ",varname)
        ncwrite(varvals,filenetcdf,varname)

    end
end

"""
All variable names and attributes.
Useful for writing NetCDF files.
"""
fieldsatts() =
    Dict("θ" => Dict("longname" => "potential temperature", "units" => "°C"),
         "σθ" => Dict("longname" => "1σ standard error in potential temperature", "units" => "°C"),
         "Sp" => Dict("longname" => "practical salinity", "units" => "PSS-78"),
         "σSp" => Dict("longname" => "1σ standard error in practical salinity", "units" => "PSS-78"),
         "δ¹⁸Ow" => Dict("longname" => "oxygen-18 to oxygen-16 ratio in seawater", "units" => "‰ VSMOW"),
         "σδ¹⁸Ow" => Dict("longname" => "1σ standard error in oxygen-18 to oxygen-16 ratio in seawater", "units" => "‰ VSMOW"),
         "PO₄" => Dict("longname" => "phosphate", "units" => "μmol/kg"),
         "σPO₄" => Dict("longname" => "1σ standard error in phosphate", "units" => "μmol/kg"),
         "qPO₄" => Dict("longname" => "local source of phosphate", "units" => "μmol/kg"),
         "NO₃" => Dict("longname" => "nitrate", "units" => "μmol/kg"),
         "σNO₃" => Dict("longname" => "1σ standard error in nitrate", "units" => "μmol/kg"),
         "O₂" => Dict("longname" => "dissolved oxygen", "units" => "μmol/kg"),
         "σO₂" => Dict("longname" => "1σ standard error in dissolved oxygen", "units" => "μmol/kg"),
         "δ¹³C" => Dict("longname" => "carbon-13 to carbon-12 ratio in DIC", "units" => "‰ PDB"),
         "σδ¹³C" => Dict("longname" => "1σ standard error fin carbon-13 to carbon-12 ratio in DIC", "units" => "‰ PDB"),
         "F₀" => Dict("longname" => "normalized mass flux out of gridcell", "units" => "(kg seawater/s)/(kg gridcell)"))

"""
Read vectors from mat file, translate to 3D,
 and save surface field to NetCDF file.
"""
function regions2nc(TMIversion,γ)

    filenetcdf = pkgdatadir("TMI_"*TMIversion*".nc")
    filemat = pkgdatadir("TMI_"*TMIversion*".mat")

    # region names
    # didn't figure out how to use an ordered dict, instead use a tuple
    list = ("GLOBAL","ANT","SUBANT",
            "NATL","NPAC","TROP","ARC",
            "MED","ROSS","WED","LAB","GIN",
            "ADEL","SUBANTATL","SUBANTPAC","SUBANTIND",
            "TROPATL","TROPPAC","TROPIND")

    regionname = Dict("GLOBAL" => "globally uniform",
                      "ANT" => "Antarctic",
                      "SUBANT" => "Subantarctic",
                      "NATL" => "North Atlantic",
                      "NPAC" => "North Pacific",
                      "TROP" => "tropical and subtropical",
                      "ARC" => "Arctic",
                      "MED" => "Mediterranean",
                      "ROSS" => "Ross Sea sector",
                      "WED" => "Weddell Sea sector",
                      "LAB" => "Labrador and Irminger Seas",
                      "GIN" => "Greenland-Iceland-Norwegian Seas",
                      "ADEL" => "Adélie Land sector",
                      "SUBANTATL" => "Atlantic-sector Subantarctic",
                      "SUBANTPAC" => "Pacific-sector Subantarctic",
                      "SUBANTIND" => "Indian-sector Subantarctic",
                      "TROPATL" => "tropical and subtropical Atlantic",
                      "TROPPAC" => "tropical and subtropical Pacific",
                      "TROPIND" => "tropical and subtropical Indian")

    matobj = matopen(filemat)
    if haskey(matobj,"d_all")
        d_all = read(matobj,"d_all")
        close(matobj)
    else
        return
    end

    # a kludge for now
    T = eltype(γ.lon)

    # iterate over all regions in d_all
    Izyx = cartesianindex(filemat)
    regions = Dict{String,Array{T,2}}()
    regionatts = Dict{String,Dict{String,String}}()

    for rr = 1:size(d_all,2)
        # 3D fields in zyx vector format
        # are changed to 3D xyz format
        d = tracerinit(d_all[:,rr],Izyx,γ.wet)

        # just save the surface 2D field
        push!(regions, list[rr] => d[:,:,1])

        push!(regionatts, list[rr] =>
         Dict("longname" => regionname[list[rr]]*" surface region", "units" => "[]"))
    end

    TMIgrids, TMIgridsatts = griddicts(γ)

    # iterate in regions Dictionary to write to NetCDF.
    for (varname,varvals) in regions
        dvarname = "d_"*varname
        nccreate(filenetcdf,dvarname,"lon",γ.lon,TMIgridsatts["lon"],"lat",γ.lat,TMIgridsatts["lat"],atts=regionatts[varname])
        println("write ",dvarname)
        ncwrite(varvals,filenetcdf,dvarname)

    end
end

function watermassmatrix2nc(TMIversion,A)

    filenetcdf = pkgdatadir("TMI_"*TMIversion*".nc")
    i, j, m = findnz(A)
    nelements = length(i)

    # add the circulation matrix: problem can't store sparse matrix.
    varname= "m"
    elementatts = Dict("longname" => "TMI sparse matrix element number")
    matts =  Dict("longname" => "TMI water-mass fraction (sparse matrix values)", "units" => "(kg source)/(kg total)")
    nccreate(filenetcdf,varname,"A_element",1:nelements,elementatts,atts=matts)
    println("write ",varname)
    ncwrite(m, filenetcdf,varname)

     varname= "Arow"
     destatts = Dict("longname" => "gridcell number of destination (row value)")
     nccreate(filenetcdf,varname,"A_element",1:nelements,elementatts,atts=destatts)
    println("write ",varname)
    ncwrite(i, filenetcdf,varname)

    varname= "Acol"
    sourceatts = Dict("longname" => "gridcell number of source (column value)")
    nccreate(filenetcdf,varname,"A_element",1:nelements,elementatts,atts=sourceatts)
    println("write ",varname)
    ncwrite(j, filenetcdf,varname)

end

"""
Save optimization parameters to NetCDF file)

Future considerations: split into 2 functions
1) read from mat
2) save to nc
"""
function optim2nc(TMIversion)

    filemat = pkgdatadir("TMI_"*TMIversion*".mat")
    filenetcdf = pkgdatadir("TMI_"*TMIversion*".nc")

    matobj = matopen(filemat)
    if haskey(matobj,"fval")
        J̃ = read(matobj,"fval")
        iteratts = Dict("longname" => "iteration number")
        Jatts =  Dict("longname" => "cost function value", "units" => "[]")

        varname = "J" # J̃ not output to screen properly
        println("write ",varname)
        nccreate(filenetcdf,varname,"iter",1:length(J̃),iteratts,atts=Jatts)
    end
    if haskey(matobj,"u")
        ũ = read(matobj,"u")
        iteratts = Dict("longname" => "control element number")
        uatts =  Dict("longname" => "control vector", "units" => "[]")

        varname = "ũ"
        println("write ",varname)
        nccreate(filenetcdf,varname,"control_element",1:length(ũ),iteratts,atts=uatts)
    end
    close(matobj)

end

"""
Save circulation matrix `L` to NetCDF file.
"""
function circulationmatrix2nc(TMIversion,L,γ)

    T = eltype(L)
    fullmatrix = false # more efficient to just save F₀, then modify A to get L
    filenetcdf = pkgdatadir("TMI_"*TMIversion*".nc")
    if !fullmatrix
        F₀ = tracerinit(γ.wet,T)
        for nd ∈ eachindex(γ.I)
            # normalized mass flux out of gridcell is found on diagonal
            F₀[γ.I[nd]] = -L[nd,nd]
        end

        TMIgrids, TMIgridsatts = griddicts(γ)
        TMIfieldsatts = fieldsatts()
        varname = "F₀"
        nccreate(filenetcdf,varname,"lon",γ.lon,TMIgridsatts["lon"],"lat",γ.lat,TMIgridsatts["lat"],"depth",γ.depth,TMIgridsatts["depth"],atts=TMIfieldsatts[varname])
        println("write ",varname)
        ncwrite(F₀,filenetcdf,varname)

    else # fullmatrix = true, makes a bigger nc file

        i, j, F = findnz(L)
        nelements = length(i)

        # add the circulation matrix: problem can't store sparse matrix.
        varname = "F"
        elementatts = Dict("longname" => "TMI sparse matrix element number")
        Fatts =  Dict("longname" => "normalized mass flux (sparse matrix values)","units" => "(kg seawater/s)/(kg gridcell)")
        nccreate(filenetcdf,varname,"L_element",1:nelements,elementatts,atts=Fatts)
        println("write ",varname)
        ncwrite(F, filenetcdf,varname)

        varname= "Lrow"
        destatts = Dict("longname" => "gridcell number of destination (row value)")
        nccreate(filenetcdf,varname,"L_element",1:nelements,elementatts,atts=destatts)
        println("write ",varname)
        ncwrite(i, filenetcdf,varname)

        varname= "Lcol"
        sourceatts = Dict("longname" => "gridcell number of source (column value)")
        nccreate(filenetcdf,varname,"L_element",1:nelements,elementatts,atts=sourceatts)
        println("write ",varname)
        ncwrite(j, filenetcdf,varname)
    end

end

"""
Save boundary matrix for transient model to NetCDF file
"""
function boundarymatrix2nc(TMIversion,B)

    filenetcdf = pkgdatadir("TMI_"*TMIversion*".nc")
    i, j, b = findnz(B)
    nelements = length(i)

    # add the circulation matrix: problem can't store sparse matrix.
    varname= "b"
    elementatts = Dict("longname" => "TMI boundary matrix element number")
    matts =  Dict("longname" => "Boundary matrix values", "units" => "[]")
    nccreate(filenetcdf,varname,"B_element",1:nelements,elementatts,atts=matts)
    println("write ",varname)
    ncwrite(b, filenetcdf,varname)

    varname= "Brow"
    destatts = Dict("longname" => "gridcell number of 3D field")
    nccreate(filenetcdf,varname,"B_element",1:nelements,elementatts,atts=destatts)
    println("write ",varname)
    ncwrite(i, filenetcdf,varname)

    varname= "Bcol"
    sourceatts = Dict("longname" => "gridcell number of surface boundary condition")
    nccreate(filenetcdf,varname,"B_element",1:nelements,elementatts,atts=sourceatts)
    println("write ",varname)
    ncwrite(j, filenetcdf,varname)

end

"""
Put grid properties (Cartesian index) into NetCDF file
"""
function grid2nc(TMIversion,γ)

    filenetcdf = pkgdatadir("TMI_"*TMIversion*".nc")

    linearindexatts = Dict("longname" => "linear index")
    nfld = length(γ.I)

    iatts =  Dict("longname" => "Cartesian index for x-direction", "units" => "[]")
    jatts =  Dict("longname" => "Cartesian index for y-direction", "units" => "[]")
    katts =  Dict("longname" => "Cartesian index for z-direction", "units" => "[]")

    varname = "i"
    nccreate(filenetcdf,varname,"linearindex",1:nfld,linearindexatts,atts=iatts)
    ncwrite(lonindex(γ.I),filenetcdf,varname)

    varname = "j"
    nccreate(filenetcdf,varname,"linearindex",1:nfld,linearindexatts,atts=jatts)
    ncwrite(latindex(γ.I),filenetcdf,varname)

    varname = "k"
    nccreate(filenetcdf,varname,"linearindex",1:nfld,linearindexatts,atts=katts)
    ncwrite(depthindex(γ.I),filenetcdf,varname)

end


text:
using HorizonSideRobots

function count_barriers(r::Robot)
    ans = 0
    move_while_can(r, (West, Sud))
    dir = Ost
    ans += process_row(r, dir, Nord)
    while !isborder(r, Nord)
        move!(r, Nord)
        dir = inverse(dir)
        ans += process_row(r, dir, Nord)
    end
    move_while_can(r, (West, Nord))

    dir = Sud
    ans += process_row(r, dir, Ost)

    while !isborder(r, Ost)
        move!(r, Ost)
        dir = inverse(dir)
        ans += process_row(r, dir, Ost)
    end

    return ans - 2
end

function process_row(r::Robot, side, toCheckSide)
    cnt = 0
    flag = false
    while !isborder(r, side) || bypass(r, side)
        if isborder(r, toCheckSide)
            if !flag
                flag = true
                cnt += 1
            end
        else
            flag = false
        end
        if !isborder(r, side)
            move!(r, side)
        end
    end
    return cnt
end

function bypass(r::Robot, dir)
    bypass_dir = nextCounterClockwise(dir)
    inverse_bypass_dir = inverse(bypass_dir)
    num_steps = 0
    while !isborder(r, bypass_dir) && isborder(r, dir)
        move!(r, bypass_dir)
        num_steps += 1
    end
    if isborder(r, dir) && isborder(r,bypass_dir)
        moves!(r, inverse_bypass_dir, num_steps)
        return false
    end
    move!(r, dir)
    moves!(r, inverse_bypass_dir, num_steps)
    return true
    end

# move while it is possible for us
function move_while_can(r::Robot, sides)
    moving = true
    # it is possible to move
    while moving
        moving = false
        for side in sides
            if !isborder(r, side)
                move!(r, side)
                moving = true
            end
        end
    end
    end

###-------------------------------------------------------------------------
function moves!(r::Robot, side::HorizonSide)
    while isborder(r, side) == false
        move!(r, side)
    end
end
function moves!(r::Robot, side::HorizonSide, num_steps::Int)
    for _ in 1:num_steps
        move!(r, side)
    end
end

function putmarkers!(r::Robot, side::HorizonSide)
    while isborder(r, side) == false
        move!(r, side)
        putmarker!(r)
    end

function move_markers!(r::Robot, side::HorizonSide)
    while ismarker(r) == true
        move!(r, side)
    end
end

left(side::HorizonSide) = HorizonSide(mod(Int(side) + 1, 4))

inverse(side::HorizonSide) = HorizonSide(mod(Int(side) + 2, 4))

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
# export tutte

"""
`tutte(G::SimpleGraph, outside)` gives `G` a Tutte embedding in which
the list of vertices in `outside` form the outer boundary. The graph
should be connected (ideally, 3-connected) and, if planar and `outside`
defines a face, the embedding will be crossing free.

`tutte(G::SimpleGraph)` assumes `G` has a rotation system in which case a
largest face will be selected to be `outside`.
"""
function _tutte(G::SimpleGraph{T}, outside::Vector{T}) where {T}
    if !issubset(outside, G.V)
        error("Some of the proposed outer vertices are not in this graph")
    end
    n = NV(G)
    VV = vlist(G)

    outside = unique(outside)  # remove dups

    lookup(v) = findfirst([v == w for w in VV])

    A = zeros(n, n)
    r = sqrt(n)

    rhs_x = zeros(n)
    rhs_y = zeros(n)

    no = length(outside)
    for j = 1:no
        theta = 2 * pi * (j - 1) / no + pi / 2
        v = outside[j]
        idx = lookup(v)
        rhs_x[idx] = r * cos(theta)
        rhs_y[idx] = r * sin(theta)
    end

    for v in G.V
        idx = lookup(v)
        A[idx, idx] = 1

        if !in(v, outside)
            d = deg(G, v)
            for w in G[v]
                j = lookup(w)
                A[idx, j] = -1 / d
            end
        end
    end

    x = A \ rhs_x
    y = A \ rhs_y

    d = Dict{T,Vector{Float64}}()
    for v in VV
        idx = lookup(v)
        d[v] = [x[idx]; y[idx]]
    end
    embed(G, d)
end

function _tutte(G::SimpleGraph)
    FF = faces(G)
    Fmax = first(FF)
    for F in FF
        if length(F) > length(Fmax)
            Fmax = F
        end
    end
    _tutte(G, first.(Fmax))
end


text:
"""
    give the error text as string
"""
function err_str(err; max_len = 10000)
    s = sprint(showerror, err, catch_backtrace())
    return length(s) > max_len ? s[1:max_len] * "\n[...]" : s
end

## ----------------------------------------------------------------------------
function _extract_dir(args...)
    isempty(args) && return ("", args)
    path = ""
    for (i, arg) in enumerate(args)
        !(arg isa Vector{<:AbstractString}) && return (path, args[i:end])
        for dir in arg
            path = joinpath(path, dir)
        end
    end
    return (path, tuple())
end