text:
showsafe(x) = (x === nothing) ? "nothing" : x

function gitcmd(path::String, gitconfig::Dict)
    cmd = ["git", "-C", path]
    for (n,v) in gitconfig
        push!(cmd, "-c")
        push!(cmd, "$n=$v")
    end
    Cmd(cmd)
end

"""
Return a `GitRepo` object for an up-to-date copy of `registry`.
"""
function get_registry(registry::String; gitconfig::Dict=Dict())
    reg_path(args...) = joinpath("registries", map(string, args)...)
    if haskey(REGISTRIES, registry)
        registry_uuid = REGISTRIES[registry]
        registry_path = reg_path(registry_uuid)
        if !ispath(registry_path)
            LibGit2.clone(registry, registry_path, branch="master")
        else
            # this is really annoying/impossible to do with LibGit2
            git = gitcmd(registry_path, gitconfig)
            run(`$git config remote.origin.url $registry`)
            run(`$git checkout -q -f master`)
            run(`$git fetch -q -P origin master`)
            run(`$git reset -q --hard origin/master`)
        end
    else
        registry_temp = mktempdir(mkpath(reg_path()))
        try
            LibGit2.clone(registry, registry_temp)
            reg = TOML.parsefile(joinpath(registry_temp, "Registry.toml"))
            registry_uuid = REGISTRIES[registry] = UUID(reg["uuid"])
            registry_path = reg_path(registry_uuid)
            rm(registry_path, recursive=true, force=true)
            mv(registry_temp, registry_path)
        finally
            rm(registry_temp, recursive=true, force=true)
        end
    end
    return GitRepo(registry_path)
end

"""
Given a remote repo URL and a git tree spec, get a `Project` object
for the project file in that tree and a hash string for the tree.
"""
# function get_project(remote_url::String, tree_spec::String)
#     # TODO?: use raw file downloads for GitHub/GitLab
#     mktempdir(mkpath("packages")) do tmp
#         # bare clone the package repo
#         @debug("bare clone the package repo")
#         repo = LibGit2.clone(remote_url, joinpath(tmp, "repo"), isbare=true)
#         tree = try
#             LibGit2.GitObject(repo, tree_spec)
#         catch err
#             err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
#             error("$remote_url: git object $(repr(tree_spec)) could not be found")
#         end
#         tree isa LibGit2.GitTree || (tree = LibGit2.peel(LibGit2.GitTree, tree))
#
#         # check out the requested tree
#         @debug("check out the requested tree")
#         tree_path = abspath(tmp, "tree")
#         GC.@preserve tree_path begin
#             opts = LibGit2.CheckoutOptions(
#                 checkout_strategy = LibGit2.Consts.CHECKOUT_FORCE,
#                 target_directory = Base.unsafe_convert(Cstring, tree_path)
#             )
#             LibGit2.checkout_tree(repo, tree, options=opts)
#         end
#
#         # look for a project file in the tree
#         @debug("look for a project file in the tree")
#         project_file = Pkg.Types.projectfile_path(tree_path)
#         project_file !== nothing && isfile(project_file) ||
#             error("$remote_url: git tree $(repr(tree_spec)) has no project file")
#
#         # parse the project file
#         @debug("parse the project file")
#         project = Pkg.Types.read_project(project_file)
#         project.name === nothing &&
#             error("$remote_url $(repr(tree_spec)): package has no name")
#         project.uuid === nothing &&
#             error("$remote_url $(repr(tree_spec)): package has no UUID")
#         project.version === nothing &&
#             error("$remote_url $(repr(tree_spec)): package has no version")
#
#         return project, string(LibGit2.GitHash(tree))
#     end
# end

"""
Write TOML data (with sorted keys).
"""
function write_toml(file::String, data::Dict)
    open(file, "w") do io
        TOML.print(io, data, sorted=true)
    end
end

function write_registry(file::String, data::Dict)
    data2 = deepcopy(data)
    delete!(data2, "packages")
    open(file, "w") do io
        TOML.print(io, data2, sorted=true)
        print(io, "\n[packages]\n")
        for key in sort(collect(keys(data["packages"])))
            name_str = string("\"", data["packages"][key]["name"], "\"")
            path_str = string("\"", escape_string(data["packages"][key]["path"]), "\"")
            path_str = replace(path_str, "\\\\" => "/")
            println(io, key, " = { name = $(name_str), path = $(path_str) }")
        end
    end
end

struct RegBranch
    name::String
    version::VersionNumber
    branch::String

    error::Union{Nothing, String}
end

function update_package_data(pkg::Pkg.Types.Project, registry_path, package_repo, package_path, tree_hash)
    # package file
    @debug("update package data: package file")
    package_info = filter(((k,v),)->!(v isa Dict), Pkg.Types.destructure(pkg))
    delete!(package_info, "version")
    package_info["repo"] = package_repo
    package_file = joinpath(package_path, "Package.toml")
    write_toml(package_file, package_info)

    # versions file
    @debug("update package data: versions file")
    versions_file = joinpath(package_path, "Versions.toml")
    versions_data = isfile(versions_file) ? TOML.parsefile(versions_file) : Dict()

    try
        versions = sort!([VersionNumber(v) for v in keys(versions_data)])
        #Base.check_new_version(versions, pkg.version)
    catch ex
        if isa(ex, ErrorException)
            return RegBranch(pkg.name, pkg.version, branch, ex.msg)
        else
            rethrow(ex)
        end
    end

    version_info = Dict{String,Any}("git-tree-sha1" => string(tree_hash))
    versions_data[string(pkg.version)] = version_info

    vnlist = sort([(VersionNumber(k), v) for (k, v) in versions_data])
    vslist = [(string(k), v) for (k, v) in vnlist]

    open(versions_file, "w") do io
        TOML.print(io, OrderedDict(vslist))
    end

    # deps file
    @debug("update package data: deps file")
    deps_file = joinpath(package_path, "Deps.toml")
    if isfile(deps_file)
        deps_data = Pkg.Compress.load(deps_file)
    else
        deps_data = Dict()
    end

    registry_file = joinpath(registry_path, "Registry.toml")
    registry_data = TOML.parsefile(registry_file)

    @debug("Verifying package name and uuid in deps")
    for (k, v) in pkg.deps
        u = string(v)
        if haskey(registry_data["packages"], u)
            name_in_reg = registry_data["packages"][u]["name"]
            if name_in_reg != k
                err = "Error in `[deps]`: UUID $u refers to package '$name_in_reg' in registry but deps file has '$k'"
                throw(err)
            end
        elseif haskey(BUILTIN_PKGS, k)
            if BUILTIN_PKGS[k] != u
                err = "Error in `[deps]`: UUID $u for package $k should be $(BUILTIN_PKGS[k])"
                throw(err)
            end
        else
#             err = "Error in `[deps]`: Package '$k' with UUID: $u not found in registry or stdlib"
#             throw(err)
        end
    end

    deps_data[pkg.version] = pkg.deps
    Pkg.Compress.save(deps_file, deps_data)

    # compat file
    @debug("check compat section")
    for (p, v) in pkg.compat
        try
            ver = Pkg.Types.semver_spec(v)
            if p == "julia" && any(map(x->!isempty(intersect(Pkg.Types.VersionRange("0-0.6"),x)), ver.ranges))
                err = "Julia version < 0.7 not allowed in `[compat]`"
                @debug(err)
                throw(err)
            end
        catch ex
            if isa(ex, ArgumentError)
                err = "Error in `[compat]`: $(ex.msg)"
                @debug(err)
                throw(err)
            else
                rethrow(ex)
            end
        end
    end

    @debug("update package data: compat file")
    compat_file = joinpath(package_path, "Compat.toml")
    if isfile(compat_file)
        compat_data = Pkg.Compress.load(compat_file)
    else
        compat_data = Dict()
    end
    compat_data[pkg.version] = Dict{String,Any}(n=>Pkg.Types.semver_spec(v) for (n,v) in pkg.compat)
    Pkg.Compress.save(compat_file, compat_data)

    return
end

function get_package_path(registry_path, pkg::Pkg.Types.Project)
    @debug("find package in registry")
    registry_file = joinpath(registry_path, "Registry.toml")
    registry_data = TOML.parsefile(registry_file)

    uuid = string(pkg.uuid)
    if haskey(registry_data["packages"], uuid)
        package_data = registry_data["packages"][uuid]
        if (package_data["name"] != pkg.name)
            err = "Changing package names not supported yet"
            @debug(err)
            throw(err)
        end
        package_path = joinpath(registry_path, package_data["path"])
    else
        @debug("Package with UUID: $uuid not found in registry, checking if UUID was changed")
        for (k, v) in registry_data["packages"]
            if v["name"] == pkg.name
                err = "Changing UUIDs is not allowed"
                @debug(err)
                throw(err)
            end
        end

        @debug("Creating directory for new package $(pkg.name)")
        first_letter = uppercase(pkg.name[1])
        relative_package_path = joinpath("$first_letter", pkg.name)
        package_path = joinpath(registry_path, relative_package_path)
        mkpath(package_path)

        registry_data["packages"][uuid] = Dict("name" => pkg.name,
                                               "path" => relative_package_path)
        write_registry(registry_file, registry_data)
    end

    return package_path
end

function commit_registry(pkg::Pkg.Types.Project, package_path, package_repo, tree_hash, git)
    @debug("commit changes")
    message = """
    New version: $(pkg.name) v$(pkg.version)

    UUID: $(pkg.uuid)
    Repo: $(package_repo)
    Tree: $(string(tree_hash))
    """
    run(`$git add --all`)
    run(`$git commit -qm $message`)
end

function get_tree_hash(package_path, gitconfig)
    # repo = LibGit2.GitRepo(package_repo)
    # commit = LibGit2.GitObject(repo, "HEAD")
    # tree = LibGit2.peel(commit)
    # tree_hash = string(LibGit2.GitHash(git_tree))

    git = gitcmd(package_path, gitconfig)
    return read(`$git log --pretty=format:%T -1`, String)
end

function get_remote_repo(package_path, gitconfig)
    git = gitcmd(package_path, gitconfig)
    remote_name = split(chomp(read(`$git remote`, String)), '\n')
    length(remote_name) > 1 && throw("Repository has multiple remotes.")
    # TODO: Handle this case?
    remote_name[1] == "" && throw("Repository does not have a remote.")
    return readchomp(`$git remote get-url $(remote_name[1])`)
end

"""
Register the package at `package_repo` / `tree_spect` in `registry`.
"""
function register(
    package_repo::String, pkg::Pkg.Types.Project, tree_hash::String;
    registry::String = DEFAULT_REGISTRY,
    push::Bool = false,
    gitconfig::Dict = Dict()
)
    # get info from package registry
    @debug("get info from package registry")
    package_repo = GitTools.normalize_url(package_repo)
    #pkg, tree_hash = get_project(package_repo, tree_spec)
    branch = "register/$(pkg.name)/v$(pkg.version)"

    # get up-to-date clone of registry
    @debug("get up-to-date clone of registry")
    registry = GitTools.normalize_url(registry)
    registry_repo = get_registry(registry; gitconfig=gitconfig)
    registry_path = LibGit2.path(registry_repo)

    clean_registry = true
    err = nothing
    try
        # branch registry repo
        @debug("branch registry repo")
        git = gitcmd(registry_path, gitconfig)
        run(`$git checkout -qf master`)
        run(`$git branch -qf $branch`)
        run(`$git checkout -qf $branch`)

        try
            package_path = get_package_path(registry_path, pkg)
            update_package_data(pkg, registry_path, package_repo, package_path, tree_hash)
            commit_registry(pkg, package_path, package_repo, tree_hash, git)
        catch err
            return RegBranch(pkg.name, pkg.version, branch, err)
        end

        # push -f branch to remote
        @debug("push -f branch to remote")
        push && run(`$git push -q -f -u origin $branch`)

        clean_registry = false
        return RegBranch(pkg.name, pkg.version, branch, nothing)
    finally
        if clean_registry
            @debug("cleaning up possibly inconsistent registry", registry_path=showsafe(registry_path), err=showsafe(err))
            rm(registry_path; recursive=true, force=true)
        end
    end
end

"""
    using <package>
    register(<package>, registry_path)

Register `package` or a new version of `package` in the registry
working copy at `registry_path`. The version number is obtained from
the version field of the package's `Project.toml`.

Note: This will only update the registry locally. Review the result
and `git push` it manually.

    register(<package>, registry_path; commit = false)

Only make the changes to the registry but do not commit.
"""
function register(package::Module, registry_path; repo = nothing, commit = true,
                  gitconfig::Dict = Dict())
    registry_path = expanduser(registry_path)
    if LibGit2.isdirty(LibGit2.GitRepo(registry_path))
        throw("Registry repo is dirty. Stash or commit files.")
    end
    git = gitcmd(registry_path, gitconfig)

    package_path = normpath(joinpath(dirname(pathof(package)), ".."))
    pkg = Pkg.Types.read_project(joinpath(package_path, "Project.toml"))
    tree_hash = get_tree_hash(package_path, gitconfig)
    package_repo = repo === nothing ? get_remote_repo(package_path, gitconfig) : repo

    @info "Registering package" package_path registry_path package_repo uuid=pkg.uuid version=pkg.version tree_hash
    clean_registry = true

    try
        package_path = get_package_path(registry_path, pkg)
        update_package_data(pkg, registry_path, package_repo, package_path, tree_hash)
        commit && commit_registry(pkg, package_path, package_repo, tree_hash, git)
        clean_registry = false
    finally
        if clean_registry
            run(`$git reset --hard`)
        end
    end

    return
end

"""
    create_registry(path, repo; description = nothing)

Create a registry at the local directory `path` and with repository
URL `repo`. Optionally add a description of the purpose of the
registry with the keyword argument `description`.

Note: This will only prepare the registry locally. Review the result
and `git push` it manually.
"""
function create_registry(name, repo; description = nothing, gitconfig::Dict = Dict())
    path = joinpath(pwd(), name)
    isdir(path) && throw("$path already exists.")
    mkpath(path)

    registry_info = Dict()
    registry_info["name"] = name
    registry_info["uuid"] = UUIDs.uuid1()
    registry_info["repo"] = repo
    description !== nothing && (registry_info["description"] = description)
    registry_info["packages"] = Dict{String, Any}()
    write_registry(joinpath(path, "Registry.toml"), registry_info)

    git = gitcmd(path, gitconfig)
    run(`$git init -q`)
    run(`$git add --all`)
    run(`$git commit -qm 'initial commit'`)
    run(`$git remote add origin $repo`)

    return path
end


text:
using LsqFit

include("utils.jl")
function dfa(z; scales=3:0.5:8, offset=-1)
    # cummulative sum of (z - <z>)
    zcum = cumsum(z .- mean(z))

    # define the rulers for each x-window
    scales = collect(scales)
    rulers = convert.(Int64,round.(2 .^ scales))

    # each rulers will have rms averaged for each window in the data
    fluct = zeros(length(rulers))

    for (k, ℓ) in enumerate(rulers)
        # x data is the same for all windows in this scale
        xdata = collect(1:ℓ)

        #each element contail the data in one window
        yvec = strided(zcum,ℓ,offset)

        #for each window one rms
        rms = zeros(length(yvec))

        for (i,ydata) in enumerate(yvec)
            fit = curve_fit((x, p) -> p[1] .+ p[2].*x, xdata, ydata, [1.0, 1.0])
            σ = stderror(fit)
            rms[i] = sqrt(mean(fit.resid.^2))
            # yfit = model(xdata, fit.param)
            # rms2 = sqrt(mean((ydata .- yfit).^2))
        end
        fluct[k] = mean(rms)
    end

    # now fit (fluct vs. scale)
    x = log10.(rulers)
    y = log10.(fluct)
    fit = curve_fit((x, p) -> p[1] .+ p[2].*x ,x, y, [.0, .0])
    h = fit.param[2]-1.0
    return  h,x,y
end

function dfa2d(z; scales=3:0.5:8, offset=-1)
    L = size(z)[1]
    hurst = 0
    n = 0
    for i=10:10:L-10
        h, x, y = dfa(z[:,i]; scales=scales, offset=offset)
        hurst += h
        n+=1
    end
    for i=10:10:L-10
        h, x, y = dfa(z[i,:]; scales=scales, offset=offset)
        hurst += h
        n+=1
    end
    hurst = hurst/n
end

ext: jl


text:
using QuWalkEC, LightGraphs
using Base.Test

@testset "have ss" begin
  @testset "provided subgraph is connected" begin
    @test !have_ss_connected(CompleteGraph(10), [5])
    @test have_ss_connected(CompleteGraph(10), [2, 6, 8])
    @test have_ss_connected(CompleteGraph(10), [2, 6])

    @test have_ss_connected(PathGraph(10), [2, 3])
    @test !have_ss_connected(PathGraph(10), [1, 2])

    @test have_ss_connected(smallgraph(:bull), [1,2,3])
    @test !have_ss_connected(smallgraph(:bull), [1,2])

    @test_throws ArgumentError have_ss_connected(CompleteGraph(10), [2, 2])
    @test_throws ArgumentError have_ss_connected(CompleteGraph(10), [2, 2, 3])
  end

  @testset "generalt subgraph" begin
    @test have_ss(PathGraph(10), [2, 3, 5, 6])
    @test !have_ss(PathGraph(10), [1, 2, 6, 7])

    @test_throws ArgumentError have_ss(PathGraph(10), [2, 2, 6, 7])
    @test_throws ArgumentError have_ss(PathGraph(10), [2, 2, 3])
  end

end

@testset "searching for EC" begin
  @test search_subgraph_ss(CompleteGraph(10)) == [1,2]
  @test search_subgraph_ss(Graph(10)) == Int[]
  @test search_subgraph_ss(PathGraph(10)) == [2, 3]
  @test search_subgraph_ss(StarGraph(10)) == []
end

max_stars_repo_name: QuantumWalks/ExceptionalConfiguration

max_issues_repo_name: QuantumWalks/ExceptionalConfiguration

max_forks_repo_name: QuantumWalks/ExceptionalConfiguration


text:
#!/usr/bin/env julia

using Luxor, Colors

fname = "/tmp/holes.pdf"
Drawing(1200, 1200, fname)

origin()
background("black")

pagetiles = Tiler(1200, 1200, 17, 19, margin=50)
for (pos, n) in pagetiles
  randomhue()
  rad = pagetiles.tilewidth
  ngon(pos, rad/2, rand(5:12), 0, :path)
  newsubpath()
  ngon(pos, rand(1:rad/3), rand(5:12), 0, :path, reversepath=true)
  fillstroke()
end

finish()
println("finished test: output in $(fname)")


text:
"""
Type encoding a rectangular region
"""
struct Rect{T}
  top::T
  left::T
  bottom::T
  right::T

  # Constructor ensures that top is above bottom etc.
  function Rect{T}(t::T, l::T, b::T, r::T) where T
    top = t > b ? b : t
    bottom = t > b ? t : b
    left = l > r ? r : l
    right = l > r ? l : r
    new(top, left, bottom, right)
  end
end
Rect(t::T, l::T, b::T, r::T) where T = Rect{T}(t, l, b, r)

top(r::Rect) = r.top
left(r::Rect) = r.left
bottom(r::Rect) = r.bottom
right(r::Rect) = r.right
width(r::Rect) = r.right - r.left
height(r::Rect) = r.bottom - r.top

"""
isoverlapping(a::Rect, b::Rect)

Return true if two rectangles are overlapping.
Warning: If one shape completely contains the other one this will not work!
"""
function isoverlapping(a::Rect, b::Rect)
  return !(right(a) < left(b) ||
           right(b) < left(a) ||
           bottom(a)< top(b) ||
           bottom(b)< top(a))
end

"""
randrect(maxx::Integer, maxy::Integer, width::Integer, height::Integer)

Return a random rectangle given size and max coordinates
"""
function randrect(maxx::Integer, maxy::Integer,
                  width::Integer, height::Integer)
  top = rand(1:(maxy-height))
  left = rand(1:(maxx-width))
  return Rect(top, left, top + height, left + width)
end

"""
Error that is thrown when the search for non-overlapping patches reaches
an upper limit of tries
"""
struct ExhaustionError{T}
  message::T
end

"""
nonoverlappingrandrect(maxx, maxy, width, height, amount; maxtries = 100)

Return a Vector of rectangles that are all non-overlapping, given size and
max coordinates.

maxtries defines how many times we want to try to find a non-overlapping
  orientation. After this we get an ExhaustionError
"""
function nonoverlappingrandrect(maxx::Integer, maxy::Integer,
                                width::Integer, height::Integer,
                                amount::Integer; maxtries = 100)
  rects = Vector{Rect}()
  for i in 1:amount
    for tries in 1:maxtries
      newrect = randrect(maxx, maxy, width, height)

      if !any(x -> isoverlapping(newrect, x), rects)
        push!(rects, newrect)
        break
      end

      tries == maxtries && throw(ExhaustionError("$max_tries was not enough"))
    end
  end
  return rects
end


text:
#!/usr/bin/env julia

# MPC Command Publisher/Controller Module Interface to the Genesis.
# This version using the Nonlinear Kinematic Bicycle Model.

###########################################
#### ROBOTOS
###########################################
using RobotOS
@rosimport genesis_path_follower.msg: state_est
@rosimport genesis_path_follower.msg: mpc_path
@rosimport std_msgs.msg: UInt8
@rosimport std_msgs.msg: Float32
rostypegen()
using genesis_path_follower.msg
using std_msgs.msg
using PyCall

###########################################
#### LOAD ROSPARAMS
###########################################
if has_param("mat_waypoints")
	mat_fname = get_param("mat_waypoints")
else
	error("No Matfile of waypoints provided!")
end

if has_param("track_using_time") && has_param("target_vel")
	track_with_time = get_param("track_using_time")
	target_vel = get_param("target_vel")
else
	error("Invalid rosparam trajectory definition: track_using_time and target_vel")
end

if has_param("scripts_dir")
	scripts_dir = get_param("scripts_dir")
else
	error("Did not provide the scripts directory!")
end

if(has_param("lat0") && has_param("lat0") && has_param("lat0"))
	lat0 = get_param("lat0")
	lon0 = get_param("lon0")
	yaw0 = get_param("yaw0")
else
	error("Invalid rosparam global origin provided!")
end

###########################################
#### Reference GPS Trajectory Module
###########################################
# Access Python modules for path processing.  Ugly way of doing it, can seek to clean this up in the future.
using PyCall
const gps_utils_loc = scripts_dir * "gps_utils"
unshift!(PyVector(pyimport("sys")["path"]), gps_utils_loc) # append the current directory to Python path
@pyimport ref_gps_traj as rgt
grt = rgt.GPSRefTrajectory(mat_filename=mat_fname, LAT0=lat0, LON0=lon0, YAW0=yaw0)

###########################################
#### MPC Controller Module with Cost Function Weights.
#### Global Variables for Callbacks/Control Loop.
###########################################
push!(LOAD_PATH, scripts_dir * "mpc_utils")
import GPSKinMPCPathFollower
const kmpc = GPSKinMPCPathFollower
kmpc.update_cost(9.0, 9.0, 10.0, 0.0, 100.0, 1000.0, 0.0, 0.0) # x,y,psi,v,da,ddf,a,df

# Reference for MPC
const t_ref = collect(0:kmpc.dt:kmpc.N*kmpc.dt)
x_ref = zeros(length(t_ref))
y_ref = zeros(length(t_ref))
psi_ref = zeros(length(t_ref))

if target_vel > 0.0
	des_speed = target_vel
else
	des_speed = 0.00
end

ref_lock = false
received_reference = false
x_curr  = 0.0
y_curr  = 0.0
psi_curr  = 0.0
v_curr  = 0.0
command_stop = false

###########################################
#### State Estimation Callback.
###########################################
function state_est_callback(msg::state_est)

	global x_curr, y_curr, psi_curr, v_curr
	global received_reference

	if ref_lock == false
		x_curr = msg.x
		y_curr = msg.y
		psi_curr = msg.psi
		v_curr = msg.v
		received_reference = true
	end
end

function pub_loop(acc_pub_obj, steer_pub_obj, mpc_path_pub_obj)
    loop_rate = Rate(20.0)
    while ! is_shutdown()
	    if ! received_reference		# Reference not received so don't use MPC yet.
	        rossleep(loop_rate)
	        continue
	    end

	    global ref_lock				# Ref lock used to ensure that get/set of state doesn't happen simultaneously.
	    ref_lock = true

		global x_curr, y_curr, psi_curr, v_curr, des_speed, command_stop

		if ! track_with_time
			# fixed velocity-based path tracking
			x_ref, y_ref, psi_ref, stop_cmd = grt[:get_waypoints](x_curr, y_curr, psi_curr, des_speed)

			if stop_cmd == true
				command_stop = true
			end

		else
			# trajectory tracking
			x_ref, y_ref, psi_ref, stop_cmd = grt[:get_waypoints](x_curr, y_curr, psi_curr)

			if stop_cmd == true
				command_stop = true
			end
		end

		# Update Model
		kmpc.update_init_cond(x_curr, y_curr, psi_curr, v_curr)
		kmpc.update_reference(x_ref, y_ref, psi_ref, des_speed)

	    ref_lock = false

		if command_stop == false
			a_opt, df_opt, is_opt, solv_time = kmpc.solve_model()

			rostm = get_rostime()
			tm_secs = rostm.secs + 1e-9 * rostm.nsecs

		    log_str = @sprintf("Solve Status: %s, Acc: %.3f, SA: %.3f, ST: %.3f", is_opt, a_opt, df_opt, solv_time)
		    loginfo(log_str)

		    if is_opt == :Optimal
				publish( acc_pub_obj,   Float32Msg(a_opt) )
				publish( steer_pub_obj, Float32Msg(df_opt) )
			end

			kmpc.update_current_input(df_opt, a_opt)
			res = kmpc.get_solver_results()

			mpc_path_msg = mpc_path()
			mpc_path_msg.header.stamp = rostm
			mpc_path_msg.solv_status  = string(is_opt)
			mpc_path_msg.solv_time = solv_time
			mpc_path_msg.xs   = res[1] 	# x_mpc
			mpc_path_msg.ys   = res[2] 	# y_mpc
			mpc_path_msg.vs   = res[3]	# v_mpc
			mpc_path_msg.psis = res[4] 	# psi_mpc
			mpc_path_msg.xr   = res[5] 	# x_ref
			mpc_path_msg.yr   = res[6] 	# y_ref
			mpc_path_msg.vr   = [res[7]]	# v_ref
			mpc_path_msg.psir = res[8] 	# psi_ref
			mpc_path_msg.df   = res[9]	# d_f
			mpc_path_msg.acc  = res[10]	# acc
			publish(mpc_path_pub_obj, mpc_path_msg)
		else
			publish( acc_pub_obj,   Float32Msg(-1.0) )
			publish( steer_pub_obj, Float32Msg(0.0) )
		end

	    rossleep(loop_rate)
	end
end

function start_mpc_node()
    init_node("dbw_mpc_pf")
    acc_pub   = Publisher("/control/accel", Float32Msg, queue_size=2)
    steer_pub = Publisher("/control/steer_angle", Float32Msg, queue_size=2)

    acc_enable_pub   = Publisher("/control/enable_accel", UInt8Msg, queue_size=2, latch=true)
    steer_enable_pub = Publisher("/control/enable_spas",  UInt8Msg, queue_size=2, latch=true)

    mpc_path_pub = Publisher("mpc_path", mpc_path, queue_size=2)
	sub_state  = Subscriber("state_est", state_est, state_est_callback, queue_size=2)

	# Start up Ipopt/Solver.
	for i=1:3
		kmpc.solve_model()
	end

	publish(acc_enable_pub, UInt8Msg(2))
	publish(steer_enable_pub, UInt8Msg(1))

    pub_loop(acc_pub, steer_pub, mpc_path_pub)
end

if ! isinteractive()
	try
    	start_mpc_node()
    catch x
    	print(x)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "teem"
version = v"0.0.0"

# Collection of sources required to build teem
sources = [
    "https://downloads.sourceforge.net/project/teem/teem/1.11.0/teem-1.11.0-src.tar.gz?r=https%3A%2F%2Fsourceforge.net%2Fprojects%2Fteem%2Ffiles%2Fteem%2F1.11.0%2Fteem-1.11.0-src.tar.gz%2Fdownload%3Fuse_mirror%3Dcfhcable%26download%3D&ts=1565376769" =>
    "a01386021dfa802b3e7b4defced2f3c8235860d500c1fa2f347483775d4c8def",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mv teem-1.11.0-src.tar.gz\?r\=https%3A%2F%2Fsourceforge.net%2Fprojects%2Fteem%2Ffiles%2Fteem%2F1.11.0%2Fteem-1.11.0-src.tar.gz%2Fdownload%3Fuse_mirror%3Dcfhcable%26download%3D\&ts\=1565376769 teem.tar.gz
tar -xvf teem.tar.gz
mkdir teem-build
cd teem-build/
cmake ../teem-1.11.0-src/
cmake  -DBUILD_SHARED_LIBS=True -DCMAKE_INSTALL_PREFIX=/workspace/destdir ../teem-1.11.0-src/
make
make install

if [ $target = "x86_64-apple-darwin14" ]; then
cd $WORKSPACE/srcdir
cd teem-build
cmake ../teem-1.11.0-src/
cmake  -DBUILD_SHARED_LIBS=True -DCMAKE_INSTALL_PREFIX=/workspace/destdir ../teem-1.11.0-src/
make
make install
fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    Linux(:i686, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libteem", :libteem)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
import Base.getindex

"""
    G = LQG(A,B,C,D, Q1, Q2, R1, R2; qQ=0, qR=0, integrator=false)
    G = LQG(sys, args...; kwargs...)

Return an LQG object that describes the closed control loop around the process `sys=ss(A,B,C,D)`
where the controller is of LQG-type. The controller is specified by weight matrices `R1,R2`
that penalizes state deviations and control signal variance respectively, and covariance
matrices `Q1,Q2` which specify state drift and measurement covariance respectively.
This constructor calls [`lqr`](@ref) and [`kalman`](@ref) and forms the closed-loop system.

If `integrator=true`, the resulting controller will have intregral action.
This is achieved by adding a model of a constant disturbance on the inputs to the system
described by `A,B,C,D`.

`qQ` and `qR` can be set to incorporate loop transfer recovery, i.e.,
```julia
L = lqr(A, B, Q1+qQ*C'C, Q2)
K = kalman(A, C, R1+qR*B*B', R2)
```

# Fields
When the LQG-object is populated by the lqg-function, the following fields have been made available
- `L` is the feedback matrix, such that `A-BL` is stable. Note that the length of the state vector (and the width of L) is increased by the number of inputs if the option `integrator=true`.
- `K` is the kalman gain such that `A-KC` is stable
- `sysc` is a dynamical system describing the controller `u=L*inv(A-BL-KC+KDL)Ky`

# Functions
Several other properties of the object are accessible with the indexing function `getindex()`
and are called with the syntax `G[:function]`. The available functions are
(some have many alternative names, separated with / )

-`G[:cl] / G[:closedloop]` is the closed-loop system, including observer, from reference to output, precompensated to have static gain 1 (`u = −Lx + lᵣr`).
-`G[:S] / G[:Sin]` Input sensitivity function
-`G[:T] / G[:Tin]` Input complementary sensitivity function
-`G[:Sout]` Output sensitivity function
-`G[:Tout]` Output complementary sensitivity function
-`G[:CS]` The transfer function from measurement noise to control signal
-`G[:DS]` The transfer function from input load disturbance to output
-`G[:lt] / G[:looptransfer] / G[:loopgain]  =  PC`
-`G[:rd] / G[:returndifference]  =  I + PC`
-`G[:sr] / G[:stabilityrobustness]  =  I + inv(PC)`
-`G[:sysc] / G[:controller]` Returns the controller as a StateSpace-system

It is also possible to access all fileds using the `G[:symbol]` syntax, the fields are `P
,Q1,Q2,R1,R2,qQ,qR,sysc,L,K,integrator`

# Example

```julia
eye(n) = Matrix{Float64}(I,n,n) # For convinience

qQ = 1
qR = 1
Q1 = 10eye(4)
Q2 = 1eye(2)
R1 = 1eye(6)
R2 = 1eye(2)

G = LQG(sys, Q1, Q2, R1, R2, qQ=qQ, qR=qR, integrator=true)

Gcl = G[:cl]
T = G[:T]
S = G[:S]
sigmaplot([S,T],exp10.(range(-3, stop=3, length=1000)))
stepplot(Gcl)
```

"""
struct LQG
    P::StateSpace
    Q1::AbstractMatrix
    Q2::AbstractMatrix
    R1::AbstractMatrix
    R2::AbstractMatrix
    qQ::Real
    qR::Real
    sysc::LTISystem
    L::AbstractMatrix
    K::AbstractMatrix
    integrator::Bool
end

# Provide some constructors
function LQG(A,B,C,D,Q1::AbstractMatrix,Q2::AbstractMatrix,R1::AbstractMatrix,R2::AbstractMatrix; qQ=0, qR=0, integrator=false)
    integrator ? _LQGi(A,B,C,D,Q1,Q2,R1,R2, qQ, qR) : _LQG(A,B,C,D,Q1,Q2,R1,R2, qQ, qR)
end # (1) Dispatches to final

function LQG(A,B,C,D,Q1::AbstractVector,Q2::AbstractVector,R1::AbstractVector,R2::AbstractVector; qQ=0, qR=0, integrator=false)
    Q1 = diagm(0 => Q1)
    Q2 = diagm(0 => Q2)
    R1 = diagm(0 => R1)
    R2 = diagm(0 => R2)
    integrator ? _LQGi(A,B,C,D,Q1,Q2,R1,R2, qQ, qR) : _LQG(A,B,C,D,Q1,Q2,R1,R2, qQ, qR)
end # (2) Dispatches to final

# (3) For conveniece of sending a sys, dispatches to (1/2)
LQG(sys::LTISystem, args...; kwargs...) = LQG(sys.A,sys.B,sys.C,sys.D,args...; kwargs...)

# This function does the actual initialization in the standard case withput integrator
function _LQG(A,B,C,D, Q1, Q2, R1, R2, qQ, qR)
    n = size(A,1)
    m = size(B,2)
    p = size(C,1)
    L = lqr(A, B, Q1+qQ*C'C, Q2)
    K = kalman(A, C, R1+qR*B*B', R2)

    # Controller system
    Ac=A-B*L-K*C+K*D*L
    Bc=K
    Cc=L
    Dc=zero(D')
    sysc = ss(Ac,Bc,Cc,Dc)

    return LQG(ss(A,B,C,D),Q1,Q2,R1,R2, qQ, qR, sysc, L, K, false)
end

# This function does the actual initialization in the integrator case
function _LQGi(A,B,C,D, Q1, Q2, R1, R2, qQ, qR)
    n = size(A,1)
    m = size(B,2)
    p = size(C,1)

    # Augment with disturbance model
    Ae = [A B; zeros(m,n+m)]
    Be = [B;zeros(m,m)]
    Ce = [C zeros(p,m)]
    De = D

    L = lqr(A, B, Q1+qQ*C'C, Q2)
    Le = [L I]
    K = kalman(Ae, Ce, R1+qR*Be*Be', R2)

    # Controller system
    Ac=Ae-Be*Le-K*Ce+K*De*Le
    Bc=K
    Cc=Le
    Dc=zero(D')
    sysc = ss(Ac,Bc,Cc,Dc)

    LQG(ss(A,B,C,D),Q1,Q2,R1,R2, qQ, qR, sysc, Le, K, true)
end

function Base.getindex(G::LQG, s)
    s == :A && return G.P.A
    s == :B && return G.P.B
    s == :C && return G.P.C
    s == :D && return G.P.D
    s == :L  && return G.L
    s == :K  && return G.K
    s == :Q1 && return G.Q1
    s == :Q2 && return G.Q2
    s == :R1 && return G.R1
    s == :R2 && return G.R2
    s == :qQ && return G.qQ
    s == :qR && return G.qR
    s ∈ [:sys, :P]  && return G.P
    s ∈ [:sysc, :controller] && return G.sysc
    s == :integrator && return G.integrator

    A = G.P.A
    B = G.P.B
    C = G.P.C
    D = G.P.D

    L = G.L
    K = G.K
    P = G.P
    sysc = G.sysc

    n = size(A,1)
    m = size(B,2)
    p = size(C,1)

    # Extract interesting values
    if G.integrator # Augment with disturbance model
        A = [A B; zeros(m,n+m)]
        B = [B;zeros(m,m)]
        C = [C zeros(p,m)]
        D = D
    end

    PC = P*sysc # Loop gain

    if s ∈ [:cl, :closedloop, :ry] # Closed-loop system
        Acl = [A-B*L B*L; zero(A) A-K*C]
        Bcl = [B; zero(B)]
        Ccl = [C zero(C)]
        Bcl = Bcl/(P.C*inv(P.B*L[:,1:n]-P.A)*P.B) # B*lᵣ # Always normalized with nominal plant static gain
        return syscl = ss(Acl,Bcl,Ccl,0)
    elseif s ∈ [:Sin, :S] # Sensitivity function
        return feedback(ss(Matrix{numeric_type(PC)}(I, m, m)),PC)
    elseif s ∈ [:Tin, :T] # Complementary sensitivity function
        return feedback(PC)
    elseif s == :Sout # Sensitivity function, output
        return feedback(ss(Matrix{numeric_type(sys_c)}(I, m, m)),sysc*P)
    elseif s == :Tout # Complementary sensitivity function, output
        return feedback(sysc*P)
    elseif s == :PS # Load disturbance to output
        return P*G[:S]
    elseif s == :CS # Noise to control signal
        return sysc*G[:S]
    elseif s ∈ [:lt, :looptransfer, :loopgain]
        return PC
    elseif s ∈ [:rd, :returndifference]
        return  ss(Matrix{numeric_type(PC)}(I, p, p)) + PC
    elseif s ∈ [:sr, :stabilityrobustness]
        return  ss(Matrix{numeric_type(PC)}(I, p, p)) + inv(PC)
    end
    error("The symbol $s does not have a function associated with it.")
end

Base.:(==)(G1::LQG, G2::LQG) = G1.K == G2.K && G1.L == G2.L && G1.P == G2.P && G1.sysc == G2.sysc

plot(G::LQG) = gangoffourplot(G)
function gangoffour(G::LQG)
    return G[:S], G[:PS], G[:CS], G[:T]
end

function gangoffourplot(G::LQG; kwargs...)
    S,D,N,T = gangoffour(G)
    f1 = sigmaplot(S, show=false, kwargs...); Plots.plot!(title="\$S = 1/(1+PC)\$")
    f2 = sigmaplot(D, show=false, kwargs...); Plots.plot!(title="\$D = P/(1+PC)\$")
    f3 = sigmaplot(N, show=false, kwargs...); Plots.plot!(title="\$N = C/(1+PC)\$")
    f4 = sigmaplot(T, show=false, kwargs...); Plots.plot!(title="\$T = PC/(1+PC\$)")
    Plots.plot(f1,f2,f3,f4)
end

# function gangoffourplot(G::LQG, args...)
#     S,D,N,T = gangoffour(G)
#     fig = subplot(n=4,nc=2)
#     Plots.plot!(fig[1,1],sigmaplot(S, args...), title="\$S = 1/(1+PC)\$")
#     Plots.plot!(fig[1,2],sigmaplot(D, args...), title="\$D = P/(1+PC)\$")
#     Plots.plot!(fig[2,1],sigmaplot(N, args...), title="\$N = C/(1+PC)\$")
#     Plots.plot!(fig[2,2],sigmaplot(T, args...), title="\$T = PC/(1+PC\$)")
#     return fig
# end


text:
#!/usr/bin/env julia --project

include("../shared_source/boilerplate.jl")
include("ThreeDimensionalCompressibleNavierStokesEquations.jl")

ClimateMachine.init()

#################
# RUN THE TESTS #
#################
@testset "$(@__FILE__)" begin

    include("refvals_bickley_jet.jl")

    ########
    # Setup physical and numerical domains
    ########
    Ωˣ = IntervalDomain(-2π, 2π, periodic = true)
    Ωʸ = IntervalDomain(-2π, 2π, periodic = true)
    Ωᶻ = IntervalDomain(-2π, 2π, periodic = true)

    first_order = DiscretizedDomain(
        Ωˣ × Ωʸ × Ωᶻ;
        elements = 32,
        polynomial_order = 1,
        overintegration_order = 1,
    )

    fourth_order = DiscretizedDomain(
        Ωˣ × Ωʸ × Ωᶻ;
        elements = 13,
        polynomial_order = 4,
        overintegration_order = 1,
    )

    grids = Dict("first_order" => first_order, "fourth_order" => fourth_order)

    ########
    # Define timestepping parameters
    ########
    start_time = 0
    end_time = 100.0
    Δt = 0.004
    method = SSPRK22Heuns

    timestepper = TimeStepper(method = method, timestep = Δt)

    callbacks = (Info(), StateCheck(10))

    ########
    # Define physical parameters and parameterizations
    ########
    parameters = (
        ϵ = 0.1,  # perturbation size for initial condition
        l = 0.5, # Gaussian width
        k = 0.5, # Sinusoidal wavenumber
        ρₒ = 1, # reference density
        cₛ = sqrt(10), # sound speed
    )

    physics = FluidPhysics(;
        advection = NonLinearAdvectionTerm(),
        dissipation = ConstantViscosity{Float64}(μ = 0, ν = 0, κ = 0),
        coriolis = nothing,
        buoyancy = nothing,
    )

    ########
    # Define initial conditions
    ########

    # The Bickley jet
    U₀(p, x, y, z) = cosh(y)^(-2)
    V₀(p, x, y, z) = 0
    W₀(p, x, y, z) = 0

    # Slightly off-center vortical perturbations
    Ψ₁(p, x, y, z) =
        exp(-(y + p.l / 10)^2 / (2 * (p.l^2))) * cos(p.k * x) * cos(p.k * y)
    Ψ₂(p, x, y, z) =
        exp(-(z + p.l / 10)^2 / (2 * (p.l^2))) * cos(p.k * y) * cos(p.k * z)

    # Vortical velocity fields (u, v, w) = (-∂ʸ, +∂ˣ, 0) Ψ₁ + (0, -∂ᶻ, +∂ʸ)Ψ₂
    u₀(p, x, y, z) =
        Ψ₁(p, x, y, z) * (p.k * tan(p.k * y) + y / (p.l^2) + 1 / (10 * p.l))
    v₀(p, x, y, z) =
        Ψ₂(p, x, y, z) * (p.k * tan(p.k * z) + z / (p.l^2) + 1 / (10 * p.l)) -
        Ψ₁(p, x, y, z) * p.k * tan(p.k * x)
    w₀(p, x, y, z) = -Ψ₂(p, x, y, z) * p.k * tan(p.k * y)
    θ₀(p, x, y, z) = sin(p.k * y)

    ρ₀(p, x, y, z) = p.ρₒ
    ρu₀(p, x...) = ρ₀(p, x...) * (p.ϵ * u₀(p, x...) + U₀(p, x...))
    ρv₀(p, x...) = ρ₀(p, x...) * (p.ϵ * v₀(p, x...) + V₀(p, x...))
    ρw₀(p, x...) = ρ₀(p, x...) * (p.ϵ * w₀(p, x...) + W₀(p, x...))
    ρθ₀(p, x...) = ρ₀(p, x...) * θ₀(p, x...)

    ρu⃗₀(p, x...) = @SVector [ρu₀(p, x...), ρv₀(p, x...), ρw₀(p, x...)]
    initial_conditions = (ρ = ρ₀, ρu = ρu⃗₀, ρθ = ρθ₀)

    for (key, grid) in grids
        @testset "$(key)" begin
            model = SpatialModel(
                balance_law = Fluid3D(),
                physics = physics,
                numerics = (flux = RoeNumericalFlux(),),
                grid = grid,
                boundary_conditions = NamedTuple(),
                parameters = parameters,
            )

            simulation = Simulation(
                model = model,
                initial_conditions = initial_conditions,
                timestepper = timestepper,
                callbacks = callbacks,
                time = (; start = start_time, finish = end_time),
            )

            ########
            # Run the model
            ########
            evolve!(
                simulation,
                model;
                refDat = getproperty(refVals, Symbol(key)),
            )
        end
    end
end


text:
using Fermi.DIIS

using TensorOperations
using LinearAlgebra
using Formatting
import Base: show

export UHF

abstract type UHFAlgorithm end

struct UHF <: AbstractHFWavefunction
    molecule::Molecule
    energy::Float64
    nocc::Int
    nvir::Int
    orbitals::UHFOrbitals
    e_conv::Float64
    d_conv::Float64
end

function UHF(x...)
    if !any(i-> i isa UHFAlgorithm, x)
        UHF(x..., get_uhf_alg())
    else
        # Print the type of arguments given for a better feedback
        args = "("
        for a in x[1:end-1]
            args *= "$(typeof(a)), "
        end
        args = args[1:end-2]*")"
        throw(FermiException("invalid arguments for UHF method: $args"))
    end
end

"""
    Fermi.HartreeFock.get_uhf_alg()

Returns a singleton type corresponding to a UHF implementation.
"""
function get_uhf_alg(N::Int = Options.get("uhf_alg"))
    try
        return get_uhf_alg(Val(N))
    catch MethodError
        throw(FermiException("implementation number $N not available for UHF."))
    end
end

struct UHFa <: UHFAlgorithm end
include("UHFa.jl")
include("UHFHelper.jl")
# And a number is assigned to the implementation
get_uhf_alg(x::Val{1}) = UHFa()

function string_repr(X::UHF)
    out = ""
    out = out*" ⇒ Fermi Unrestricted Hartree--Fock Wave function\n"
    out = out*" ⋅ Basis:                  $(X.orbitals.basis)\n"
    out = out*" ⋅ Energy:                 $(X.energy)\n"
    out = out*" ⋅ Occ. Spatial Orbitals:  $(X.nocc)\n"
    out = out*" ⋅ Vir. Spatial Orbitals:  $(X.nvir)\n"
    out = out*"Convergence: "
    out = out*"ΔE => $(format("{:1.2e}",abs(X.e_conv)))"
    out = out*" Dᵣₘₛ => $(format("{:1.2e}",abs(X.d_conv)))"
    return out
end

function show(io::IO, ::MIME"text/plain", X::UHF)
    print(string_repr(X))
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
for test_file in readdir(joinpath(dirname(@__FILE__), "..", "examples"))
    if test_file[end-2:end] == ".jl"
        include(joinpath(dirname(@__FILE__), "..", "examples", test_file))
    end
end


text:
# A RulePtr tracks a location within a rule
#==
 The index value has a pythonic meaning of
 being "between" two values. So an index of
 0 means 'prior to the first element'
 and an index of length = 'after the last
 element' (and therefore satisfied). So
 in Python the rule starts at zero and
 finishes at length (even though the
 expansion only goes to coordinate length-1)

 We emulate this in Julia by interpreting
 the index accordingly as zero-based.

==#
struct RulePtr
    rule::Rule
    index::Int
    RulePtr(r::Rule,i) = begin
        @assert i <= length(r.expansion) "$i beyond end of $r"
        new(r,i)
        end
end

Base.show(io::IO,r::RulePtr) = begin
    if r.index == 0
        before = []
        after = r.rule.expansion
    elseif r.index < length(r.rule.expansion) #emulate python
        before = r.rule.expansion[1:r.index]
        after = r.rule.expansion[r.index+1:end]
    else
        before = r.rule.expansion
        after = []
    end
    print(io,"<$(r.rule.origin.name) : $(join(before," ")) * $(join(after," "))>")
end

# For sorted displays
Base.isless(rp1::RulePtr,rp2::RulePtr) = begin
    rp1.rule.origin.name < rp2.rule.origin.name ||
        rp1.rule.origin.name == rp2.rule.origin.name && rp1.index < rp2.index
end

next(r::RulePtr) = r.rule.expansion[r.index+1]

advance(r::RulePtr,sym) = begin
    @assert next(r) == sym
    RulePtr(r.rule,r.index + 1)
end

# Have we satisfied the rule by getting to the end?
is_satisfied(r::RulePtr) = r.index == length(r.rule.expansion)

Base.isequal(r1::RulePtr,r2::RulePtr) = r1.rule == r2.rule && r1.index == r2.index
Base.hash(r1::RulePtr,h::UInt64) = hash(r1.rule,hash(r1.index,h))

# state generation ensures no duplicate LR0ItemSets
struct LR0ItemSet
    kernel::Set{RulePtr}
    closure::Set{RulePtr}
    transitions::Dict
    lookaheads::DefaultDict
end

LR0ItemSet(kernel,closure) = LR0ItemSet(Set(kernel),Set(closure),Dict(),DefaultDict{Any,Set}())

Base.show(io::IO,l::LR0ItemSet) = begin
    joinstr = ", "
    println(io,"{$(join(l.kernel,joinstr)) | $(join(l.closure,joinstr))}")
end

# For displaying sorted lists of states
Base.isless(lr1::LR0ItemSet,lr2::LR0ItemSet) = begin
    first(sort(collect(lr1.kernel))) < first(sort(collect(lr2.kernel)))
end

update_set(set1,set2) = begin
    if isempty(set2) || length(set1) > length(set2)
        return false
    end
    old_len = length(set1)
    union!(set1,set2)
    return old_len < length(set1)
end
"""Calculate FOLLOW sets.

Adapted from: http://lara.epfl.ch/w/cc09:algorithm_for_first_and_follow_sets

"""

calculate_sets(rules) = begin
    symbols = union(Set(sym for rule in rules for sym in rule.expansion), Set(rule.origin for rule in rules))
    # foreach grammar rule X ::= Y(1) ... Y(k)
    # if k=0 or {Y(1),...,Y(k)} subset of NULLABLE then
    #   NULLABLE = NULLABLE union {X}
    # for i = 1 to k
    #   if i=1 or {Y(1),...,Y(i-1)} subset of NULLABLE then
    #     FIRST(X) = FIRST(X) union FIRST(Y(i))
    #   for j = i+1 to k
    #     if i=k or {Y(i+1),...Y(k)} subset of NULLABLE then
    #       FOLLOW(Y(i)) = FOLLOW(Y(i)) union FOLLOW(X)
    #     if i+1=j or {Y(i+1),...,Y(j-1)} subset of NULLABLE then
    #       FOLLOW(Y(i)) = FOLLOW(Y(i)) union FIRST(Y(j))
    # until none of NULLABLE,FIRST,FOLLOW changed in last iteration
    NULLABLE = Set{LarkSymbol}()
    FIRST = Dict{LarkSymbol,Set{LarkSymbol}}()
    FOLLOW = Dict{LarkSymbol,Set{LarkSymbol}}()
    for sym in symbols
        FIRST[sym] = if is_terminal(sym) Set([sym]) else Set() end
        FOLLOW[sym] = Set()
    end

    # Calculate NULLABLE and FIRST
    changed = true
    while changed
        changed = false
        for rule in rules
            if Set(rule.expansion) <= NULLABLE
                if update_set(NULLABLE,Set([rule.origin]))
                    changed = true
                end
            end
            for (i,sym) in enumerate(rule.expansion)
                if i==1 || Set(rule.expansion[1:i-1]) <= NULLABLE
                    if update_set(FIRST[rule.origin],FIRST[sym])
                        changed = true
                    end
                elseif i > 1 break end
            end
        end
    end

    # Calculate FOLLOW
    changed = true
    while changed
        changed = false

        for rule in rules
            for (i,sym) in enumerate(rule.expansion)
                if i == length(rule.expansion) || Set(rule.expansion[i+1:end]) <= NULLABLE
                    if update_set(FOLLOW[sym], FOLLOW[rule.origin])
                        changed = true
                    end
                end

                for j in i+1:length(rule.expansion)
                    if j==i+1 || Set(rule.expansion[i+1:j-1]) <= NULLABLE
                        if update_set(FOLLOW[sym], FIRST[rule.expansion[j]])
                            changed = true
                        end
                    end
                end
            end
        end
    end

    #==
    println("\nFIRST: $(length(FIRST)) items")
    for (i,j) in FIRST
        println("$i:$j")
    end
    println("\nFOLLOW: $(length(FOLLOW)) items")
    for (i,j) in FOLLOW
        println("$i:$j")
    end
    println("\nNULLABLE: $NULLABLE") ==#
    return FIRST,FOLLOW, NULLABLE
end

abstract type GrammarAnalyzer end

init_analyser!(g::GrammarAnalyzer,parser_conf;debug=false) = begin
    g.debug = debug
    root_rules = Dict([start => Rule(NonTerminal("\$root_" * start), [NonTerminal(start), Terminal("\$END")]) for start in parser_conf.start])
    rules = vcat(parser_conf.rules, collect(values(root_rules)))
    g.rules_by_origin = classify(rules, key = r -> r.origin)

    if length(rules) != length(Set(rules))
        duplicates = filter(x->count(y->x==y,rules)>1,rules)
        throw(GrammarError("Rules defined twice: $duplicates"))
    end

    for r in rules
        for sym in r.expansion
            if !(is_terminal(sym) || sym in keys(g.rules_by_origin))
                throw(GrammarError("Using an undefined rule: $sym"))
            end
        end
    end

    g.start_states = Dict([start=>expand_rule(g,root_rule.origin) for (start,root_rule) in root_rules])
    g.end_states = Dict([start=>Set([RulePtr(root_rule,length(root_rule.expansion))]) for (start,root_rule) in root_rules])
    lr0_root_rules = Dict([start=>Rule(NonTerminal("\$root_" * start),[NonTerminal(start)])
                           for start in parser_conf.start])
    lr0_rules = vcat(parser_conf.rules,collect(values(lr0_root_rules)))
    @assert length(lr0_rules) == length(unique(lr0_rules))
    g.lr0_rules_by_origin = classify(lr0_rules,key=r->r.origin)

    # Cache RulePtr(r,0) in r (no duplicate RulePtr objects)
    g.lr0_start_states = Dict([start=>LR0ItemSet([RulePtr(root_rule,0)],expand_rule(g,root_rule.origin,g.lr0_rules_by_origin))
                               for (start,root_rule) in lr0_root_rules])
    g.FIRST, g.FOLLOW, g.NULLABLE = calculate_sets(rules)
    return g
end

"""Returns all init_ptrs accessible by source_rule (recursive)"""
expand_rule(g::GrammarAnalyzer,source_rule, rules_by_origin=nothing) = begin
    if rules_by_origin === nothing
        rules_by_origin = g.rules_by_origin
    end

    init_ptrs = Set{RulePtr}()
    _expand_rule(_rule) = begin
        Channel() do rule_chan
            #@assert !_rule.is_term _rule
            for r in rules_by_origin[_rule]
                init_ptr = RulePtr(r,0)
                push!(init_ptrs,init_ptr)

                if !isempty(r.expansion)
                    # Next symbol from this
                    new_r = next(init_ptr)
                    if !is_terminal(new_r)
                        # Non-terminal, send it back
                        put!(rule_chan,new_r)
                    end
                end
            end
        end
    end

    _ = collect(bfs([source_rule], _expand_rule))

    return init_ptrs
end


text:
using Discord
using JSON
include("base.jl")

#management_dict = Dict()

#@show management_dict
function set_channel_id(client::Client,m::Discord.Message,channel_name,channel_mention)
    global management_dict
    if "<#" != channel_mention[1:2]
        return reply(client,m,"`$(channel_mention)` は適切に入力されていません．リンク付きのチャンネルを入力してください")
    end
    channel_mention_id = parse(Int,channel_mention[3:length(channel_mention)-1])
    input_channel_id = string(m.channel_id)

    if !haskey(management_dict,input_channel_id)
        management_dict[input_channel_id] = Dict("guild_id" => m.guild_id)
    end

    if channel_name  == "残凸把握板"
        management_dict[input_channel_id]["remain_attack_channel"] = channel_mention_id
        save_management_dict(management_dict)
        return reply(client,m,"残凸把握版を$(channel_mention)に登録しました")

    elseif channel_name == "タスキル把握板"
        management_dict[input_channel_id]["task_kill_channel"] = channel_mention_id
        save_management_dict(management_dict)
        return reply(client,m,"タスキル把握版を$(channel_mention)に登録しました")

    elseif channel_name == "予約確認板"
        management_dict[input_channel_id]["output_channel"] = channel_mention_id
        save_management_dict(management_dict)
        return reply(client,m,"予約確認版を$(channel_mention)に登録しました")
    else
        return reply(client,m,"チャンネルの種類が適切に入力されておりません(`$(channel_name)`)．`残凸把握板` `予約確認板`のいずれかを入力してください")
    end
end

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle


text:
abstract Decoder
abstract DeterministicDecoder <: Decoder
abstract QuasiRandomDecoder <: Decoder
abstract RandomDecoder <: Decoder
# abstract EMDecoder <: QuasiRandomDecoder
# abstract GibbsDecoder <: RandomDecoder
# abstract ParticleFilterDecoder <: RandomDecoder

type CorDecoder <: DeterministicDecoder
end

abstract MCMLDecoder <: QuasiRandomDecoder

function decode{T<:Real}(d::CorDecoder, x::Vector{T}, y::Vector{T})
  dec = dot(x, y)
  if dec > 0
    return 1
  elseif dec < 0
    return -1
  elseif dec == 0
    return -2
  end
end


text:
using ServiceMeshLearning, Flux

hyperparams = Dict(pairs((;
    # Set params
    tag = "test_simple",
    env = :single,

    # Optimized params
    lr_alpha = 5.0e-5,
    target_entropy = 7.5,
    frames = 2,
    γ = 0.93,
    τ = 0.01,
    start_steps = 50000,
    update_after = 1000,
    replay_size = 500000,
    hidden_units_policy = 20,
    hidden_units_value = 120,
    hidden_layers_policy = 2,
    hidden_layers_value = 4,
    batch_size = 100,
    update_freq = 8,
    actfun = elu,
    reward_scaling = 0.5,
    closing_cost = 0.0,

    # Always
    logging = true,
    verbose = true,
    seed = 24837293,
    timesteps = 20_000_000,
    start_policy = :RandomAgent,
)))

# Run RL + baselines locally with logging and different seed

for agent in [:SAC, :SimpleAgent, :K8sAgent, :OracleAgent]
    run_experiment(;
        hyperparams...,
        alg=agent,
    )
end

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
struct Vector3D{T} <: FieldVector{3, T}
    x̂::T
    ŷ::T
    ẑ::T
end

struct PseudoVector3D{T} <: FieldVector{3, T}
    x̂ŷ::T
    x̂ẑ::T
    ŷẑ::T
end

struct MultiVector3D{T} <: FieldVector{8, T}
    ŝ::T
    x̂::T
    ŷ::T
    ẑ::T
    x̂ŷ::T
    x̂ẑ::T
    ŷẑ::T
    x̂ŷẑ::T
end

Base.:*(s::T, ::Scalar) where {T <: Real}         = MultiVector3D(s, zeros(T, 7)...)
Base.:*(s::Real, b::VectorBasis3D)                =       Vector3D((i ≠ b ? 0 : s for i ∈ instances(VectorBasis3D))...)
Base.:*(s::Real, b::PseudoVectorBasis3D)          = PseudoVector3D((i ≠ b ? 0 : s for i ∈ instances(PseudoVectorBasis3D))...)
Base.:*(s::T, ::PseudoScalar3D) where {T <: Real} = MultiVector3D(zeros(T, 7)..., s)

for (Type, basis) ∈ [:Vector3D => instances(VectorBasis3D), :PseudoVector3D => instances(PseudoVectorBasis3D), :MultiVector3D => _uvecs]
    quote
        show(io::IO, ::MIME"text/plain", a::$Type) = show(io, a)
        function show(io::IO, a::$Type)
            compact = get(io, :compact, false)
            vec = [zip(a, $basis)...][map(!iszero, a)]
            isempty(vec) && return
            coord, uv = vec[begin]
            show(io, coord); show(io, uv)

            for (coord, uv) ∈ vec[begin+1:end]
                iszero(coord) && continue
                if signbit(coord) && !isnan(coord)
                    print(io, compact ? "-" : " - ")
                    if isa(coord,Signed) && !isa(coord,BigInt) && coord == typemin(typeof(coord))
                        show(io, -widen(coord))
                    else
                        show(io, -coord)
                    end
                else
                    print(io, compact ? "+" : " + ")
                    show(io, coord)
                end
                if !(isa(coord,Integer) && !isa(coord,Bool) || isa(coord,AbstractFloat) && isfinite(coord))
                    print(io, "*")
                end
                print(io, uv)
            end
        end
    end |> eval
end

@inline function _apply(f::Function, mv::MultiVector3D, v::Union{Vector3D, PseudoVector3D}, range::AbstractVector)
    arr = mv |> Vector
    arr[range] = f(@view(arr[range]), v)
    arr |> MultiVector3D
end

@inline function _apply(f::Function, mv::MultiVector3D, z::Complex, range::AbstractVector)
    arr = mv |> Vector
    arr[range] = f(@view(arr[range]), reim(z) |> Vector)
    arr |> MultiVector3D
end

@inline function _apply(f::Function, mv::MultiVector3D, r::Real, i::Int)
    arr = mv |> Vector
    arr[i] = f(arr[i], r)
    arr |> MultiVector3D
end

+(mv::MultiVector3D, v::Vector3D) = _apply(+, mv, v, 2:4)
-(mv::MultiVector3D, v::Vector3D) = _apply(-, mv, v, 2:4)

+(mv::MultiVector3D, v::PseudoVector3D) = _apply(+, mv, v, 5:7)
-(mv::MultiVector3D, v::PseudoVector3D) = _apply(-, mv, v, 5:7)

+(mv::MultiVector3D, r::Real) = _apply(+, mv, r, mv |> firstindex)
-(mv::MultiVector3D, r::Real) = _apply(-, mv, r, mv |> firstindex)

+(mv::MultiVector3D, z::Complex) = _apply(+, mv, z, mv .|> [firstindex, lastindex])
-(mv::MultiVector3D, z::Complex) = _apply(-, mv, z, mv .|> [firstindex, lastindex])

+(v::Union{Real, Complex, Vector3D, PseudoVector3D}, mv::MultiVector3D) = +(mv, v)
-(v::Union{Real, Complex, Vector3D, PseudoVector3D}, mv::MultiVector3D) = -(mv, v)

+(v::Vector3D, pv::PseudoVector3D) = [false, v...,  pv..., false] |> MultiVector3D
-(v::Vector3D, pv::PseudoVector3D) = [false, v..., -pv..., false] |> MultiVector3D

+(v::Vector3D, r::Real) = [ r, v..., zeros(Bool, 4)...] |> MultiVector3D
-(v::Vector3D, r::Real) = [-r, v..., zeros(Bool, 4)...] |> MultiVector3D

+(v::Vector3D, z::Complex) = [ real(z), v..., zeros(Bool, 3)...,  imag(z)] |> MultiVector3D
-(v::Vector3D, z::Complex) = [-real(z), v..., zeros(Bool, 3)..., -imag(z)] |> MultiVector3D

@inline +(vv::Union{PseudoVector3D, Real, Complex}, v::Vector3D) = +(v, vv)
@inline -(vv::Union{PseudoVector3D, Real, Complex}, v::Vector3D) = -(v, vv)

+(pv::PseudoVector3D, r::Real) = [ r, zeros(Bool, 3)..., pv..., false] |> MultiVector3D
-(pv::PseudoVector3D, r::Real) = [-r, zeros(Bool, 3)..., pv..., false] |> MultiVector3D

+(pv::PseudoVector3D, z::Complex) = [ real(z), zeros(Bool, 3)..., pv...,  imag(z)] |> MultiVector3D
-(pv::PseudoVector3D, z::Complex) = [-real(z), zeros(Bool, 3)..., pv..., -imag(z)] |> MultiVector3D

@inline +(v::Union{Real, Complex}, pv::PseudoVector3D) = +( pv, v)
@inline -(v::Union{Real, Complex}, pv::PseudoVector3D) = +(-pv, v)

@inline +(v1::Union{MultiVector3D, Vector3D, PseudoVector3D, Real, Complex},
          v2::Union{MultiVector3D, Vector3D, PseudoVector3D, Real, Complex},
          vs::Union{MultiVector3D, Vector3D, PseudoVector3D, Real, Complex}...) = v1 + +(v2, vs...)
@inline -(v1::Union{MultiVector3D, Vector3D, PseudoVector3D, Real, Complex},
          v2::Union{MultiVector3D, Vector3D, PseudoVector3D, Real, Complex},
          vs::Union{MultiVector3D, Vector3D, PseudoVector3D, Real, Complex}...) = v1 - +(v2, vs...)

max_stars_repo_name: Samuele-Colombo/Multivectors

max_issues_repo_name: Samuele-Colombo/Multivectors

max_forks_repo_name: Samuele-Colombo/Multivectors


text:
module PropBiology

    using LinearAlgebra
    using Interpolations

    # include("vlm_kernel.jl")
    # include("Airfoil.jl")
    include("Simulation.jl")

end # module

max_stars_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/PropBiology


text:
"""
    dndlnMh(pk, z, Δm, Ωm, lnMh; t10MF = false)

This function returns a value of the comoving halo number density per log mass interval in units of h^3/Mpc^3 at a specified value of the log(halo mass).

# Arguments
- `pk`(k): a function which returns a linear matter power spectrum with the argument k being the comoving wavenumber.
   - `pk` is in units of (Mpc/h)^3.
- `z::Real`: redshift.
- Δm: halo overdensity with respect to the mean mass density of the Universe
- `Ωm::Real`: present-day matter density parameter. When excluding the neutrino contribution, use Ωcb (baryon + cold dark matter density parameter).
- `lnMh::Real`: log(halo mass)
   - halo mass is in units of M⊙/h.

# Optional keyword arguments
- `t10MF::Bool=false`: if `true`, use `tinker10MF` for the halo multiplicity function. If `false` (the default), use `tinker08MF`.
"""
function dndlnMh(pk, z, Δm, Ωm, lnMh; t10MF = false)
   ρc = 2.775e11 # in units of h^2 M⊙/Mpc^3
   Rh = cbrt(exp(lnMh) * 3 / 4π / ρc / Ωm) # in units of Mpc/h
   σ2 = sigma2(pk, Rh)
   dlnσ2dlnRh = Rh * dsigma2dR(pk, Rh) / σ2
   lnν = 2 * log(1.6865) - log(σ2)
   if t10MF
      MF = tinker10MF(lnν, z, Δm)
   else
      MF = tinker08MF(lnν, z, Δm)
   end
   dndlnMh = -dlnσ2dlnRh * MF / 4π / Rh^3 # in units of h^3 Mpc^-3
end

"""
    bdndlnMh(pk, z, Δm, Ωm, lnMh; t10MF = false)

This function returns a value of the linear halo bias times the comoving halo number density per log mass interval in units of h^3/Mpc^3 at a specified value of the log(halo mass).

# Arguments
- `pk`(k): a function which returns a linear matter power spectrum with the argument k being the comoving wavenumber.
   - `pk` is in units of (Mpc/h)^3.
- `z::Real`: redshift.
- Δm: halo overdensity with respect to the mean mass density of the Universe
- `Ωm::Real`: present-day matter density parameter. When excluding the neutrino contribution, use Ωcb (baryon + cold dark matter density parameter).
- `lnMh::Real`: log(halo mass)
   - halo mass is in units of M⊙/h.

# Optional keyword arguments
- `t10MF::Bool=false`: if `true`, use `tinker10MF` for the halo multiplicity function. If `false` (the default), use `tinker08MF`.
"""
function bdndlnMh(pk, z, Δm, Ωm, lnMh; t10MF = false)
   ρc = 2.775e11 # in units of h^2 M⊙/Mpc^3
   Rh = cbrt(exp(lnMh) * 3 / 4π / ρc / Ωm) # in units of Mpc/h
   σ2 = sigma2(pk, Rh)
   dlnσ2dlnRh = Rh * dsigma2dR(pk, Rh) / σ2
   lnν = 2 * log(1.6865) - log(σ2)
   if t10MF
      MF = tinker10MF(lnν, z, Δm)
   else
      MF = tinker08MF(lnν, z, Δm)
   end
   # Linear halo bias, Equation (6) of Tinker et al., Astrophysical Journal, 724, 878–886 (2010)
   # Note that ν(here) is [ν(Tinker)]^2; thus, 0.5*lnν(here) = lnν(Tinker)
   y = log10(Δm)
   A = 1 + 0.24y * exp(-(4 / y)^4)
   a = 0.44y - 0.88
   C = 0.019 + 0.107y + 0.19 * exp(-(4 / y)^4)
   B, b, c = 0.183, 1.5, 2.4
   δc = 1.686
   bν =
      1 - A * exp(0.5 * a * lnν) / (exp(0.5 * a * lnν) + δc^a) +
      B * exp(0.5 * b * lnν) +
      C * exp(0.5 * c * lnν)
   bdndlnMh = -bν * dlnσ2dlnRh * MF / 4π / Rh^3 # in units of h^3 Mpc^-3
end


text:
"""
    transport!(ts, stt, acc)
1. filter `stt.dispatch_queue` for dispatch to be executed at current timestep `ts` and \
remove filtered dispatch from `stt.dispatch_queue`
2. update `stt.current_stock` given the executed dispatch at the current timestep `ts` \
and add the executed dispatch to `acc.executed_dispatch`
"""
function transport!(ts::Int, stt::States, acc::Accumulators)
    to_execute = MetaDigraph{locper}()
    for a in Iterators.filter(p -> src(p).per == ts, arcs(stt.dispatch_queue))
        add_arc!(to_execute, a)
        set_props!(to_execute, a, stt.dispatch_queue[a])
        rem_arc!(stt.dispatch_queue, a)
    end

    for a in arcs(to_execute)
        for k in keys(to_execute[a][:flow])
            stt.current_stock[
                stt.current_stock.id.==src(a).loc.&&stt.current_stock.pecahan.==k,
                :value
            ] .-= to_execute[a][:flow][k] # stock reduction at source
            stt.current_stock[
                stt.current_stock.id.==tgt(a).loc.&&stt.current_stock.pecahan.==k,
                :value
            ] .+= to_execute[a][:flow][k] # stock addition at target
        end

        add_arc!(acc.executed_dispatch, a)
        set_props!(acc.executed_dispatch, a, to_execute[a])
    end

    return nothing
end

transport!(sim::Simulation) = transport!(sim.t, sim.stt, sim.acc)


text:
# This file was generated, do not modify it. # hide
X_uni = select(X, :LStat) # only a single feature
mach_uni = machine(mdl, X_uni, y)
fit!(mach_uni)


text:
using Memoize

"""
Store functional properties of univariate functions over a domain interval
"""
struct UVFuncProp
    convex::Bool
    concave::Bool
    monotone::Bool
    lb::Float64             #lower bound on the variable range (domain) over which the properties hold
    ub::Float64             #upper bound on the variable range (domain) over which the properties hold
end

"""
Type to store a univariate function together with its properties
"""
mutable struct UVFunc
    f::Function
    prop::Array{UVFuncProp, 1}         # properties are stored in an array that represent intervals of domain over which a property holds
end

#=
"""
Compute convexity of a univariate function over a given domain

# Input:
- `f`: univariate function
- `dom`: a tuple of (lowerbound, upperbound, variable type)
"""
function compute_convexity(f::Function, dom::Tuple{Float64, Float64, Symbol})::Bool
end

"""
Compute monotonicity of a univariate function over a given domain

# Input:
- `f`: univariate function
- `dom`: a tuple of (lowerbound, upperbound, variable type)
"""
function compute_convexity(f::Function, dom::Tuple{Float64, Float64, Symbol})::Bool
end
=#

#=
"""
Identify the functional properties of the univariate functions for each variable in the constraint.

# Input:
- `func`
"""
function identify_property!(func::UVFunc)
    # each identified property over an interval is added at the field `prop` of func
end
=#

"""
Update the functional properties of the univariate function using the information given as input.

# Input:
- `func`
- `p`: Array of interval properties
"""
function update_property!(func::UVFunc, p::Array{UVFuncProp, 1})
    func.prop = copy(p)
end

"""
Compute a lower bound (minimum) of a univariate function over a domain

# Input:
- `f`: univariate function
- `prop`: array of intervals in the function domain, where each interval has functional properties
- `dom`: a tuple of (lowerbound, upperbound, variable type)
"""
function uvf_lb(f::Function, prop::Array{UVFuncProp, 1}, dom::Tuple{Float64, Float64, Symbol})::Float64
    # prop stores consequitive intervals on variable domain together with function properties over that interval
    lb = -Inf
    # find the intervals that cover the input domain
    for i=1:length(prop)
        if (dom[1] >= prop[i].lb) && (dom[1] <= prop[i].ub)     # the position of lb of domain is found,
            lb = Inf
            j = i
            while (dom[2] >= prop[j].lb)
                lb1 = max(prop[j].lb, dom[1])           # finds the startpoint of interval
                ub1 = min(prop[j].ub, dom[2])           # finds the endpoint of interval
                # if the function is monotone or concave, evaluate the lower and upper bound only
                if prop[j].monotone || prop[j].concave
                    if (dom[3] == :Int)
                        lb0 = min(f(ceil(lb1)), f(floor(ub1)))
                    else
                        lb0 = min(f(lb1), f(ub1))
                    end
                # if the variable is integer and the range is not too large, an enumeration yeilds an exact minimum over domain
                elseif (dom[3] == :Int) && (ub1 - lb1 <= 1000)
                    lb0 = Inf
                    for j in ceil(Int,lb1):floor(Int,ub1)
                        lb0 = min(f(float(j)), lb0)
                    end
                #TODO: if the function is convex a lowerbound is computed using univariate minimization packages and then using convexity inequality
                #TODO: otherwise, a univariate minimization should be run using optim.jl, howerver there is no guarantee for lower bound!
                #TODO: one way is to decompose the univariate function into separable underestimator (if it is a product), or serveral simple functions with above properties, and take lb of each
                else
                    lb0 = -Inf
                end
                lb = min(lb, lb0)
                j = j+1
                if j > length(prop)
                    break
                end
            end
            break
        end
    end
    return lb
end

"""
Compute an upper bound (maximum) of a univariate function over a domain

# Input:
- `f`: univariate function
- `prop`: array of intervals in the function domain, where each interval has functional properties
- `dom`: a tuple of (lowerbound, upperbound, variable type)
"""
function uvf_ub(f::Function, prop::Array{UVFuncProp, 1}, dom::Tuple{Float64, Float64, Symbol})::Float64
    # prop stores consequitive intervals on variable domain together with function properties over that interval
    ub = Inf
    # find the intervals that cover the input domain
    for i=1:length(prop)
        if (dom[1] >= prop[i].lb) && (dom[1] <= prop[i].ub)     # the position of lb of domain is found,
            ub = -Inf
            j = i
            while (dom[2] >= prop[j].lb)
                lb1 = max(prop[j].lb, dom[1])           # finds the startpoint of interval
                ub1 = min(prop[j].ub, dom[2])           # finds the endpoint of interval
                # if the function is monotone or concave, evaluate the lower and upper bound only
                if prop[j].monotone || prop[j].convex
                    if (dom[3] == :Int)
                        ub0 = max(f(ceil(lb1)), f(floor(ub1)))
                    else
                        ub0 = max(f(lb1), f(ub1))
                    end
                # if the variable is integer and the range is not too large, an enumeration yeilds an exact minimum over domain
            elseif (dom[3] == :Int) && (ub1 - lb1 <= 1000)
                    ub0 = -Inf
                    for j in ceil(Int,lb1):floor(Int,ub1)
                        lb0 = max(f(float(j)), ub0)
                    end
                #TODO: if the function is convex a lowerbound is computed using univariate minimization packages and then using convexity inequality
                #TODO: otherwise, a univariate minimization should be run using optim.jl, howerver there is no guarantee for lower bound!
                #TODO: one way is to decompose the univariate function into separable underestimator (if it is a product), or serveral simple functions with above properties, and take lb of each
                else
                    ub0 = Inf
                end
                ub = max(ub, ub0)
                j = j+1
                if j > length(prop)
                    break
                end
            end
            break
        end
    end
    return ub
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
@testset "Cost Matrix/String" begin
    @test distance_matrix("ab", "aab", deletion = 2) == 1                # Should go with INSERTION instead of DELETION.
    @test distance_matrix("ab", "aab", deletion = 2, insertion = 2) == 2 # Could go with either.
    @test distance_matrix("abc", "bbc", substitution = 2) == 2           # DELETION + INSERTION = SUBSTITUTION
end

@testset "Cost Matrix/Array" begin
    @test distance_matrix([ 1 2 ], [ 1 1 2 ], deletion = 2) == 1                # Should go with INSERTION instead of DELETION.
    @test distance_matrix([ 1 2 ], [ 1 1 2 ], deletion = 2, insertion = 2) == 2 # Could go with either.
    @test distance_matrix([ 1 2 3 ], [ 2 2 3 ], substitution = 2) == 2          # DELETION + INSERTION = SUBSTITUTION
end

max_stars_repo_name: jlegare/LevenshteinToolkit

max_issues_repo_name: jlegare/LevenshteinToolkit

max_forks_repo_name: jlegare/LevenshteinToolkit


text:
abstract type Direction end
struct Forward <: Direction end
struct Up <: Direction end
struct Down <: Direction end
const Command = Tuple{<:Direction,<:Integer}

const directions = Dict("forward" => Forward(), "up" => Up(), "down" => Down())

# Load data
function load(::Day{2}, path)
    return map(readlines(path)) do line
        dir, incr = split(line, " ")
        directions[dir], parse(Int, incr)
    end
end

mutable struct Submarine{T<:Integer}
    depth::T
    hpos::T
    aim::T
end
Submarine() = Submarine(0, 0, 0)

# Part 1
move_submarine!(s::Submarine, c::Command) = move_submarine!(s, c...)
move_submarine!(s::Submarine, ::Forward, incr) = (s.hpos += incr)
move_submarine!(s::Submarine, ::Down, incr) = (s.depth += incr)
move_submarine!(s::Submarine, ::Up, incr) = (s.depth -= incr)

# Part 2
move_submarine2!(s::Submarine, c::Command) = move_submarine2!(s, c...)
move_submarine2!(s::Submarine, ::Down, incr) = (s.aim += incr)
move_submarine2!(s::Submarine, ::Up, incr) = (s.aim -= incr)
function move_submarine2!(s::Submarine, ::Forward, incr)
    s.hpos += incr
    s.depth += s.aim * incr
    return nothing
end

function solve1(::Day{2}, data)
    sub = Submarine()
    for command in data
        move_submarine!(sub, command)
    end
    return sub.depth * sub.hpos
end

function solve2(::Day{2}, data)
    sub = Submarine()
    for command in data
        move_submarine2!(sub, command)
    end
    return sub.depth * sub.hpos
end

testresult1(::Day{2}) = 150
testresult2(::Day{2}) = 900


text:
include("util.jl")

import Base.Broadcast: BroadcastFunction

# Function definition for every possible broadcast op involving a Tensor
# Needed to keep track of computation for backward pass
for (F, S) in ((Tensor, Tensor), (Tensor, Any), (Any, Tensor))
    @eval begin
        # Generic broadcast function for computations involving a Tensor
        function Broadcast.broadcasted(op::Any, x::$F, y::$S)
            Tensor(broadcast(op, data(x), data(y)), op=BroadcastFunction{typeof(op)}, prev=(x, y))
        end
    end
end

max_stars_repo_name: leonvol/magrad

max_issues_repo_name: leonvol/magrad

max_forks_repo_name: leonvol/magrad


text:
################################################################################
# fmpz_mat -> nemo file
# use as include(...)
################################################################################
function toNemo(io::IOStream, A::fmpz_mat; name = "A")
  println(io, name, " = MatrixSpace(FlintZZ, ", nrows(A), ", ", ncols(A), ")([")
  for i = 1:nrows(A)
    for j = 1:ncols(A)
      print(io, A[i,j])
      if j < ncols(A)
        print(io, " ")
      end
    end
    if i<nrows(A)
      println(io, ";")
    end
  end
  println(io, "]);")
  println(io, "println(\"Loaded ", name, "\")")
end

function toNemo(s::String, A::fmpz_mat)
  f = open(s, "w")
  toNemo(f, A)
  close(f)
end

################################################################################
# fmpz_mat -> magma file
# use as read(...)
################################################################################
function toMagma(io::IOStream, A::fmpz_mat; name = "A")
  println(io, name, " := Matrix(Integers(), ", nrows(A), ", ", ncols(A), ", [")
  for i = 1:nrows(A)
    for j = 1:ncols(A)
      print(io, A[i,j])
      if j < ncols(A)
        print(io, ", ")
      end
    end
    if i<nrows(A)
      println(io, ",")
    end
  end
  println(io, "]);")
  println(io, "\"Loaded ", name, "\";")
end

function toMagma(s::String, A::fmpz_mat)
  f = open(s, "w")
  toMagma(f, A)
  close(f)
end

################################################################################
# SMat -> magma file
# use as read(...)
################################################################################
function toMagma(io::IOStream, A::SMat; name = "A")
  println(io, name, " := SparseMatrix(Integers(), ", nrows(A), ", ", ncols(A), ", [")
  for i = 1:nrows(A)
    for xx = 1:length(A.rows[i].pos)
      print(io, "<", i, ", ", A.rows[i].pos[xx], ", ", A.rows[i].values[xx], ">")
      if xx < length(A.rows[i].pos) || i<nrows(A)
        print(io, ", ")
      end
    end
    println(io, "")
  end
  println(io, "]);")
  println(io, "\"Loaded ", name, "\";")
end

function toMagma(s::String, A::SMat)
  f = open(s, "w")
  toMagma(f, A)
  close(f)
end


text:
blockrowsupport(_, A, k) = blockaxes(A,2)
""""
    blockrowsupport(A, k)

gives an iterator containing the possible non-zero blocks in the k-th block-row of A.
"""
blockrowsupport(A, k) = blockrowsupport(MemoryLayout(A), A, k)
blockrowsupport(A) = blockrowsupport(A, blockaxes(A,1))

blockcolsupport(_, A, j) = Block.(colsupport(blocks(A), Int.(j)))

""""
    blockcolsupport(A, j)

gives an iterator containing the possible non-zero blocks in the j-th block-column of A.
"""
blockcolsupport(A, j) = blockcolsupport(MemoryLayout(A), A, j)
blockcolsupport(A) = blockcolsupport(A, blockaxes(A,2))

blockcolstart(A...) = first(blockcolsupport(A...))
blockcolstop(A...) = last(blockcolsupport(A...))
blockrowstart(A...) = first(blockrowsupport(A...))
blockrowstop(A...) = last(blockrowsupport(A...))

for Func in (:blockcolstart, :blockcolstop, :blockrowstart, :blockrowstop)
    @eval $Func(A, i::Block{1}) = $Func(A, Int(i))
end

abstract type AbstractBlockLayout <: MemoryLayout end
struct BlockLayout{ArrLay,BlockLay} <: AbstractBlockLayout end

function colsupport(::BlockLayout, A, j)
    JR = Int(findblock(axes(A,2), minimum(j))):Int(findblock(axes(A,2), maximum(j)))
    KR = colsupport(blocks(A), JR)
    axes(A,1)[Block.(KR)]
end

function rowsupport(::BlockLayout, A, k)
    KR = Int(findblock(axes(A,1), minimum(k))):Int(findblock(axes(A,1), maximum(k)))
    JR = rowsupport(blocks(A), KR)
    axes(A,2)[Block.(JR)]
end

similar(M::MulAdd{<:AbstractBlockLayout,<:AbstractBlockLayout}, ::Type{T}, axes) where {T,N} =
    similar(BlockArray{T}, axes)

@inline MemoryLayout(::Type{<:PseudoBlockArray{T,N,R}}) where {T,N,R} = MemoryLayout(R)
@inline MemoryLayout(::Type{<:BlockArray{T,N,R}}) where {T,N,D,R<:AbstractArray{D,N}} =
    BlockLayout{typeof(MemoryLayout(R)),typeof(MemoryLayout(D))}()

sublayout(::BlockLayout{MLAY,BLAY}, ::Type{<:NTuple{N,BlockSlice1}}) where {MLAY,BLAY,N} = BLAY()
sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:NTuple{N,<:BlockSlice{BlockRange{1,Tuple{II}}}}}) where {N,MLAY,BLAY,II} =
    BlockLayout{typeof(sublayout(MLAY(),NTuple{N,II})), BLAY}()
sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:Tuple{BlockSlice1,<:BlockSlice{BlockRange{1,Tuple{II}}}}}) where {MLAY,BLAY,II} =
    BlockLayout{typeof(sublayout(MLAY(),Tuple{Int,II})), BLAY}()
sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:Tuple{BlockSlice{BlockRange{1,Tuple{II}}},BlockSlice1}}) where {MLAY,BLAY,II} =
    BlockLayout{typeof(sublayout(MLAY(),Tuple{II,Int})), BLAY}()
# This might need modification: no guarantee axes(BL,1) == axes(MLAY,1) so Slice might not be right here
sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:Tuple{Sl,<:BlockSlice{BlockRange{1,Tuple{II}}}}}) where {MLAY,BLAY,Sl<:Slice,II} =
    BlockLayout{typeof(sublayout(MLAY(),Tuple{Sl,II})), BLAY}()
sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:Tuple{<:BlockSlice{BlockRange{1,Tuple{II}}},Sl}}) where {MLAY,BLAY,Sl<:Slice,II} =
    BlockLayout{typeof(sublayout(MLAY(),Tuple{II,Sl})), BLAY}()
sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:Tuple{Sl1,Sl2}}) where {MLAY,BLAY,Sl1<:Slice,Sl2<:Slice,II} =
    BlockLayout{typeof(sublayout(MLAY(),Tuple{Sl1,Sl2})), BLAY}()

# materialize views, used for `getindex`
sub_materialize(::AbstractBlockLayout, V, _) = BlockArray(V)
sub_materialize(::AbstractBlockLayout, V, ::Tuple{<:BlockedUnitRange}) = BlockArray(V)
sub_materialize(::AbstractBlockLayout, V, ::Tuple{<:BlockedUnitRange,<:BlockedUnitRange}) = BlockArray(V)
sub_materialize(::AbstractBlockLayout, V, ::Tuple{<:AbstractUnitRange,<:BlockedUnitRange}) = BlockArray(V)
sub_materialize(::AbstractBlockLayout, V, ::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) = BlockArray(V)

# if it's not a block layout, best to use PseudoBlockArray to take advantage of strideness
sub_materialize(_, V, ::Tuple{<:BlockedUnitRange}) = PseudoBlockArray(V)
sub_materialize(_, V, ::Tuple{<:BlockedUnitRange,<:BlockedUnitRange}) = PseudoBlockArray(V)
sub_materialize(_, V, ::Tuple{<:AbstractUnitRange,<:BlockedUnitRange}) = PseudoBlockArray(V)
sub_materialize(_, V, ::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) = PseudoBlockArray(V)

# Special for FillArrays.jl

# special case for fill blocks
LinearAlgebra.fill!(V::SubArray{T,1,<:BlockArray,<:Tuple{BlockSlice1}}, x) where T =
    fill!(view(parent(V), parentindices(V)[1].block), x)

FillArrays.getindex_value(V::SubArray{T,1,<:BlockArray,<:Tuple{BlockSlice1}}) where T =
    FillArrays.getindex_value(view(parent(V), block(parentindices(V)[1])))

sub_materialize(::ArrayLayouts.AbstractFillLayout, V, ax::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) =
    Fill(FillArrays.getindex_value(V), ax)
sub_materialize(::ArrayLayouts.OnesLayout, V, ax::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) =
    Ones{eltype(V)}(ax)
sub_materialize(::ArrayLayouts.ZerosLayout, V, ax::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) =
    Zeros{eltype(V)}(ax)
sub_materialize(::ArrayLayouts.AbstractFillLayout, V, ax::Tuple{<:AbstractUnitRange,<:BlockedUnitRange}) =
    Fill(FillArrays.getindex_value(V), ax)
sub_materialize(::ArrayLayouts.OnesLayout, V, ax::Tuple{<:AbstractUnitRange,<:BlockedUnitRange}) =
    Ones{eltype(V)}(ax)
sub_materialize(::ArrayLayouts.ZerosLayout, V, ax::Tuple{<:AbstractUnitRange,<:BlockedUnitRange}) =
    Zeros{eltype(V)}(ax)

conjlayout(::Type{T}, ::BlockLayout{MLAY,BLAY}) where {T<:Complex,MLAY,BLAY} = BlockLayout{MLAY,typeof(conjlayout(T,BLAY()))}()
conjlayout(::Type{T}, ::BlockLayout{MLAY,BLAY}) where {T<:Real,MLAY,BLAY} = BlockLayout{MLAY,BLAY}()

transposelayout(::BlockLayout{MLAY,BLAY}) where {MLAY,BLAY} =
    BlockLayout{typeof(transposelayout(MLAY())),typeof(transposelayout(BLAY()))}()

#############
# copyto!
#############

function _copyto!(_, ::AbstractBlockLayout, dest::AbstractVector, src::AbstractVector)
    if !blockisequal(axes(dest), axes(src))
        # impose block structure
        copyto!(PseudoBlockArray(dest, axes(src)), src)
        return dest
    end

    @inbounds for K = blockaxes(src,1)
        copyto!(view(dest,K), view(src,K))
    end
    dest
end

function _copyto!(_, ::AbstractBlockLayout, dest::AbstractMatrix, src::AbstractMatrix)
    if !blockisequal(axes(dest), axes(src))
        # impose block structure
        copyto!(PseudoBlockArray(dest, axes(src)), src)
        return dest
    end

    @inbounds for J = blockaxes(src,2),K = blockaxes(src,1)
        copyto!(view(dest,K,J), view(src,K,J))
    end
    dest
end

#############
# BLAS overrides
#############

function materialize!(M::MatMulVecAdd{<:AbstractBlockLayout,<:AbstractStridedLayout,<:AbstractStridedLayout})
    α, A, x_in, β, y_in = M.α, M.A, M.B, M.β, M.C
    if length(x_in) != size(A,2) || length(y_in) != size(A,1)
        throw(DimensionMismatch())
    end

    # impose block structure
    y = PseudoBlockArray(y_in, (axes(A,1),))
    x = PseudoBlockArray(x_in, (axes(A,2),))

    _fill_lmul!(β, y)

    for J = blockaxes(A,2)
        for K = blockcolsupport(A,J)
            muladd!(α, view(A,K,J), view(x,J), one(α), view(y,K))
        end
    end
    y_in
end

function _block_muladd!(α, A, X, β, Y)
    _fill_lmul!(β, Y)
    for J = blockaxes(X,2), N = blockcolsupport(X,J), K = blockcolsupport(A,N)
        muladd!(α, view(A,K,N), view(X,N,J), one(α), view(Y,K,J))
    end
    Y
end

mul_blockscompatible(A, B, C) = blockisequal(axes(A,2), axes(B,1)) &&
    blockisequal(axes(A,1), axes(C,1)) &&
    blockisequal(axes(B,2), axes(C,2))

function materialize!(M::MatMulMatAdd{<:AbstractBlockLayout,<:AbstractBlockLayout,<:AbstractBlockLayout})
    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
    if mul_blockscompatible(A,B,C)
        _block_muladd!(α, A, B, β, C)
    else # use default
        materialize!(MulAdd{UnknownLayout,UnknownLayout,UnknownLayout}(α, A, B, β, C))
    end
end

function materialize!(M::MatMulMatAdd{<:AbstractBlockLayout,<:AbstractBlockLayout,<:AbstractColumnMajor})
    α, A, X, β, Y_in = M.α, M.A, M.B, M.β, M.C
    Y = PseudoBlockArray(Y_in, (axes(A,1), axes(Y_in,2)))
    _block_muladd!(α, A, X, β, Y)
    Y_in
end

function materialize!(M::MatMulMatAdd{<:AbstractBlockLayout,<:AbstractColumnMajor,<:AbstractColumnMajor})
    α, A, X_in, β, Y_in = M.α, M.A, M.B, M.β, M.C
    _fill_lmul!(β, Y_in)
    X = PseudoBlockArray(X_in, (axes(A,2), axes(X_in,2)))
    Y = PseudoBlockArray(Y_in, (axes(A,1), axes(Y_in,2)))
    _block_muladd!(α, A, X, β, Y)
    Y_in
end

function materialize!(M::MatMulMatAdd{<:AbstractColumnMajor,<:AbstractBlockLayout,<:AbstractColumnMajor})
    α, A_in, X, β, Y_in = M.α, M.A, M.B, M.β, M.C
    _fill_lmul!(β, Y_in)
    A = PseudoBlockArray(A_in, (axes(A_in,1),axes(X,1)))
    Y = PseudoBlockArray(Y_in, (axes(Y_in,1),axes(X,2)))
    _block_muladd!(α, A, X, β, Y)
    Y_in
end

####
# Triangular
####

@inline hasmatchingblocks(A) = blockisequal(axes(A)...)

triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:AbstractBlockLayout} = Tri{ML}()

_triangular_matrix(::Val{'U'}, ::Val{'N'}, A) = UpperTriangular(A)
_triangular_matrix(::Val{'L'}, ::Val{'N'}, A) = LowerTriangular(A)
_triangular_matrix(::Val{'U'}, ::Val{'U'}, A) = UnitUpperTriangular(A)
_triangular_matrix(::Val{'L'}, ::Val{'U'}, A) = UnitLowerTriangular(A)

function _matchingblocks_triangular_mul!(::Val{'U'}, UNIT, A::AbstractMatrix{T}, dest) where T
    # impose block structure
    b = PseudoBlockArray(dest, (axes(A,1),))

    for K = blockaxes(A,1)
        b_2 = view(b, K)
        Ũ = _triangular_matrix(Val('U'), UNIT, view(A, K,K))
        materialize!(Lmul(Ũ, b_2))
        JR = (K+1):last(blockrowsupport(A,K))
        if !isempty(JR)
            muladd!(one(T), view(A, K, JR), view(b,JR), one(T), b_2)
        end
    end
    dest
end

function _matchingblocks_triangular_mul!(::Val{'L'}, UNIT, A::AbstractMatrix{T}, dest) where T
    # impose block structure
    b = PseudoBlockArray(dest, (axes(A,1),))

    N = blocksize(A,1)

    for K = N:-1:1
        b_2 = view(b, Block(K))
        L̃ = _triangular_matrix(Val('L'), UNIT, view(A, Block(K,K)))
        materialize!(Lmul(L̃, b_2))
        JR = blockrowstart(A,K):Block(K-1)
        if !isempty(JR)
            muladd!(one(T), view(A, Block(K), JR), view(b,JR), one(T), b_2)
        end
    end

    dest
end

@inline function materialize!(M::MatLmulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractBlockLayout},
                                   <:AbstractStridedLayout}) where {UPLO,UNIT}
    U,x = M.A,M.B
    T = eltype(M)
    @boundscheck size(U,1) == size(x,1) || throw(BoundsError())
    if hasmatchingblocks(U)
        _matchingblocks_triangular_mul!(Val(UPLO), Val(UNIT), triangulardata(U), x)
    else # use default
        x_1 = PseudoBlockArray(copy(x), (axes(U,2),))
        x_2 = PseudoBlockArray(x, (axes(U,1),))
        _block_muladd!(one(T), U, x_1, zero(T), x_2)
    end
end

for UNIT in ('U', 'N')
    @eval begin
        @inline function materialize!(M::MatLdivVec{<:TriangularLayout{'U',$UNIT,<:AbstractBlockLayout},
                                        <:AbstractStridedLayout})
            U,dest = M.A,M.B
            T = eltype(dest)

            A = triangulardata(U)
            if !hasmatchingblocks(A) # Use default for now
                return materialize!(Ldiv{TriangularLayout{'U',$UNIT,UnknownLayout},
                                         typeof(MemoryLayout(dest))}(U, dest))
            end

            @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())

            # impose block structure
            b = PseudoBlockArray(dest, (axes(A,1),))

            N = blocksize(A,1)

            for K = N:-1:1
                b_2 = view(b, Block(K))
                Ũ = _triangular_matrix(Val('U'), Val($UNIT), view(A, Block(K,K)))
                materialize!(Ldiv(Ũ, b_2))

                if K ≥ 2
                    KR = blockcolstart(A, K):Block(K-1)
                    V_12 = view(A, KR, Block(K))
                    b̃_1 = view(b, KR)
                    muladd!(-one(T), V_12, b_2, one(T), b̃_1)
                end
            end

            dest
        end

        @inline function materialize!(M::MatLdivVec{<:TriangularLayout{'L',$UNIT,<:AbstractBlockLayout},
                                        <:AbstractStridedLayout})
            L,dest = M.A, M.B
            T = eltype(dest)
            A = triangulardata(L)
            if !hasmatchingblocks(A) # Use default for now
                return materialize!(Ldiv{TriangularLayout{'L',$UNIT,UnknownLayout},
                                         typeof(MemoryLayout(dest))}(L, dest))
            end

            @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())

            # impose block structure
            b = PseudoBlockArray(dest, (axes(A,1),))

            N = blocksize(A,1)

            for K = 1:N
                b_2 = view(b, Block(K))
                L̃ = _triangular_matrix(Val('L'), Val($UNIT), view(A, Block(K,K)))
                materialize!(Ldiv(L̃, b_2))

                if K < N
                    KR = Block(K+1):blockcolstop(A, K)
                    V_12 = view(A, KR, Block(K))
                    b̃_1 = view(b, KR)
                    muladd!(-one(T), V_12, b_2, one(T), b̃_1)
                end
            end

            dest
        end
    end
end

# For now, use PseudoBlockArray
_inv(::AbstractBlockLayout, axes, A) = BlockArray(inv(PseudoBlockArray(A)))

for op in (:exp, :log, :sqrt)
    @eval begin
       $op(A::PseudoBlockMatrix) = PseudoBlockMatrix($op(A.blocks), axes(A))
    end
end


text:
# State Constructors
export pure_state, mixed_state
export bloch_qubit_state
export computational_basis_states
export bell_states, generalized_bell_states

# State ensemble constructors
export planar_symmetric_qubit_states, mirror_symmetric_qubit_states
export trine_qubit_states, bb84_qubit_states, sic_qubit_states

"""
    pure_state( ψ :: Ket; atol=ATOL :: Float64 ) :: State
    pure_state( ψ :: Bra; atol=ATOL :: Float64 ) :: State

A [`State`](@ref) is "pure" if it is rank-one (see [`is_pure`](@ref)).
A rank-one density matrix is constructed by taking the outer-product of a
[`Ket`](@ref) (or [`Bra`](@ref)) with itself, ``|\\psi\\rangle\\langle \\psi| = \\rho``.
This method alternatively accepts `Vector` inputs.

```julia
pure_state( ψ :: Vector ) :: State
pure_state( ψ :: Adjoint{T,Vector{T}} where T; atol=ATOL :: Float64 ) :: State
```
"""
function pure_state(ket::Ket; atol=ATOL :: Float64) :: State
    State(ket.ψ*ket.ψ', atol=atol)
end
function pure_state(bra::Bra; atol=ATOL :: Float64) :: State
    State(bra.ψ'*bra.ψ, atol=atol)
end
function pure_state(ψ::Vector{<:Number}; atol=ATOL :: Float64) :: State
    State(ψ*ψ', atol=atol)
end
function pure_state(
    ψ::Adjoint{T,Vector{T}} where T; atol=ATOL :: Float64
) :: State
    State(ψ'*ψ, atol=atol)
end

"""
    mixed_state(
        priors :: Probabilities,
        states :: Vector{<:State};
        atol=ATOL :: Float64
    ) :: State

Constructs the statistical mixture (weighted average) of quantum states.
This method can also be used `Vector` and `Matrix` types which satisfy the requirements
of [`is_probability_distribution`](@ref) and [`is_density_matrix`](@ref) respectively.

```julia
mixed_state(
    priors :: AbstractVector,
    states :: Vector{<:AbstractMatrix}
) :: State
```

A `DomainError` is thrown if the provided `Vector`s of `priors` or `states` do
not satisfy their respective requirements.
"""
function mixed_state(
        priors :: Probabilities,
        states :: Vector{<:State};
        atol=ATOL :: Float64
) :: State
    State(sum(priors .* states); atol=atol)
end
function mixed_state(
    priors::AbstractVector,
    states::Vector{<:AbstractMatrix};
    atol=ATOL :: Float64
) :: State
    if !is_probability_distribution(priors)
        throw(DomainError(priors, "Input `priors` is not a valid probability distribution."))
    elseif findfirst(!is_density_matrix, states) !== nothing
        throw(DomainError(states, "Input `states` is not a set of valid density matrices."))
    end
    State(sum(priors .* states); atol=atol)
end

"""
    computational_basis_states( dim :: Int64 ) :: Vector{State}

The set of density matrices constructed from the [`computational_basis_vectors`](@ref)
with dimension, `dim`.
"""
function computational_basis_states(dim::Int64)::Vector{State}
    pure_state.(computational_basis_vectors(dim))
end

"""
    bell_states() :: Vector{State{Float64}}

The Bell basis density matrices. See [`bell_kets`](@ref) for more details.
"""
bell_states() :: Vector{State{Float64}} = pure_state.(bell_kets())

"""
    generalized_bell_states( dim :: Int64 ) :: Vector{State{ComplexF64}

The generalized Bell basis density matrices. See  [`generalized_bell_kets`](@ref) for more details.
A `DomainError` is thrown if `dim ≥ 2` is not satisfied.
"""
generalized_bell_states(dim :: Int64) :: Vector{State{ComplexF64}} = pure_state.(generalized_bell_kets(dim))

"""
    planar_symmetric_qubit_states( n :: Int64 ) :: Vector{State{Float64}}

Constructs a set of pure `State`s oriented symmetrically in the x-z-plane.
See [`planar_symmetric_qubit_kets`](@ref) for details.
"""
function planar_symmetric_qubit_states(n :: Int64) :: Vector{State{Float64}}
    pure_state.(planar_symmetric_qubit_kets(n))
end

"""
    mirror_symmetric_qubit_states( θ ::  Real ) :: Vector{State{Float64}}

Returns a set of 3 mirror symmetric qubit density matrices. The first
state is ``|0\\rangle\\langle 0|`` the other two are symmetric about the  ``|0\\rangle`` axis.
See [`mirror_symmetric_qubit_kets`](@ref).

Input:
* `θ ∈ [0,π/2]`: the hilbert space angle between ``|0\\rangle`` and ``|\\psi_{2/3}\\rangle``.
"""
function mirror_symmetric_qubit_states(θ :: Real)::Vector{State{Float64}}
    pure_state.(mirror_symmetric_qubit_kets(θ))
end

"""
Returns the qubit density matrix for quantum state described by a coordinate on
bloch sphere. See [`bloch_qubit_ket`](@ref)

Spherical Coordinates:

States on the surface of bloch sphere may be described by spherical coordinates.

```julia
bloch_qubit_state(θ::Real, ϕ::Real) :: State{Complex{Float64}}
```

* `θ ∈ [0, π]`: polar angle (w.r.t z-axis).
* `ϕ ∈ [0, 2π]`: azimuthal angle (x-y plane)

Cartesian Coordinates:

States within the volume of bloch sphere may be described in cartesian coordinates.

```julia
bloch_qubit_state(x::Real, y::Real, z::Real) :: State{Complex{Float64}}
```

 * where `x`, `y`, and `z` are constrained to the unit sphere, `0 <= norm([x,y,z]) <= 1`.

 A `DomainError` is thrown if the coordinates `(x,y,z)` do  not adhere to constraints.
"""
function bloch_qubit_state(θ::Real,ϕ::Real) :: State{Complex{Float64}}
    pure_state(bloch_qubit_ket(θ,ϕ))
end
function bloch_qubit_state(θ::Real) :: State{Float64}
    pure_state(bloch_qubit_ket(θ))
end
function bloch_qubit_state(x :: Real, y :: Real, z :: Real) :: State{Complex{Float64}}
    r = [x,y,z]

    if !(0 <= norm(r) <= 1)
        throw(DomainError(r, "norm([x,y,z]) is not in valid range"))
    end

    paulis_xyz = [[0 1;1 0],[0 -im;im 0],[1 0;0 -1]]

    ρ = 0.5*([1 0;0 1] + sum(r .* paulis_xyz))

    State(ρ)
end

"""
    trine_qubit_states() :: Vector{State{Float64}}

Returns the qubit trine states in density matrix form. See [`trine_qubit_kets`](@ref).
"""
trine_qubit_states() :: Vector{State{Float64}} = pure_state.(trine_qubit_kets())

"""
    sic_qubit_states :: Vector{State{Complex{Float64}}}

The quadruplet of symmetric informationally complete (SIC) qubits. The qubits
are the vertices of a tetrahedron inscribed on bloch sphere. See [`sic_qubit_kets`](@ref).
"""
sic_qubit_states() :: Vector{State{Complex{Float64}}} = pure_state.(sic_qubit_kets())

"""
    bb84_qubit_states :: Vector{State}

The quadruplet of qubits used in the BB84 Quantum Key Distribution protocol. The
states are ``|0\\rangle\\langle 0|``, ``|1\\rangle\\langle 1|``, ``|+\\rangle\\langle +|``, and ``|- \\rangle\\langle -|``.
See [`bb84_qubit_kets`](@ref).
"""
bb84_qubit_states() :: Vector{State{Float64}} = pure_state.(bb84_qubit_kets())


text:
using LinearAlgebra, Zygote, Plots

const n = 100
S = circshift(Matrix{Float64}(I, n, n),(1,0))

function loss(W)
    norm(W*S-S*W)/norm(W)
end

function step!(W;lr=0.003)
    # computing current loss and backprop
    current_loss, back_loss = pullback(w -> loss(w),W)
    # computing gradient
    grads = back_loss(1)[1]
    # updating W
    W .-= lr .*grads
end

W = randn(n,n)
W ./= norm(W)

heatmap(W,clims=(-0.03,0.03),legend=:none,axis=nothing)

for i=1:1000
    step!(W)
end

heatmap(W,clims=(-0.03,0.03),legend=:none,axis=nothing)