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

abstract type ClusterManager end

mutable struct WorkerConfig
    # Common fields relevant to all cluster managers
    io::Nullable{IO}
    host::Nullable{AbstractString}
    port::Nullable{Integer}

    # Used when launching additional workers at a host
    count::Nullable{Union{Int, Symbol}}
    exename::Nullable{AbstractString}
    exeflags::Nullable{Cmd}

    # External cluster managers can use this to store information at a per-worker level
    # Can be a dict if multiple fields need to be stored.
    userdata::Nullable{Any}

    # SSHManager / SSH tunnel connections to workers
    tunnel::Nullable{Bool}
    bind_addr::Nullable{AbstractString}
    sshflags::Nullable{Cmd}
    max_parallel::Nullable{Integer}

    # Used by Local/SSH managers
    connect_at::Nullable{Any}

    process::Nullable{Process}
    ospid::Nullable{Integer}

    # Private dictionary used to store temporary information by Local/SSH managers.
    environ::Nullable{Dict}

    # Connections to be setup depending on the network topology requested
    ident::Nullable{Any}      # Worker as identified by the Cluster Manager.
    # List of other worker idents this worker must connect with. Used with topology T_CUSTOM.
    connect_idents::Nullable{Array}

    # Run multithreaded blas on worker
    enable_threaded_blas::Nullable{Bool}

    function WorkerConfig()
        wc = new()
        for n in 1:length(WorkerConfig.types)
            T = eltype(fieldtype(WorkerConfig, n))
            setfield!(wc, n, Nullable{T}())
        end
        wc
    end
end

@enum WorkerState W_CREATED W_CONNECTED W_TERMINATING W_TERMINATED
mutable struct Worker
    id::Int
    del_msgs::Array{Any,1}
    add_msgs::Array{Any,1}
    gcflag::Bool
    state::WorkerState
    c_state::Condition      # wait for state changes
    ct_time::Float64        # creation time

    r_stream::IO
    w_stream::IO
    w_serializer::ClusterSerializer  # writes can happen from any task hence store the
                                     # serializer as part of the Worker object
    manager::ClusterManager
    config::WorkerConfig
    version::Nullable{VersionNumber}  # Julia version of the remote process

    function Worker(id::Int, r_stream::IO, w_stream::IO, manager::ClusterManager;
                             version=Nullable{VersionNumber}(), config=WorkerConfig())
        w = Worker(id)
        w.r_stream = r_stream
        w.w_stream = buffer_writes(w_stream)
        w.w_serializer = ClusterSerializer(w.w_stream)
        w.manager = manager
        w.config = config
        w.version = version
        set_worker_state(w, W_CONNECTED)
        register_worker_streams(w)
        w
    end

    function Worker(id::Int)
        @assert id > 0
        if haskey(map_pid_wrkr, id)
            return map_pid_wrkr[id]
        end
        w=new(id, [], [], false, W_CREATED, Condition(), time())
        register_worker(w)
        w
    end

    Worker() = Worker(get_next_pid())
end

function set_worker_state(w, state)
    w.state = state
    notify(w.c_state; all=true)
end

function check_worker_state(w::Worker)
    if w.state == W_CREATED
        if PGRP.topology == :all_to_all
            # Since higher pids connect with lower pids, the remote worker
            # may not have connected to us yet. Wait for some time.
            timeout =  worker_timeout() - (time() - w.ct_time)
            timeout <= 0 && error("peer $(w.id) has not connected to $(myid())")

            @schedule (sleep(timeout); notify(w.c_state; all=true))
            wait(w.c_state)
            w.state == W_CREATED && error("peer $(w.id) didn't connect to $(myid()) within $timeout seconds")
        else
            error("peer $(w.id) is not connected to $(myid()). Topology : " * string(PGRP.topology))
        end
    end
end

## process group creation ##

mutable struct LocalProcess
    id::Int
    bind_addr::AbstractString
    bind_port::UInt16
    cookie::AbstractString
    LocalProcess() = new(1)
end

function disable_threaded_libs()
    BLAS.set_num_threads(1)
end

worker_timeout() = parse(Float64, get(ENV, "JULIA_WORKER_TIMEOUT", "60.0"))

## worker creation and setup ##

# The entry point for julia worker processes. does not return. Used for TCP transport.
# Cluster managers implementing their own transport will provide their own.
# Argument is descriptor to write listening port # to.
start_worker(cookie::AbstractString) = start_worker(STDOUT, cookie)
function start_worker(out::IO, cookie::AbstractString)
    # we only explicitly monitor worker STDOUT on the console, so redirect
    # stderr to stdout so we can see the output.
    # at some point we might want some or all worker output to go to log
    # files instead.
    # Currently disabled since this caused processes to spin instead of
    # exit when process 1 shut down. Don't yet know why.
    #redirect_stderr(STDOUT)

    init_worker(cookie)
    interface = IPv4(LPROC.bind_addr)
    if LPROC.bind_port == 0
        (actual_port,sock) = listenany(interface, UInt16(9009))
        LPROC.bind_port = actual_port
    else
        sock = listen(interface, LPROC.bind_port)
    end
    @schedule while isopen(sock)
        client = accept(sock)
        process_messages(client, client, true)
    end
    print(out, "julia_worker:")  # print header
    print(out, "$(dec(LPROC.bind_port))#") # print port
    print(out, LPROC.bind_addr)
    print(out, '\n')
    flush(out)
    # close STDIN; workers will not use it
    #close(STDIN)

    disable_nagle(sock)

    if ccall(:jl_running_on_valgrind,Cint,()) != 0
        println(out, "PID = $(getpid())")
    end

    try
        # To prevent hanging processes on remote machines, newly launched workers exit if the
        # master process does not connect in time.
        # TODO : Make timeout configurable.
        check_master_connect()
        while true; wait(); end
    catch err
        print(STDERR, "unhandled exception on $(myid()): $(err)\nexiting.\n")
    end

    close(sock)
    exit(0)
end

function redirect_worker_output(ident, stream)
    @schedule while !eof(stream)
        line = readline(stream)
        if startswith(line, "\tFrom worker ")
            # STDOUT's of "additional" workers started from an initial worker on a host are not available
            # on the master directly - they are routed via the initial worker's STDOUT.
            println(line)
        else
            println("\tFrom worker $(ident):\t$line")
        end
    end
end

# The default TCP transport relies on the worker listening on a free
# port available and printing its bind address and port.
# The master process uses this to connect to the worker and subsequently
# setup a all-to-all network.
function read_worker_host_port(io::IO)
    while true
        conninfo = readline(io)
        bind_addr, port = parse_connection_info(conninfo)
        if bind_addr != ""
            return bind_addr, port
        end
    end
end

function parse_connection_info(str)
    m = match(r"^julia_worker:(\d+)#(.*)", str)
    if m !== nothing
        (m.captures[2], parse(Int16, m.captures[1]))
    else
        ("", Int16(-1))
    end
end

"""
    init_worker(cookie::AbstractString, manager::ClusterManager=DefaultClusterManager())

Called by cluster managers implementing custom transports. It initializes a newly launched
process as a worker. Command line argument `--worker` has the effect of initializing a
process as a worker using TCP/IP sockets for transport.
`cookie` is a [`cluster_cookie`](@ref).
"""
function init_worker(cookie::AbstractString, manager::ClusterManager=DefaultClusterManager())
    # On workers, the default cluster manager connects via TCP sockets. Custom
    # transports will need to call this function with their own manager.
    global cluster_manager
    cluster_manager = manager

    # Since our pid has yet to be set, ensure no RemoteChannel / Future  have been created or addprocs() called.
    assert(nprocs() <= 1)
    assert(isempty(PGRP.refs))
    assert(isempty(client_refs))

    # System is started in head node mode, cleanup related entries
    empty!(PGRP.workers)
    empty!(map_pid_wrkr)

    cluster_cookie(cookie)
    nothing
end

# The main function for adding worker processes.
# `manager` is of type ClusterManager. The respective managers are responsible
# for launching the workers. All keyword arguments (plus a few default values)
# are available as a dictionary to the `launch` methods
#
# Only one addprocs can be in progress at any time
#
const worker_lock = ReentrantLock()

"""
    addprocs(manager::ClusterManager; kwargs...) -> List of process identifiers

Launches worker processes via the specified cluster manager.

For example, Beowulf clusters are supported via a custom cluster manager implemented in
the package `ClusterManagers.jl`.

The number of seconds a newly launched worker waits for connection establishment from the
master can be specified via variable `JULIA_WORKER_TIMEOUT` in the worker process's
environment. Relevant only when using TCP/IP as transport.
"""
function addprocs(manager::ClusterManager; kwargs...)
    lock(worker_lock)
    try
        addprocs_locked(manager::ClusterManager; kwargs...)
    finally
        unlock(worker_lock)
    end
end

function addprocs_locked(manager::ClusterManager; kwargs...)
    params = merge(default_addprocs_params(), AnyDict(kwargs))
    topology(Symbol(params[:topology]))

    # References to launched workers, filled when each worker is fully initialized and
    # has connected to all nodes.
    launched_q = Int[]   # Asynchronously filled by the launch method

    # The `launch` method should add an object of type WorkerConfig for every
    # worker launched. It provides information required on how to connect
    # to it.
    launched = WorkerConfig[]
    launch_ntfy = Condition()

    # call manager's `launch` is a separate task. This allows the master
    # process initiate the connection setup process as and when workers come
    # online
    t_launch = @schedule launch(manager, params, launched, launch_ntfy)

    @sync begin
        while true
            if isempty(launched)
                istaskdone(t_launch) && break
                @schedule (sleep(1); notify(launch_ntfy))
                wait(launch_ntfy)
            end

            if !isempty(launched)
                wconfig = shift!(launched)
                let wconfig=wconfig
                    @async setup_launched_worker(manager, wconfig, launched_q)
                end
            end
        end
    end

    wait(t_launch)      # catches any thrown errors from the launch task

    # Since all worker-to-worker setups may not have completed by the time this
    # function returns to the caller, send the complete list to all workers.
    # Useful for nprocs(), nworkers(), etc to return valid values on the workers.
    all_w = workers()
    for pid in all_w
        remote_do(set_valid_processes, pid, all_w)
    end

    sort!(launched_q)
end

function set_valid_processes(plist::Array{Int})
    for pid in setdiff(plist, workers())
        myid() != pid && Worker(pid)
    end
end

default_addprocs_params() = AnyDict(
    :topology => :all_to_all,
    :dir      => pwd(),
    :exename  => joinpath(JULIA_HOME, julia_exename()),
    :exeflags => ``,
    :enable_threaded_blas => false)

function setup_launched_worker(manager, wconfig, launched_q)
    pid = create_worker(manager, wconfig)
    push!(launched_q, pid)

    # When starting workers on remote multi-core hosts, `launch` can (optionally) start only one
    # process on the remote machine, with a request to start additional workers of the
    # same type. This is done by setting an appropriate value to `WorkerConfig.cnt`.
    cnt = get(wconfig.count, 1)
    if cnt === :auto
        cnt = get(wconfig.environ)[:cpu_cores]
    end
    cnt = cnt - 1   # Removing self from the requested number

    if cnt > 0
        launch_n_additional_processes(manager, pid, wconfig, cnt, launched_q)
    end
end

function launch_n_additional_processes(manager, frompid, fromconfig, cnt, launched_q)
    @sync begin
        exename = get(fromconfig.exename)
        exeflags = get(fromconfig.exeflags, ``)
        cmd = `$exename $exeflags`

        new_addresses = remotecall_fetch(launch_additional, frompid, cnt, cmd)
        for address in new_addresses
            (bind_addr, port) = address

            wconfig = WorkerConfig()
            for x in [:host, :tunnel, :sshflags, :exeflags, :exename, :enable_threaded_blas]
                setfield!(wconfig, x, getfield(fromconfig, x))
            end
            wconfig.bind_addr = bind_addr
            wconfig.port = port

            let wconfig=wconfig
                @async begin
                    pid = create_worker(manager, wconfig)
                    remote_do(redirect_output_from_additional_worker, frompid, pid, port)
                    push!(launched_q, pid)
                end
            end
        end
    end
end

function create_worker(manager, wconfig)
    # only node 1 can add new nodes, since nobody else has the full list of address:port
    assert(LPROC.id == 1)

    # initiate a connect. Does not wait for connection completion in case of TCP.
    w = Worker()

    (r_s, w_s) = connect(manager, w.id, wconfig)
    w = Worker(w.id, r_s, w_s, manager; config=wconfig)
    # install a finalizer to perform cleanup if necessary
    finalizer(w, (w)->if myid() == 1 manage(w.manager, w.id, w.config, :finalize) end)

    # set when the new worker has finshed connections with all other workers
    ntfy_oid = RRID()
    rr_ntfy_join = lookup_ref(ntfy_oid)
    rr_ntfy_join.waitingfor = myid()

    # Start a new task to handle inbound messages from connected worker in master.
    # Also calls `wait_connected` on TCP streams.
    process_messages(w.r_stream, w.w_stream, false)

    # send address information of all workers to the new worker.
    # Cluster managers set the address of each worker in `WorkerConfig.connect_at`.
    # A new worker uses this to setup an all-to-all network if topology :all_to_all is specified.
    # Workers with higher pids connect to workers with lower pids. Except process 1 (master) which
    # initiates connections to all workers.

    # Connection Setup Protocol:
    # - Master sends 16-byte cookie followed by 16-byte version string and a JoinPGRP message to all workers
    # - On each worker
    #   - Worker responds with a 16-byte version followed by a JoinCompleteMsg
    #   - Connects to all workers less than its pid. Sends the cookie, version and an IdentifySocket message
    #   - Workers with incoming connection requests write back their Version and an IdentifySocketAckMsg message
    # - On master, receiving a JoinCompleteMsg triggers rr_ntfy_join (signifies that worker setup is complete)

    join_list = []
    if PGRP.topology == :all_to_all
        # need to wait for lower worker pids to have completed connecting, since the numerical value
        # of pids is relevant to the connection process, i.e., higher pids connect to lower pids and they
        # require the value of config.connect_at which is set only upon connection completion
        for jw in PGRP.workers
            if (jw.id != 1) && (jw.id < w.id)
                (jw.state == W_CREATED) && wait(jw.c_state)
                push!(join_list, jw)
            end
        end

    elseif PGRP.topology == :custom
        # wait for requested workers to be up before connecting to them.
        filterfunc(x) = (x.id != 1) && isdefined(x, :config) && (get(x.config.ident) in get(wconfig.connect_idents, []))

        wlist = filter(filterfunc, PGRP.workers)
        while length(wlist) < length(get(wconfig.connect_idents, []))
            sleep(1.0)
            wlist = filter(filterfunc, PGRP.workers)
        end

        for wl in wlist
            (wl.state == W_CREATED) && wait(wl.c_state)
            push!(join_list, wl)
        end
    end

    all_locs = map(x -> isa(x, Worker) ? (get(x.config.connect_at, ()), x.id) : ((), x.id, true), join_list)
    send_connection_hdr(w, true)
    join_message = JoinPGRPMsg(w.id, all_locs, PGRP.topology, get(wconfig.enable_threaded_blas, false))
    send_msg_now(w, MsgHeader(RRID(0,0), ntfy_oid), join_message)

    @schedule manage(w.manager, w.id, w.config, :register)
    wait(rr_ntfy_join)
    lock(client_refs) do
        delete!(PGRP.refs, ntfy_oid)
    end

    return w.id
end

# Called on the first worker on a remote host. Used to optimize launching
# of multiple workers on a remote host (to leverage multi-core)

additional_io_objs=Dict()
function launch_additional(np::Integer, cmd::Cmd)
    io_objs = Vector{Any}(np)
    addresses = Vector{Any}(np)

    for i in 1:np
        io, pobj = open(pipeline(detach(cmd), stderr=STDERR), "r")
        io_objs[i] = io
    end

    for (i,io) in enumerate(io_objs)
        (host, port) = read_worker_host_port(io)
        addresses[i] = (host, port)
        additional_io_objs[port] = io
    end

    return addresses
end

function redirect_output_from_additional_worker(pid, port)
    io = additional_io_objs[port]
    redirect_worker_output("$pid", io)
    delete!(additional_io_objs, port)
    nothing
end

function check_master_connect()
    timeout = worker_timeout()
    # If we do not have at least process 1 connect to us within timeout
    # we log an error and exit, unless we're running on valgrind
    if ccall(:jl_running_on_valgrind,Cint,()) != 0
        return
    end
    @schedule begin
        start = time()
        while !haskey(map_pid_wrkr, 1) && (time() - start) < timeout
            sleep(1.0)
        end

        if !haskey(map_pid_wrkr, 1)
            print(STDERR, "Master process (id 1) could not connect within $timeout seconds.\nexiting.\n")
            exit(1)
        end
    end
end

"""
    Base.cluster_cookie() -> cookie

Returns the cluster cookie.
"""
cluster_cookie() = LPROC.cookie

"""
    Base.cluster_cookie(cookie) -> cookie

Sets the passed cookie as the cluster cookie, then returns it.
"""
function cluster_cookie(cookie)
    # The cookie must be an ASCII string with length <=  HDR_COOKIE_LEN
    assert(isascii(cookie))
    assert(length(cookie) <= HDR_COOKIE_LEN)

    cookie = rpad(cookie, HDR_COOKIE_LEN)

    LPROC.cookie = cookie
    cookie
end

let next_pid = 2    # 1 is reserved for the client (always)
    global get_next_pid
    function get_next_pid()
        retval = next_pid
        next_pid += 1
        retval
    end
end

mutable struct ProcessGroup
    name::AbstractString
    workers::Array{Any,1}
    refs::Dict                  # global references
    topology::Symbol

    ProcessGroup(w::Array{Any,1}) = new("pg-default", w, Dict(), :all_to_all)
end
const PGRP = ProcessGroup([])

function topology(t)
    assert(t in [:all_to_all, :master_slave, :custom])
    if (PGRP.topology==t) || ((myid()==1) && (nprocs()==1)) || (myid() > 1)
        PGRP.topology = t
    else
        error("Workers with Topology $(PGRP.topology) already exist. Requested Topology $(t) cannot be set.")
    end
    t
end

get_bind_addr(pid::Integer) = get_bind_addr(worker_from_id(pid))
get_bind_addr(w::LocalProcess) = LPROC.bind_addr
function get_bind_addr(w::Worker)
    if isnull(w.config.bind_addr)
        if w.id != myid()
            w.config.bind_addr = remotecall_fetch(get_bind_addr, w.id, w.id)
        end
    end
    get(w.config.bind_addr)
end

# globals
const LPROC = LocalProcess()
const HDR_VERSION_LEN=16
const HDR_COOKIE_LEN=16
const map_pid_wrkr = Dict{Int, Union{Worker, LocalProcess}}()
const map_sock_wrkr = ObjectIdDict()
const map_del_wrkr = Set{Int}()

# cluster management related API
"""
    myid()

Get the id of the current process.
"""
myid() = LPROC.id

"""
    nprocs()

Get the number of available processes.
"""
nprocs() = length(PGRP.workers)

"""
    nworkers()

Get the number of available worker processes. This is one less than `nprocs()`. Equal to
`nprocs()` if `nprocs() == 1`.
"""
function nworkers()
    n = nprocs()
    n == 1 ? 1 : n-1
end

"""
    procs()

Returns a list of all process identifiers.
"""
procs() = Int[x.id for x in PGRP.workers]

"""
    procs(pid::Integer)

Returns a list of all process identifiers on the same physical node.
Specifically all workers bound to the same ip-address as `pid` are returned.
"""
function procs(pid::Integer)
    if myid() == 1
        if (pid == 1) || (isa(map_pid_wrkr[pid].manager, LocalManager))
            Int[x.id for x in filter(w -> (w.id==1) || (isa(w.manager, LocalManager)), PGRP.workers)]
        else
            ipatpid = get_bind_addr(pid)
            Int[x.id for x in filter(w -> get_bind_addr(w) == ipatpid, PGRP.workers)]
        end
    else
        remotecall_fetch(procs, 1, pid)
    end
end

"""
    workers()

Returns a list of all worker process identifiers.
"""
function workers()
    allp = procs()
    if nprocs() == 1
       allp
    else
       filter(x -> x != 1, allp)
    end
end

"""
    rmprocs(pids...; waitfor=typemax(Int))

Removes the specified workers. Note that only process 1 can add or remove
workers.

Argument `waitfor` specifies how long to wait for the workers to shut down:
    - If unspecified, `rmprocs` will wait until all requested `pids` are removed.
    - An `ErrorException` is raised if all workers cannot be terminated before
      the requested `waitfor` seconds.
    - With a `waitfor` value of 0, the call returns immediately with the workers
      scheduled for removal in a different task. The scheduled `Task` object is
      returned. The user should call `wait` on the task before invoking any other
      parallel calls.
"""
function rmprocs(pids...; waitfor=typemax(Int))
    # Only pid 1 can add and remove processes
    if myid() != 1
        throw(ErrorException("only process 1 can add and remove processes"))
    end

    pids = vcat(pids...)
    if waitfor == 0
        t = @schedule _rmprocs(pids, typemax(Int))
        yield()
        return t
    else
        _rmprocs(pids, waitfor)
        # return a dummy task object that user code can wait on.
        return @schedule nothing
    end
end

function _rmprocs(pids, waitfor)
    lock(worker_lock)
    try
        rmprocset = []
        for p in vcat(pids...)
            if p == 1
                warn("rmprocs: process 1 not removed")
            else
                if haskey(map_pid_wrkr, p)
                    w = map_pid_wrkr[p]
                    set_worker_state(w, W_TERMINATING)
                    kill(w.manager, p, w.config)
                    push!(rmprocset, w)
                end
            end
        end

        start = time()
        while (time() - start) < waitfor
            all(w -> w.state == W_TERMINATED, rmprocset) && break
            sleep(min(0.1, waitfor - (time() - start)))
        end

        unremoved = [wrkr.id for wrkr in filter(w -> w.state != W_TERMINATED, rmprocset)]
        if length(unremoved) > 0
            estr = string("rmprocs: pids ", unremoved, " not terminated after ", waitfor, " seconds.")
            throw(ErrorException(estr))
        end
    finally
        unlock(worker_lock)
    end
end

"""
    ProcessExitedException()

After a client Julia process has exited, further attempts to reference the dead child will
throw this exception.
"""
ProcessExitedException()
mutable struct ProcessExitedException <: Exception end

worker_from_id(i) = worker_from_id(PGRP, i)
function worker_from_id(pg::ProcessGroup, i)
    if in(i, map_del_wrkr)
        throw(ProcessExitedException())
    end
    if !haskey(map_pid_wrkr,i)
        if myid() == 1
            error("no process with id $i exists")
        end
        w = Worker(i)
        map_pid_wrkr[i] = w
    else
        w = map_pid_wrkr[i]
    end
    w
end

"""
    Base.worker_id_from_socket(s) -> pid

A low-level API which given a `IO` connection or a `Worker`,
returns the `pid` of the worker it is connected to.
This is useful when writing custom [`serialize`](@ref) methods for a type,
which optimizes the data written out depending on the receiving process id.
"""
function worker_id_from_socket(s)
    w = get(map_sock_wrkr, s, nothing)
    if isa(w,Worker)
        if s === w.r_stream || s === w.w_stream
            return w.id
        end
    end
    if isa(s,IOStream) && fd(s)==-1
        # serializing to a local buffer
        return myid()
    end
    return -1
end

register_worker(w) = register_worker(PGRP, w)
function register_worker(pg, w)
    push!(pg.workers, w)
    map_pid_wrkr[w.id] = w
end

function register_worker_streams(w)
    map_sock_wrkr[w.r_stream] = w
    map_sock_wrkr[w.w_stream] = w
end

deregister_worker(pid) = deregister_worker(PGRP, pid)
function deregister_worker(pg, pid)
    pg.workers = filter(x -> !(x.id == pid), pg.workers)
    w = pop!(map_pid_wrkr, pid, nothing)
    if isa(w, Worker)
        if isdefined(w, :r_stream)
            pop!(map_sock_wrkr, w.r_stream, nothing)
            if w.r_stream != w.w_stream
                pop!(map_sock_wrkr, w.w_stream, nothing)
            end
        end

        if myid() == 1
            # Notify the cluster manager of this workers death
            manage(w.manager, w.id, w.config, :deregister)
            if PGRP.topology != :all_to_all
                for rpid in workers()
                    try
                        remote_do(deregister_worker, rpid, pid)
                    catch
                    end
                end
            end
        end
    end
    push!(map_del_wrkr, pid)

    # delete this worker from our remote reference client sets
    ids = []
    tonotify = []
    lock(client_refs) do
        for (id,rv) in pg.refs
            if in(pid,rv.clientset)
                push!(ids, id)
            end
            if rv.waitingfor == pid
                push!(tonotify, (id,rv))
            end
        end
        for id in ids
            del_client(pg, id, pid)
        end

        # throw exception to tasks waiting for this pid
        for (id,rv) in tonotify
            notify_error(rv.c, ProcessExitedException())
            delete!(pg.refs, id)
        end
    end
end

function interrupt(pid::Integer)
    assert(myid() == 1)
    w = map_pid_wrkr[pid]
    if isa(w, Worker)
        manage(w.manager, w.id, w.config, :interrupt)
    end
end

"""
    interrupt(pids::Integer...)

Interrupt the current executing task on the specified workers. This is equivalent to
pressing Ctrl-C on the local machine. If no arguments are given, all workers are interrupted.
"""
interrupt(pids::Integer...) = interrupt([pids...])

"""
    interrupt(pids::AbstractVector=workers())

Interrupt the current executing task on the specified workers. This is equivalent to
pressing Ctrl-C on the local machine. If no arguments are given, all workers are interrupted.
"""
function interrupt(pids::AbstractVector=workers())
    assert(myid() == 1)
    @sync begin
        for pid in pids
            @async interrupt(pid)
        end
    end
end

function disable_nagle(sock)
    # disable nagle on all OSes
    ccall(:uv_tcp_nodelay, Cint, (Ptr{Void}, Cint), sock.handle, 1)
    @static if is_linux()
        # tcp_quickack is a linux only option
        if ccall(:jl_tcp_quickack, Cint, (Ptr{Void}, Cint), sock.handle, 1) < 0
            warn_once("Networking unoptimized ( Error enabling TCP_QUICKACK : ", Libc.strerror(Libc.errno()), " )")
        end
    end
end

function check_same_host(pids)
    if myid() != 1
        return remotecall_fetch(check_same_host, 1, pids)
    else
        # We checkfirst if all test pids have been started using the local manager,
        # else we check for the same bind_to addr. This handles the special case
        # where the local ip address may change - as during a system sleep/awake
        if all(p -> (p==1) || (isa(map_pid_wrkr[p].manager, LocalManager)), pids)
            return true
        else
            first_bind_addr = get(map_pid_wrkr[pids[1]].config.bind_addr)
            return all(p -> (p != 1) && (get(map_pid_wrkr[p].config.bind_addr) == first_bind_addr), pids[2:end])
        end
    end
end

function terminate_all_workers()
    myid() != 1 && return

    if nprocs() > 1
        try
            rmprocs(workers(); waitfor=5.0)
        catch _ex
            warn("Forcibly interrupting busy workers")
            # Might be computation bound, interrupt them and try again
            interrupt(workers())
            try
                rmprocs(workers(); waitfor=5.0)
            catch _ex2
                warn("Unable to terminate all workers")
            end
        end
    end
end

# initialize the local proc network address / port
function init_bind_addr()
    opts = JLOptions()
    if opts.bindto != C_NULL
        bind_to = split(unsafe_string(opts.bindto), ":")
        bind_addr = string(parse(IPAddr, bind_to[1]))
        if length(bind_to) > 1
            bind_port = parse(Int,bind_to[2])
        else
            bind_port = 0
        end
    else
        bind_port = 0
        try
            bind_addr = string(getipaddr())
        catch
            # All networking is unavailable, initialize bind_addr to the loopback address
            # Will cause an exception to be raised only when used.
            bind_addr = "127.0.0.1"
        end
    end
    global LPROC
    LPROC.bind_addr = bind_addr
    LPROC.bind_port = UInt16(bind_port)
end

function init_parallel()
    start_gc_msgs_task()
    atexit(terminate_all_workers)

    init_bind_addr()

    # start in "head node" mode, if worker, will override later.
    global PGRP
    global LPROC
    LPROC.id = 1
    cluster_cookie(randstring(HDR_COOKIE_LEN))
    assert(isempty(PGRP.workers))
    register_worker(LPROC)
end


text:
@warn("This script is not intended to be run, it is meant as a reference for further example scripts.")

# Auckland, with single replication.
# This is the same as the example script: JEMSS/example/example.jl.
# To animate (after sim is loaded), run: animate!(sim)
include("$(JEMSS.jemssDir)/example/cities/auckland/urban/1/single/example.jl")

# Auckland, with multiple replications.
# Multiple replications cannot be animated.
include("$(JEMSS.jemssDir)/example/cities/auckland/urban/1/reps/example.jl")

# Small and artificially generated example; sim loads relatively quickly.
include("$(JEMSS.jemssDir)/example/cities/generated/small/1/example.jl")

ext: jl


text:
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# speed2.jl
#
# Runs some JuMP benchmarks to test for speed-related regressions.
# Test time to build model in memory.
#############################################################################

using JuMP

function pMedian(numFacility::Int,numCustomer::Int,numLocation::Int)
    srand(10)
    customerLocations = [rand(1:numLocation) for a = 1:numCustomer ]

    tic()
    m = Model()

    # Facility locations
    @variable(m, 0 <= s[1:numLocation] <= 1)

    # Aux. variable: x_a,i = 1 iff the closest facility to a is at i
    @variable(m, 0 <= x[1:numLocation,1:numCustomer] <= 1)

    # Objective: min distance
    @objective(m, Max, sum{abs(customerLocations[a]-i)*x[i,a], a = 1:numCustomer, i = 1:numLocation} )

    # Constraints
    for a in 1:numCustomer
        # Subject to linking x with s
        for i in 1:numLocation
            @constraint(m, x[i,a] - s[i] <= 0)
        end
        # Subject to one of x must be 1
        @constraint(m, sum{x[i,a],i=1:numLocation} == 1 )
    end

    # Subject to must allocate all facilities
    @constraint(m, sum{s[i],i=1:numLocation} == numFacility )
    buildTime = toq()

    tic()
    JuMP.build(m)
    writeTime = toq()

    return buildTime, writeTime
end

function cont5(n)
    m = n
    n1 = n-1
    m1 = m-1
    dx = 1/n
    T = 1.58
    dt = T/m
    h2 = dx^2
    a = 0.001
    yt = [0.5*(1 - (j*dx)^2) for j=0:n]

    tic()
    mod = Model()
    @variable(mod,  0 <= y[0:m,0:n] <= 1)
    @variable(mod, -1 <= u[1:m] <= 1)
    @objective(mod, Min, 0.25*dx*( (y[m,0] - yt[1])^2 +
       2*sum{ (y[m,j]-yt[j+1])^2, j=1:n1} + (y[m,n]-yt[n+1])^2) +
       0.25*a*dt*(2*sum{u[i]^2,i=1:m1} + u[m]^2))

    # PDE
    for i = 0:m1
        for j = 1:n1
            @constraint(mod, h2*(y[i+1,j] - y[i,j]) == 0.5*dt*(y[i,j-1] - 2*y[i,j] + y[i,j+1] + y[i+1,j-1] - 2*y[i+1,j] + y[i+1,j+1]) )
        end
    end

    # IC
    for j = 0:n
        @constraint(mod, y[0,j] == 0)
    end

    # BC
    for i = 1:m
        @constraint(mod, y[i,2]   - 4*y[i,1]  + 3*y[i,0] == 0)
        @constraint(mod, y[i,n-2] - 4*y[i,n1] + 3*y[i,n] == (2*dx)*(u[i] - y[i,n]))
    end
    buildTime = toq()

    tic()
    JuMP.build(mod)
    writeTime = toq()

    return buildTime, writeTime
end

function RunTests()
    # Pmedian
    pmedian_build = Float64[]
    pmedian_write = Float64[]
    for runs = 1:9
        bt, wt = pMedian(100,100,5000)
        push!(pmedian_build, bt)
        push!(pmedian_write, wt)
    end
    sort!(pmedian_build)
    sort!(pmedian_write)
    print("PMEDIAN BUILD MIN=",minimum(pmedian_build),"  MED=",pmedian_build[5],"\n")
    print("PMEDIAN INTRN MIN=",minimum(pmedian_write),"  MED=",pmedian_write[5],"\n")

    # Cont5
    cont5_build = Float64[]
    cont5_write = Float64[]
    for runs = 1:9
        bt, wt = cont5(500)
        push!(cont5_build, bt)
        push!(cont5_write, wt)
    end
    sort!(cont5_build)
    sort!(cont5_write)
    print("CONT5 BUILD   MIN=",minimum(cont5_build),"  MED=",cont5_build[5],"\n")
    print("CONT5 INTRN   MIN=",minimum(cont5_write),"  MED=",cont5_write[5],"\n")

end

RunTests()


text:
# ~ bound checking functions ~
checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Integer) =
    (0 < k ≤ size(A, 1) && 0 < j ≤ size(A, 2) || throw(BoundsError(A, (k,j))))

checkbounds(A::AbstractDerivativeOperator, kr::AbstractRange, j::Integer) =
    (checkbounds(A, first(kr), j); checkbounds(A,  last(kr), j))

checkbounds(A::AbstractDerivativeOperator, k::Integer, jr::AbstractRange) =
    (checkbounds(A, k, first(jr)); checkbounds(A, k,  last(jr)))

checkbounds(A::AbstractDerivativeOperator, kr::AbstractRange, jr::AbstractRange) =
    (checkbounds(A, kr, first(jr)); checkbounds(A, kr,  last(jr)))

checkbounds(A::AbstractDerivativeOperator, k::Colon, j::Integer) =
    (0 < j ≤ size(A, 2) || throw(BoundsError(A, (size(A,1),j))))

checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Colon) =
    (0 < k ≤ size(A, 1) || throw(BoundsError(A, (k,size(A,2)))))

@inline function getindex(A::AbstractDerivativeOperator, i::Int, j::Int)
    @boundscheck checkbounds(A, i, j)
    # TODO: Implement the lazy version
    if use_winding(A)
        return Array(A)[i,j]
    end

    bpc = A.boundary_point_count
    N = A.len
    bsl = A.boundary_stencil_length
    slen = A.stencil_length
    if bpc > 0 && 1<=i<=bpc
        if j > bsl
            return 0
        else
            return A.low_boundary_coefs[i][j]
        end
    elseif bpc > 0 && (N-bpc)<i<=N
        if j < N+2-bsl
            return 0
        else
            return A.high_boundary_coefs[bpc-(N-i)][bsl-(N+2-j)]
        end
    else
        if j < i-bpc || j > i+slen-bpc-1
            return 0
        else
            return A.stencil_coefs[j-i + 1 + bpc]
        end
    end
end

# scalar - colon - colon
@inline getindex(A::AbstractDerivativeOperator, ::Colon, ::Colon) = Array(A)

@inline function getindex(A::AbstractDerivativeOperator, ::Colon, j)
    T = eltype(A.stencil_coefs)
    v = zeros(T, A.len)
    v[j] = one(T)
    copyto!(v, A*v)
    return v
end

@inline function getindex(A::AbstractDerivativeOperator, i, ::Colon)
    @boundscheck checkbounds(A, i, 1)
    T = eltype(A.stencil_coefs)
    v = zeros(T, A.len+2)

    bpc = A.boundary_point_count
    N = A.len
    bsl = A.boundary_stencil_length
    slen = A.stencil_length

    if bpc > 0 && 1<=i<=bpc
        v[1:bsl] .= A.low_boundary_coefs[i]
    elseif bpc > 0 && (N-bpc)<i<=N
         v[1:bsl]  .= A.high_boundary_coefs[i-(N-1)]
    else
        if use_winding(A)
            v[i-bpc+slen:i-bpc+2*slen-1] .= A.stencil_coefs
        else
            v[i-bpc:i-bpc+slen-1] .= A.stencil_coefs
        end
    end
    return v
end

# UnitRanges
@inline function getindex(A::AbstractDerivativeOperator, rng::UnitRange{Int}, ::Colon)
    m = BandedMatrix(A)
    return m[rng, cc]
end

@inline function getindex(A::AbstractDerivativeOperator, ::Colon, rng::UnitRange{Int})
    m = BandedMatrix(A)
    return m[rnd, cc]
end

@inline function getindex(A::AbstractDerivativeOperator, r::Int, rng::UnitRange{Int})
    m = A[r, :]
    return m[rng]
end

@inline function getindex(A::AbstractDerivativeOperator, rng::UnitRange{Int}, c::Int)
    m = A[:, c]
    return m[rng]
end

@inline function getindex(A::AbstractDerivativeOperator{T}, rng::UnitRange{Int}, cng::UnitRange{Int}) where T
    return BandedMatrix(A)[rng,cng]
end

# Base.length(A::AbstractDerivativeOperator) = A.stencil_length
Base.ndims(A::AbstractDerivativeOperator) = 2
Base.size(A::AbstractDerivativeOperator) = (A.len, A.len + 2)
Base.size(A::AbstractDerivativeOperator,i::Integer) = size(A)[i]
Base.length(A::AbstractDerivativeOperator) = reduce(*, size(A))

#=
    For the evenly spaced grid we have a symmetric matrix
=#
Base.transpose(A::DerivativeOperator) = A
Base.adjoint(A::DerivativeOperator) = A
LinearAlgebra.issymmetric(::DerivativeOperator) = true

#=
    Fallback methods that use the full representation of the operator
=#
Base.exp(A::AbstractDerivativeOperator{T}) where T = exp(convert(A))
Base.:\(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A \ convert(Array,B)
Base.:\(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = Array(A) \ B
Base.:/(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A / convert(Array,B)
Base.:/(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = Array(A) / B

#=
    The Inf opnorm can be calculated easily using the stencil coeffiicents, while other opnorms
    default to compute from the full matrix form.
=#
function LinearAlgebra.opnorm(A::DerivativeOperator, p::Real=2)
    if p == Inf
        sum(abs.(A.stencil_coefs)) / A.dx^A.derivative_order
    else
        opnorm(BandedMatrix(A), p)
    end
end

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

get_type(::AbstractDerivativeOperator{T}) where {T} = T

function *(M::AbstractMatrix,A::AbstractDerivativeOperator)
    y = zeros(promote_type(eltype(A),eltype(M)), size(M,1), size(A,2))
    LinearAlgebra.mul!(y, M, BandedMatrix(A))
    return y
end

function *(A::AbstractDerivativeOperator,B::AbstractDerivativeOperator)
    return BandedMatrix(A)*BandedMatrix(B)
end

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

function *(coeff_func::Function, A::DerivativeOperator{T,N,Wind}) where {T,N,Wind}
    coefficients = A.coefficients === nothing ? Vector{T}(undef,A.len) : A.coefficients
    DerivativeOperator{T,N,Wind,typeof(A.dx),typeof(A.stencil_coefs),
                       typeof(A.low_boundary_coefs),typeof(coefficients),
                       typeof(coeff_func)}(
        A.derivative_order, A.approximation_order,
        A.dx, A.len, A.stencil_length,
        A.stencil_coefs,
        A.boundary_stencil_length,
        A.boundary_point_count,
        A.low_boundary_coefs,
        A.high_boundary_coefs,coefficients,coeff_func
        )
end

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

function DiffEqBase.update_coefficients!(A::AbstractDerivativeOperator,u,p,t)
    if A.coeff_func !== nothing
        A.coeff_func(A.coefficients,u,p,t)
    end
end

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

(L::DerivativeOperator)(u,p,t) = L*u
(L::DerivativeOperator)(du,u,p,t) = mul!(du,L,u)


text:
using RegisterPOI
using Test


text:
@time using Dubins
@time using Dubins
@time using Dubins
# using Makie

@time d = Dubins.dubins_shortest_path([0., 0., 1.], [1., 1., 2.], 20.)
@time d = Dubins.dubins_shortest_path([0., 0., 1.], [1., 1., 2.], 20.)
@time d = Dubins.dubins_shortest_path([0., 0., 1.], [1., 1., 2.], 20.)


text:
module BeerTest

include("main.jl")

export main, initialize

end


text:
module CLIFramework

export CommandBundle, add!, dispatch

type CommandBundle
    commands::Dict{AbstractString, Function}
    main::Function
    CommandBundle(f::Function) = new(Dict{AbstractArray, Function}(), f)
end

"""
Assign a function to a command (will be passed the remainter of ARGS
    as an argument)
"""
add!(f::Function, cb::CommandBundle, key::AbstractString) = cb.commands[key] = f

function dispatch(cb::CommandBundle)
    if isempty(ARGS)
        cb.main(ARGS)
    elseif ARGS[1] in keys(cb.commands)
        cb.commands[ARGS[1]](ARGS[2:end])
    else
        cb.main(ARGS)
    end
end

end  # module CLIFramework


text:
#=
This code can be used to replicate the test results from
A.N. Riseth - "Objective acceleration for unconstrained optimization", 2017.
https://arxiv.org/abs/1710.05200
=#
module OACCEL2017

using OptimTestProblems.MultivariateProblems
using AccelerationBenchmark
using Optim, LineSearches
using JLD, CSV

const MVP = MultivariateProblems
const UP = MVP.UnconstrainedProblems

function defaulttestsetup()
    # Line search parameters taken from De Sterck
    # MATLAB files for Poblano toolbox.
    # The De Sterck gtol=1e-2 is *very* small. Typical values are 1e-1 for CG methods and 0.9 for quasi-Newton
    # In the updated paper we use gtol=0.1.
    ls = MoreThuente(f_tol = 1e-4, gtol = 1e-1, x_tol = 1e-15,
                     alphamin = 1e-15, alphamax = 1e15, maxfev = 20)
    lsstatic = Static(alpha = 1.0)
    agscaled = InitialStatic(alpha = 1e-4, scaled = true)
    ag = InitialStatic(alpha=1.0, scaled=false)

    gdls = GradientDescent(alphaguess = ag, linesearch = ls)
    gdst = GradientDescent(alphaguess = agscaled, linesearch = lsstatic)

    solvers = [
        LBFGS(alphaguess = ag, linesearch = ls, m = 5, scaleinvH0=true),
        NGMRES(alphaguess = ag, linesearch = ls, wmax = 20, nlprecon = gdls),
        NGMRES(alphaguess = ag, linesearch = ls, wmax = 20, nlprecon = gdst),
        OACCEL(alphaguess = ag, linesearch = ls, wmax = 20, nlprecon = gdls),
        OACCEL(alphaguess = ag, linesearch = ls, wmax = 20, nlprecon = gdst),
    ]
    solvernames = ["L-BFGS",
                   "N-GMRES-A", "N-GMRES-B",
                   "O-ACCEL-A", "O-ACCEL-B",
                   ]
    stoptype = :FunctionTolRelative
    stoptol  = 1e-10
    #stoptype = :GradientTolRelative
    #stoptol  = 1e-8
    timelog  = 0 # timelog = 1 can cause very inaccurate timings due to compilation and garbage collection
    maxiter  = 2000
    timelimit = NaN
    TestSetup(solvers,solvernames,stoptype,stoptol,timelog,maxiter,timelimit)
end

"Problem A in  AN Riseth - *Objective acceleration for unconstrained optimization*, 2017."
function createA(N::Int, seednum::Int)
    probname = "A-$N-$seednum"
    prob = UP._quadraticproblem(N)

    rnd = RandomizeInitialx([seednum])
    randomizeproblem!(prob, rnd, 1)
    return prob, probname
end

"Problem B in  AN Riseth - *Objective acceleration for unconstrained optimization*, 2017."
function createB(N::Int, seednum::Int)
    probname = "B-$N-$seednum"
    prob = UP._paraboloidproblem(N)

    rnd = RandomizeInitialx([seednum])
    randomizeproblem!(prob, rnd, 1)
    return prob, probname
end

"Problem C in  AN Riseth - *Objective acceleration for unconstrained optimization*, 2017."
function createC(N::Int, seednum::Int)
    probname = "C-$N-$seednum"
    prob = UP._paraboloidproblem(N; mat=eye(N))

    rnd = RandomizeInitialxMat([seednum])
    randomizeproblem!(prob, rnd, 1)
    return prob, probname
end

"Problem D in  AN Riseth - *Objective acceleration for unconstrained optimization*, 2017."
function createD(N::Int, seednum::Int)
    probname = "D-$N-$seednum"
    prob = UP._extrosenbrockproblem(N)

    rnd = RandomizeInitialx([seednum])
    randomizeproblem!(prob, rnd, 1)
    return prob, probname
end

"Problem E in  AN Riseth - *Objective acceleration for unconstrained optimization*, 2017."
function createE(N::Int, seednum::Int)
    probname = "E-$N-$seednum"
    prob = UP._extpowellproblem(N)

    rnd = RandomizeInitialx([seednum])
    randomizeproblem!(prob, rnd, 1)
    return prob, probname
end

"""
Problem F in  AN Riseth - *Objective acceleration for unconstrained optimization*, 2017.

WARNING: `_trigonometricproblem` follows Moré et al. - *Testing unconstrained optimization software*, 1981,
which has a sign difference compared to Riseth.
"""
function createF(N::Int, seednum::Int)
    probname = "F-$N-$seednum"
    prob = UP._trigonometricproblem(N)

    rnd = RandomizeInitialx([seednum])
    randomizeproblem!(prob, rnd, 1)
    return prob, probname
end

"""
Problem G in  AN Riseth - *Objective acceleration for unconstrained optimization*, 2017.

WARNING: This function does not have a closed-form minimum value,
and the minimum depends on `N`.
"""
function createG(N::Int, seednum::Int)
    probname = "G-$N-$seednum"
    prob = UP._penfunIproblem(N)
    rnd = RandomizeInitialx([seednum])
    randomizeproblem!(prob, rnd, 1)
    return prob, probname
end

macro createproblem(fun, N, seednum)
    eval(Symbol(:create,fun))(eval(N),eval(seednum))
end

"Helper function for pmap in createrunsrandomized."
function getoruns(seed::Int, N::Int, cfun, ts::TestSetup, showname::Bool)
    prob, probname = cfun(N,seed)
    showname && println("$probname")
    createruns(prob, probname,
               ts.solvers, ts.solvernames,
               ts.stoptol, ts.stoptype, ts.timelog, ts.maxiter;
               time_limit = ts.timelimit)
end

"""
Run the same problem family multiple times with different randomized conditions.

Return a vector or `OptimizationRun`s
"""
function createrunsrandomized(fun::Symbol, N::Int, seeds::AbstractArray,
                              ts::TestSetup = defaulttestsetup();
                              showname::Bool = false)
    cfun = eval(Symbol(:create,fun))
    np = length(seeds)

    oruns = reduce(append!, pmap(k->getoruns(seeds[k],N,cfun,ts,showname), 1:np))

    return oruns
end

"""
Run problems and create tables with the corresponding measures, to be used
to investigate performance and plot performance profiles.
"""
function runmany(funs::AbstractVector, seeds::AbstractArray{Int},
                 ts::TestSetup = defaulttestsetup();
                 savejld::Bool = true,
                 savecsv::Bool = true,
                 savebase::AbstractString = "data/",
                 showname::Bool = true)
    # funs is a collection of Tuple{Symbol,Int}
    @assert savejld || savecsv
    if myid() == 1 && !isdir(savebase)
        mkdir(savebase)
    end
    for fN in funs
        oruns = createrunsrandomized(fN..., seeds, ts; showname=showname)
        mdf = createmeasuredataframe(oruns)
        if savejld
            myid() == 1 && save(savebase*"/$(fN[1])-$(fN[2]).jld", "mdf", mdf)
        end
        if savecsv
            myid() == 1 && CSV.write(savebase*"/$(fN[1])-$(fN[2]).csv", mdf)
        end
    end
end

end


text:
export CircularArrayBuffer, capacity, isfull

using ReinforcementLearningBase

"""
    CircularArrayBuffer{T}(d::Integer...) -> CircularArrayBuffer{T, N}

`CircularArrayBuffer` uses a `N`-dimention `Array` of size `d` to serve as a buffer for
`N-1`-dimention `Array`s with the same size.

# Examples

```julia-repl
julia> b = CircularArrayBuffer{Float64}(2, 2, 3)
2×2×0 CircularArrayBuffer{Float64,3}

julia> capacity(b)
3

julia> length(b)
0

julia> push!(b, [1. 1.; 2. 2.])
2×2×1 CircularArrayBuffer{Float64,3}:
[:, :, 1] =
 1.0  1.0
 2.0  2.0

julia> b
2×2×1 CircularArrayBuffer{Float64,3}:
[:, :, 1] =
 1.0  1.0
 2.0  2.0

julia> length(b)
4

julia> nframes(cb::CircularArrayBuffer) = cb.length
nframes (generic function with 1 method)

julia> nframes(b)
1

julia> ones(2,2)
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

julia> 3 .* ones(2,2)
2×2 Array{Float64,2}:
 3.0  3.0
 3.0  3.0

julia> 3 * ones(2,2)
2×2 Array{Float64,2}:
 3.0  3.0
 3.0  3.0

julia> b = CircularArrayBuffer{Float64}(2, 2, 3)
2×2×0 CircularArrayBuffer{Float64,3}

julia> capacity(b)
3

julia> nframes(b)
0

julia> push!(b, 1 .* ones(2,2))
2×2×1 CircularArrayBuffer{Float64,3}:
[:, :, 1] =
 1.0  1.0
 1.0  1.0

julia> b
2×2×1 CircularArrayBuffer{Float64,3}:
[:, :, 1] =
 1.0  1.0
 1.0  1.0

julia> nframes(b)
1

julia> for i in 2:4
           push!(b, i .* ones(2,2))
       end

julia> b
2×2×3 CircularArrayBuffer{Float64,3}:
[:, :, 1] =
 2.0  2.0
 2.0  2.0

[:, :, 2] =
 3.0  3.0
 3.0  3.0

[:, :, 3] =
 4.0  4.0
 4.0  4.0

julia> isfull(b)
true

julia> nframes(b)
3

julia> size(b)
(2, 2, 3)
```
"""
mutable struct CircularArrayBuffer{T,N} <: AbstractArray{T,N}
    buffer::Array{T,N}
    first::Int
    length::Int
    step_size::Int

    function CircularArrayBuffer{T}(d::Integer...) where {T}
        N = length(d)
        new{T,N}(Array{T}(undef, d...), 1, 0, N == 1 ? 1 : *(d[1:end-1]...))
    end
end

Base.IndexStyle(::CircularArrayBuffer) = IndexLinear()
Base.size(cb::CircularArrayBuffer{<:Any,N}, i::Integer) where {N} =
    i == N ? cb.length : size(cb.buffer, i)
Base.size(cb::CircularArrayBuffer{<:Any,N}) where {N} = ntuple(M -> size(cb, M), N)
Base.getindex(cb::CircularArrayBuffer{T,N}, i::Int) where {T,N} =
    getindex(cb.buffer, _buffer_index(cb, i))
Base.setindex!(cb::CircularArrayBuffer{T,N}, v, i::Int) where {T,N} =
    setindex!(cb.buffer, v, _buffer_index(cb, i))
capacity(cb::CircularArrayBuffer{T,N}) where {T,N} = size(cb.buffer, N)
isfull(cb::CircularArrayBuffer) = cb.length == capacity(cb)
Base.isempty(cb::CircularArrayBuffer) = cb.length == 0

@inline function _buffer_index(cb::CircularArrayBuffer, i::Int)
    ind = (cb.first - 1) * cb.step_size + i
    if ind > length(cb.buffer)
        ind - length(cb.buffer)
    else
        ind
    end
end

@inline function _buffer_frame(cb::CircularArrayBuffer, i::Int)
    n = capacity(cb)
    idx = cb.first + i - 1
    if idx > n
        idx - n
    else
        idx
    end
end

_buffer_frame(cb::CircularArrayBuffer, I::Vector{Int}) = map(i -> _buffer_frame(cb, i), I)

function Base.empty!(cb::CircularArrayBuffer)
    cb.first = 1
    cb.length = 0
    cb
end

"""
    update!(cb::CircularArrayBuffer{T,N}, data::AbstractArray)

`update!` the last frame of `cb` with data.
"""
function RLBase.update!(cb::CircularArrayBuffer{T,N}, data::AbstractArray) where {T,N}
    select_last_dim(cb.buffer, _buffer_frame(cb, cb.length)) .= data
    cb
end

function RLBase.update!(cb::CircularArrayBuffer{T,1}, data) where {T}
    cb.buffer[_buffer_frame(cb, cb.length)] = data
    cb
end

function Base.push!(cb::CircularArrayBuffer, data)
    # length(data) == cb.step_size || throw(ArgumentError("length of , $(cb.step_size) != $(length(data))"))
    push!(cb, missing)
    update!(cb, data)
    cb
end

function Base.push!(cb::CircularArrayBuffer{T,N}, ::Missing) where {T,N}
    if cb.length == capacity(cb)
        cb.first = (cb.first == capacity(cb) ? 1 : cb.first + 1)
    else
        cb.length += 1
    end
    cb
end

function Base.pop!(cb::CircularArrayBuffer)
    res = select_last_frame(cb)
    if cb.length <= 0
        throw(ArgumentError("buffer must be non-empty"))
    else
        cb.length -= 1
    end
    res
end

frame_type(::CircularArrayBuffer{T,N}) where {T,N} = Array{T,N - 1}
frame_type(::CircularArrayBuffer{T,1}) where {T} = T


text:
function gen_newton_functions(solver::Solver)
    n,m,N = get_sizes(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    c_function!, c_jacobian!, c_labels, cI!, cE! = generate_constraint_functions(solver.obj)
    f = solver.fd
    s_x = reshape(1:(N*n),n,N)
    s_u = reshape(N*n .+ (1:(m*(N-1))),m,N-1)
    s_λ = (N*n + (N-1)*m) .+ (1:((N-1)*pE + pE_N))
    s_μ = (N*n + (N-1)*m + length(s_λ)) .+ (1:((N-1)*pI + pI_N))
    s_ν = (N*n + (N-1)*m + length(s_λ) + length(s_μ)) .+ (1:N*n)
    s_s = (N*n + (N-1)*m + length(s_λ) + length(s_μ) + N*n) .+ (1:length(s_μ))
    sze = (X=s_x, U=s_u, λ=s_λ, μ=s_μ, ν=s_ν, s=s_s)

    function newton_cost(V::Vector)
        X,U,λ,μ,ν,s = unpackV(V)
        J = cost(solver,X,U) + λ'cE(V) + μ'*(cI(V) + 0.5*s.^2) + ν'd(V)
        return J
    end

    function cI(V)
        X,U,λ,μ,ν,s = unpackV(V)
        C = zeros(eltype(V),pI,N-1)
        for k = 1:N-1
            cI!(view(C,1:pI,k),X[:,k],U[:,k])
        end
        CN = zeros(eltype(V),pI_N)
        cI!(CN,X[:,N])
        CI = [vec(C); CN]
        return CI
    end

    function cE(V)
        X,U,λ,μ,ν,s = unpackV(V)
        C = zeros(eltype(V),pE,N-1)
        for k = 1:N-1
            cE!(view(C,1:pE,k),X[:,k],U[:,k])
        end
        CN = zeros(eltype(V),pE_N)
        cE!(CN,X[:,N])
        CE = [vec(C); CN]
        return CE
    end

    function d(V)
        X,U,λ,μ,ν,s = unpackV(V)
        D = zeros(eltype(V),n,N)
        D[:,1] = X[:,1] - solver.obj.x0
        for k = 2:N
            f(view(D,1:n,k),X[:,k-1],U[:,k-1])
            D[:,k] -= X[:,k]
        end
        return vec(D)
    end

    function active_set(V,eps=1e-2)
        X,U,λ,μ,ν,s = unpackV(V)
        a = ones(Bool,length(V))
        ci = cI(V)
        c_inds = ones(Bool,length(ci))
        c_inds = ci .>= -eps

        a[sze.μ] = c_inds
        a[sze.s] = c_inds
        return a
    end

    function packV(X,U,λ,μ,ν,s)
        V = [vec(X);vec(U);λ;μ;ν;s]
    end

    function unpackV(V)
        X = V[sze.X]
        U = V[sze.U]
        λ = V[sze.λ]
        μ = V[sze.μ]
        ν = V[sze.ν]
        s = V[sze.s]
        return X,U,λ,μ,ν,s
    end
    return newton_cost, packV, unpackV, cI, cE, d, active_set, sze
end

model,obj = Dynamics.dubinscar_parallelpark
solver = Solver(model,obj,N=5)
solver.opts.constraint_tolerance = 1e-4
n,m,N = get_sizes(solver)
U0 = rand(m,N-1)
res,stats = solve(solver,U0)
cost(solver,res)
λ_update_default!(res,solver);
@assert max_violation(res) < solver.opts.constraint_tolerance
max_violation(res)

newton_cost, packV, unpackV, cI, cE, d, active_set = gen_newton_functions(solver)
max_c(V) = max(maximum(cI(V)),norm(cE(V),Inf),norm(d(V),Inf))
p,pI,pE = get_num_constraints(solver)
pN,pI_N,pE_N = get_num_terminal_constraints(solver)

X = to_array(res.X)
U = to_array(res.U)
μ = vcat([k < N ? res.λ[k][1:pI] : res.λ[k][1:pI_N] for k = 1:N]...)
λ = vcat([k < N ? res.λ[k][pI .+ (1:pE)] : res.λ[k][pI_N .+ (1:pE_N)] for k = 1:N]...)
ν = vec(to_array(res.s))
s = vcat([k < N ? res.C[k][1:pI] : res.C[k][1:pI_N] for k = 1:N]...)
s = sqrt.(2*max.(0,-s))

V = packV(X,U,λ,μ,ν,s)

J1 = newton_cost(V)
g = ForwardDiff.gradient(newton_cost,V)
H = ForwardDiff.hessian(newton_cost,V)
# V = line_search(V,H,g)

a = active_set(V,1e-3)
rank(H)
size(H)
rank(H[a,a])
size(H[a,a])
ρ = 1.0
newton_cost(V - (H + ρ*I)\g)
max_c(V - (H + ρ*I)\g)

X,U,λ,μ,ν,s == unpackV(V)
ν
# newton_cost(V)
# d(V)
# norm(d(V),Inf)
# X,U,λ,μ,ν,s == unpackV(V)
# d(V)
#
# function line_search(V,H,g)
#     J0 = newton_cost(V)
#     α = 1.
#     V_ = V - α*(H\g)
#     J = newton_cost(V_)
#     iter = 1
#     while J > J0
#         α /= 2
#         V_ = V - α*(H\g)
#         J = newton_cost(V_)
#         iter += 1
#         @show J, α
#         iter > 10 ? break : nothing
#     end
#     @show J0-J
#     return V_
# end

max_c(V)
findmax(cI(V))
cI(V)
X,U,λ,μ,ν,s == unpackV(V)
μ[8]
norm(μ.*s,Inf)
d(V)

# Block move
opts = TrajectoryOptimization.SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-4
opts.cost_tolerance_intermediate = 1e-4
opts.constraint_tolerance = 1e-4
opts.square_root = false
opts.active_constraint_tolerance = 0.0
opts.outer_loop_update_type = :default
opts.live_plotting = false

# Block move
model, obj = TrajectoryOptimization.Dynamics.double_integrator
u_min = -0.2
u_max = 0.2
obj_con = TrajectoryOptimization.ConstrainedObjective(obj, tf=5.0,use_xf_equality_constraint=true,u_min=u_min, u_max=u_max)#, x_max=x_max)

solver = TrajectoryOptimization.Solver(model,obj_con,integration=:rk4,N=5,opts=opts)
U = rand(solver.model.m, solver.N)

res, stats = TrajectoryOptimization.solve(solver,U)
@assert max_violation(res) < opts.constraint_tolerance
max_violation(res)

newton_cost, packV, unpackV, cI, cE, d, active_set = gen_newton_functions(solver)
max_c(V) = max(maximum(cI(V)),norm(cE(V),Inf),norm(d(V),Inf))
p,pI,pE = get_num_constraints(solver)
pN,pI_N,pE_N = get_num_terminal_constraints(solver)

X = to_array(res.X)
U = to_array(res.U)
μ = vec(to_array(res.λ[1:solver.N-1]))
λ = vec(res.λ[solver.N])
ν = vec(to_array(res.s))
s = vcat([k < N ? res.C[k][1:pI] : nothing for k = 1:solver.N-1]...)
s = sqrt.(2*max.(0,-s))
V = packV(X,U,λ,μ,ν,s)
cost(solver,res)
max_c(V)
newton_cost(V)
g = ForwardDiff.gradient(newton_cost,V)
H = ForwardDiff.hessian(newton_cost,V)
a = active_set(V,1.0)
rank(H[a,a])
size(H[a,a])
V[a] = V[a] - (H[a,a] + 0I)\g[a]
# V = V - H\g
newton_cost(V)
max_c(V)

X,U,λ,μ,ν,s == unpackV(V)
ν

d(V)


text:
"""
    metropolis!(model::Ising)

Implementation of the Metropolis algorithm for the Ising model.

Arguments
* `model`: Ising type
"""
function metropolis!(model::Ising)
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)
    ΔE = [-8, -4, 0, 4, 8]
    P = exp(-model.params.β) .^ ΔE

    for t = 1:model.params.n_sweeps
        for j = 1:model.params.L, i = 1:model.params.L
            k = Int(model.σ[i, j] * sum_of_neighbors(model, (i, j)) / 2 + 3)

            (ΔE[k] <= 0 || rand(model.rng) < P[k]) && (model.σ[i, j] *= -1)
        end

        t > t₀ && update_observables!(model)
    end

    compute_observables_statistics!(model)
    return model
end

"""
    metropolis!(model::Potts)

Implementation of the Metropolis algorithm for the Potts model.

Arguments
* `model`: Potts type
"""
function metropolis!(model::Potts)
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)
    ΔE = collect(-4:4)
    P = exp(-model.params.β) .^ ΔE

    for t = 1:model.params.n_sweeps
        for j = 1:model.params.L, i = 1:model.params.L
            counts = compute_counts(model, (i, j))
            old_spin = model.σ[i, j]
            new_spin = rand(model.rng, UnitRange{Int8}(1, model.params.q))
            k = (counts[old_spin] - counts[new_spin]) + 5

            (ΔE[k] <= 0 || rand(model.rng) < P[k]) && (model.σ[i, j] = new_spin)
        end

        t > t₀ && update_observables!(model)
    end

    compute_observables_statistics!(model)
    return model
end

"""
    metropolis!(model::XY)

Implementation of the Metropolis algorithm for the XY model.

Arguments
* `model`: XY type
"""
function metropolis!(model::XY)
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)
    eᵝ = exp(-model.params.β)

    for t = 1:model.params.n_sweeps
        for j = 1:model.params.L, i = 1:model.params.L
            old_spin = model.σ[i, j]
            new_spin = random_XYVector(model.rng)
            ΔE = dot(old_spin .- new_spin, sum_of_neighbors(model, (i, j)))

            (ΔE <= 0 || rand(model.rng) < eᵝ^ΔE) && (model.σ[i, j] = new_spin)
        end

        t > t₀ && update_observables!(model)
    end

    compute_observables_statistics!(model)
    return model
end

"""
    heat_bath!(model::Ising)

Implementation of the heat bath algorithm for the Ising model.

Arguments
* `model`: Ising type
"""
function heat_bath!(model::Ising)
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)
    ΔE = [-8, -4, 0, 4, 8]
    P = exp(-model.params.β) .^ ΔE

    for t = 1:model.params.n_sweeps
        for j = 1:model.params.L, i = 1:model.params.L
            k = Int(model.σ[i, j] * sum_of_neighbors(model, (i, j)) / 2 + 3)

            (rand(model.rng) < P[k] / (P[k] + 1)) && (model.σ[i, j] *= -1)
        end

        t > t₀ && update_observables!(model)
    end

    compute_observables_statistics!(model)
    return model
end

"""
    heat_bath!(model::Potts)

Implementation of the heat bath algorithm for the Potts model.

Arguments
* `model`: Potts type
"""
function heat_bath!(model::Potts)
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)

    for t = 1:model.params.n_sweeps
        for j = 1:model.params.L, i = 1:model.params.L
            counts = compute_counts(model, (i, j))
            old_spin = model.σ[i, j]
            ps = exp.(model.params.β .* counts)
            ps ./= sum(ps)
            R = rand(model.rng)
            P = 0.0

            for (new_spin, p) in enumerate(ps)
                P += p
                R < P && (model.σ[i, j] = new_spin; break)
            end
        end

        t > t₀ && update_observables!(model)
    end

    compute_observables_statistics!(model)
    return model
end

"""
    heat_bath!(model::XY)

Implementation of the heat bath algorithm for the XY model.

Arguments
* `model`: XY type
"""
function heat_bath!(model::XY)
    t₀ = floor(Int, model.params.cutoff * model.params.n_sweeps)
    eᵝ = exp(-model.params.β)

    for t = 1:model.params.n_sweeps
        for j = 1:model.params.L, i = 1:model.params.L
            old_spin = model.σ[i, j]
            new_spin = random_XYVector(model.rng)
            ΔE = dot(old_spin .- new_spin, sum_of_neighbors(model, (i, j)))

            (rand(model.rng) < eᵝ^ΔE / (eᵝ^ΔE + 1)) && (model.σ[i, j] = new_spin)
        end

        t > t₀ && update_observables!(model)
    end

    compute_observables_statistics!(model)
    return model
end


text:
#
# Optimizations:
#  - replace leftmap, rightmap with sorted vectors
#

###########################
# SiteOp                  #
###########################

struct SiteOp
  name::String
  site::Int
end

name(s::SiteOp) = s.name
site(s::SiteOp) = s.site

Base.show(io::IO,s::SiteOp) = print(io,"\"$(name(s))\"($(site(s)))")

Base.:(==)(s1::SiteOp,s2::SiteOp) = (s1.site==s2.site && s1.name==s2.name)

function Base.isless(s1::SiteOp,s2::SiteOp)::Bool
  if site(s1) < site(s2)
    return true
  end
  return name(s1) < name(s2)
end

###########################
# OpTerm                  #
###########################

const OpTerm = Vector{SiteOp}
mult(t1::OpTerm,t2::OpTerm) = isempty(t2) ? t1 : vcat(t1,t2)

function isfermionic(t::OpTerm,sites)::Bool
  p = +1
  for op in t
    if has_fermion_string(sites[site(op)],name(op))
      p *= -1
    end
  end
  return (p == -1)
end

###########################
# MPOTerm                 #
###########################

mutable struct MPOTerm
  coef::ComplexF64
  ops::OpTerm
end
coef(op::MPOTerm) = op.coef
ops(op::MPOTerm) = op.ops

function Base.:(==)(t1::MPOTerm,
                    t2::MPOTerm)
  return coef(t1) ≈ coef(t2) && ops(t1) == ops(t2)
end

function Base.isless(t1::MPOTerm, t2::MPOTerm)::Bool
  if !isapprox(coef(t1),coef(t2))
    ct1 = coef(t1)
    ct2 = coef(t2)
    #"lexicographic" ordering on  complex numbers
    return real(ct1) < real(ct2) ||
           (real(ct1) == real(ct2) && imag(ct1) < imag(ct2))
  end
  return ops(t1) < ops(t2)
end

function MPOTerm(c::Number,
                 op1::String,
                 i1::Int)
  return MPOTerm(convert(ComplexF64,c),[SiteOp(op1,i1)])
end

function MPOTerm(c::Number,
                 op1::String,i1::Int,
                 op2::String,i2::Int)
  return MPOTerm(convert(ComplexF64,c),[SiteOp(op1,i1),SiteOp(op2,i2)])
end

function MPOTerm(c::Number,
                 op1::String,i1::Int,
                 op2::String,i2::Int,
                 ops...)
  vop = OpTerm(undef,2+div(length(ops),2))
  vop[1] = SiteOp(op1,i1)
  vop[2] = SiteOp(op2,i2)
  for n = 1:div(length(ops),2)
    vop[2+n] = SiteOp(ops[2*n-1],ops[2*n])
  end
  return MPOTerm(convert(ComplexF64,c),vop)
end

#function MPOTerm(c::Number,
#                 ops::OpTerm)
#  return MPOTerm(convert(ComplexF64,c),ops)
#end

function Base.show(io::IO,
                   op::MPOTerm)
  c = coef(op)
  if iszero(imag(c))
    print(io,"$(real(c)) ")
  elseif iszero(real(c))
    print(io,"$(imag(c))im ")
  else
    print(io,"($c) ")
  end
  for o in ops(op)
    print(io,"\"$(name(o))\"($(site(o))) ")
  end
end

############################
## AutoMPO                 #
############################

"""
An AutoMPO stores a collection of
operator terms, to be later summed
together into an MPO by calling
the function `MPO` on the AutoMPO object.
Each term is a product of local operators
specified by names such as "Sz" or "N",
times an optional coefficient which
can be real or complex.

Which local operator names are available
is determined by the function `op`
associated with the TagType defined by
special Index tags, such as "S=1/2","S=1",
"Fermion", and "Electron".
"""
struct AutoMPO
  data::Vector{MPOTerm}
  AutoMPO(terms::Vector{MPOTerm}) = new(terms)
end

"""
    AutoMPO()

Construct an empty AutoMPO
"""
AutoMPO() = AutoMPO(Vector{MPOTerm}())

data(ampo::AutoMPO) = ampo.data

Base.:(==)(ampo1::AutoMPO,
           ampo2::AutoMPO) = data(ampo1) == data(ampo2)

Base.copy(ampo::AutoMPO) = AutoMPO(copy(data(ampo)))

Base.size(ampo::AutoMPO) = size(data(ampo))

"""
    add!(ampo::AutoMPO,
         op1::String, i1::Int)

    add!(ampo::AutoMPO,
         coef::Number,
         op1::String, i1::Int)

    add!(ampo::AutoMPO,
         op1::String, i1::Int,
         op2::String, i2::Int,
         ops...)

    add!(ampo::AutoMPO,
         coef::Number,
         op1::String, i1::Int,
         op2::String, i2::Int,
         ops...)

    +(ampo:AutoMPO, term::Tuple)

Add a single- or multi-site operator
term to the AutoMPO `ampo`. Each operator
is specified by a name (String) and a
site number (Int). The second version
accepts a real or complex coefficient.

The `+` operator version of this function
accepts a tuple with entries either
(String,Int,String,Int,...) or
(Number,String,Int,String,Int,...)
where these tuple values are the same
as valid inputs to the `add!` function.
For inputting a very large number of
terms (tuples) to an AutoMPO, consider
using the broadcasted operator `.+=`
which avoids reallocating the AutoMPO
after each addition.

# Examples
```julia
ampo = AutoMPO()

add!(ampo,"Sz",2,"Sz",3)

ampo += ("Sz",3,"Sz",4)

ampo += (0.5,"S+",4,"S-",5)

ampo .+= (0.5,"S+",5,"S-",6)
```
"""
function add!(ampo::AutoMPO,
              coef::Number,
              op::String, i::Int)
  push!(data(ampo),MPOTerm(coef,op,i))
  return
end

add!(ampo::AutoMPO,op::String, i::Int) = add!(ampo,1.0,op,i)

function add!(ampo::AutoMPO,
              coef::Number,
              op1::String, i1::Int,
              op2::String, i2::Int)
  push!(data(ampo),MPOTerm(coef,op1,i1,op2,i2))
  return
end

add!(ampo::AutoMPO,
     op1::String, i1::Int,
     op2::String, i2::Int) = add!(ampo,1.0,op1,i1,op2,i2)

function add!(ampo::AutoMPO,
              coef::Number,
              op1::String, i1::Int,
              op2::String, i2::Int,
              ops...)
  push!(data(ampo), MPOTerm(coef, op1, i1, op2, i2, ops...))
  return ampo
end

add!(ampo::AutoMPO,
     op1::String, i1::Int,
     op2::String, i2::Int,
     ops...) = add!(ampo, 1.0, op1, i1, op2, i2, ops...)

"""
    subtract!(ampo::AutoMPO,
              op1::String, i1::Int,
              op2::String, i2::Int,
              ops...)

    subtract!(ampo::AutoMPO,
              coef::Number,
              op1::String, i1::Int,
              op2::String, i2::Int,
              ops...)

Subtract a multi-site operator term
from the AutoMPO `ampo`. Each operator
is specified by a name (String) and a
site number (Int). The second version
accepts a real or complex coefficient.
"""
subtract!(ampo::AutoMPO,
          op1::String, i1::Int,
          op2::String, i2::Int,
          ops...) = add!(ampo, -1.0, op1, i1, op2, i2, ops...)

function subtract!(ampo::AutoMPO,
                   coef::Number,
                   op1::String, i1::Int,
                   op2::String, i2::Int,
                   ops...)
  push!(data(ampo), -MPOTerm(coef, op1, i1, op2, i2, ops...))
  return ampo
end

Base.:-(t::MPOTerm) = MPOTerm(-coef(t), ops(t))

function Base.:+(ampo::AutoMPO,
                 term::Tuple)
  ampo_plus_term = copy(ampo)
  add!(ampo_plus_term, term...)
  return ampo_plus_term
end

function Base.:-(ampo::AutoMPO,
                 term::Tuple)
  ampo_plus_term = copy(ampo)
  subtract!(ampo_plus_term, term...)
  return ampo_plus_term
end

#
# ampo .+= ("Sz",1) syntax using broadcasting
#

struct AutoMPOStyle <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{<:AutoMPO}) = AutoMPOStyle()

struct AutoMPOAddTermStyle <: Broadcast.BroadcastStyle end

Base.broadcastable(ampo::AutoMPO) = ampo

Base.BroadcastStyle(::AutoMPOStyle,
                    ::Broadcast.Style{Tuple}) = AutoMPOAddTermStyle()

Broadcast.instantiate(bc::Broadcast.Broadcasted{AutoMPOAddTermStyle}) = bc

function Base.copyto!(ampo,
                      bc::Broadcast.Broadcasted{AutoMPOAddTermStyle,
                                                <:Any,
                                                typeof(+)})
  add!(ampo, bc.args[2]...)
  return ampo
end

#
# ampo .-= ("Sz",1) syntax using broadcasting
#

function Base.copyto!(ampo,
                      bc::Broadcast.Broadcasted{AutoMPOAddTermStyle,
                                                <:Any,
                                                typeof(-)})
  subtract!(ampo, bc.args[2]...)
  return ampo
end

function Base.show(io::IO,
                   ampo::AutoMPO)
  println(io,"AutoMPO:")
  for term in data(ampo)
    println(io,"  $term")
  end
end

##################################
# MatElem (simple sparse matrix) #
##################################

struct MatElem{T}
  row::Int
  col::Int
  val::T
end

#function Base.show(io::IO,m::MatElem)
#  print(io,"($(m.row),$(m.col),$(m.val))")
#end

function toMatrix(els::Vector{MatElem{T}})::Matrix{T} where {T}
  nr = 0
  nc = 0
  for el in els
    nr = max(nr,el.row)
    nc = max(nc,el.col)
  end
  M = zeros(T,nr,nc)
  for el in els
    M[el.row,el.col] = el.val
  end
  return M
end

function Base.:(==)(m1::MatElem{T},m2::MatElem{T})::Bool where {T}
  return (m1.row==m2.row && m1.col==m2.col && m1.val==m2.val)
end

function Base.isless(m1::MatElem{T},m2::MatElem{T})::Bool where {T}
  if m1.row != m2.row
    return m1.row < m2.row
  elseif m1.col != m2.col
    return m1.col < m2.col
  end
  return m1.val < m2.val
end

struct QNMatElem{T}
  rowqn::QN
  colqn::QN
  row::Int
  col::Int
  val::T
end

function Base.:(==)(m1::QNMatElem{T},m2::QNMatElem{T})::Bool where {T}
  return (m1.row==m2.row && m1.col==m2.col && m1.val==m2.val && m1.rowqn==m2.rowqn && m1.colqn==m2.colqn)
end

function Base.isless(m1::QNMatElem{T},m2::QNMatElem{T})::Bool where {T}
  if m1.rowqn != m2.rowqn
    return m1.rowqn < m2.rowqn
  elseif m1.colqn != m2.colqn
    return m1.colqn < m2.colqn
  elseif m1.row != m2.row
    return m1.row < m2.row
  elseif m1.col != m2.col
    return m1.col < m2.col
  end
  return m1.val < m2.val
end

function posInLink!(linkmap::Vector{OpTerm},
                    op::OpTerm)::Int
  isempty(op) && return -1
  for n=1:length(linkmap)
    (linkmap[n]==op) && return n
  end
  push!(linkmap,op)
  return length(linkmap)
end

function determineValType(terms::Vector{MPOTerm})
  for t in terms
    (!isreal(coef(t))) && return ComplexF64
  end
  return Float64
end

function computeSiteProd(sites,
                         ops::OpTerm)::ITensor
  i = ops[1].site
  T = op(sites[i],ops[1].name)
  for j=2:length(ops)
    (ops[j].site != i) && error("Mismatch of site number in computeSiteProd")
    opj = op(sites[i],ops[j].name)
    T = product(T, opj)
  end
  return T
end

function remove_dups!(v::Vector{T}) where {T}
  N = length(v)
  (N==0) && return
  sort!(v)
  n = 1
  u = 2
  while u <= N
    while u < N && v[u]==v[n]
      u += 1
    end
    if v[u] != v[n]
      v[n+1] = v[u]
      n += 1
    end
    u += 1
  end
  resize!(v,n)
  return
end #remove_dups!

function svdMPO(ampo::AutoMPO,
                sites;
                kwargs...)::MPO

  mindim::Int = get(kwargs,:mindim,1)
  maxdim::Int = get(kwargs,:maxdim,10000)
  cutoff::Float64 = get(kwargs,:cutoff,1E-13)

  N = length(sites)

  ValType = determineValType(data(ampo))

  Vs = [Matrix{ValType}(undef,1,1) for n=1:N]
  tempMPO = [MatElem{MPOTerm}[] for n=1:N]

  crosses_bond(t::MPOTerm,n::Int) = (ops(t)[1].site <= n <= ops(t)[end].site)

  rightmap = OpTerm[]
  next_rightmap = OpTerm[]

  for n=1:N

    leftbond_coefs = MatElem{ValType}[]

    leftmap = OpTerm[]
    for term in data(ampo)
      crosses_bond(term,n) || continue

      left::OpTerm   = filter(t->(t.site < n),ops(term))
      onsite::OpTerm = filter(t->(t.site == n),ops(term))
      right::OpTerm  = filter(t->(t.site > n),ops(term))

      bond_row = -1
      bond_col = -1
      if !isempty(left)
        bond_row = posInLink!(leftmap,left)
        bond_col = posInLink!(rightmap,mult(onsite,right))
        bond_coef = convert(ValType,coef(term))
        push!(leftbond_coefs,MatElem(bond_row,bond_col,bond_coef))
      end

      A_row = bond_col
      A_col = posInLink!(next_rightmap,right)
      site_coef = 1.0+0.0im
      if A_row == -1
        site_coef = coef(term)
      end
      if isempty(onsite)
        if isfermionic(right,sites)
          #println("Putting F on site $n")
          push!(onsite,SiteOp("F",n))
        else
          #println("Putting Id on site $n")
          push!(onsite,SiteOp("Id",n))
        end
      end
      el = MatElem(A_row,A_col,MPOTerm(site_coef,onsite))
      push!(tempMPO[n],el)
    end
    rightmap = next_rightmap
    next_rightmap = OpTerm[]

    remove_dups!(tempMPO[n])

    if n > 1 && !isempty(leftbond_coefs)
      M = toMatrix(leftbond_coefs)
      U,S,V = svd(M)
      P = S.^2
      truncate!(P;maxdim=maxdim,cutoff=cutoff,mindim=mindim)
      tdim = length(P)
      nc = size(M,2)
      Vs[n-1] = Matrix{ValType}(V[1:nc,1:tdim])
    end

  end

  llinks = Vector{Index{Int}}(undef, N+1)
  llinks[1] = Index(2, "Link,l=0")

  H = MPO(sites)

  # Constants which define MPO start/end scheme
  rowShift = 2
  colShift = 2
  startState = 2
  endState = 1

  for n=1:N
    VL = Matrix{ValType}(undef,1,1)
    if n > 1
      VL = Vs[n-1]
    end
    VR = Vs[n]
    tdim = size(VR,2)

    llinks[n+1] = Index(2+tdim,"Link,l=$n")

    finalMPO = Dict{OpTerm,Matrix{ValType}}()

    ll = llinks[n]
    rl = llinks[n+1]

    idTerm = [SiteOp("Id",n)]
    finalMPO[idTerm] = zeros(ValType,dim(ll),dim(rl))
    idM = finalMPO[idTerm]
    idM[1,1] = 1.0
    idM[2,2] = 1.0

    defaultMat() = zeros(ValType,dim(ll),dim(rl))

    for el in tempMPO[n]
      A_row = el.row
      A_col = el.col
      t = el.val
      (abs(coef(t)) > eps()) || continue

      M = get!(finalMPO,ops(t),defaultMat())

      ct = convert(ValType,coef(t))
      if A_row==-1 && A_col==-1 #onsite term
        M[startState,endState] += ct
      elseif A_row==-1 #term starting on site n
        for c=1:size(VR,2)
          z = ct*VR[A_col,c]
          M[startState,colShift+c] += z
        end
      elseif A_col==-1 #term ending on site n
        for r=1:size(VL,2)
          z = ct*conj(VL[A_row,r])
          M[rowShift+r,endState] += z
        end
      else
        for r=1:size(VL,2),c=1:size(VR,2)
          z = ct*conj(VL[A_row,r])*VR[A_col,c]
          M[rowShift+r,colShift+c] += z
        end
      end
    end

    s = sites[n]
    #H[n] = emptyITensor(dag(s),s',ll,rl)
    H[n] = emptyITensor(Any)
    for (op,M) in finalMPO
      T = itensor(M,ll,rl)
      H[n] += T*computeSiteProd(sites,op)
    end

  end

  L = emptyITensor(llinks[1])
  L[startState] = 1.0

  R = emptyITensor(llinks[N+1])
  R[endState] = 1.0

  H[1] *= L
  H[N] *= R

  return H
end #svdMPO

function qn_svdMPO(ampo::AutoMPO,
                   sites;
                   kwargs...)::MPO

  mindim::Int = get(kwargs,:mindim,1)
  maxdim::Int = get(kwargs,:maxdim,10000)
  cutoff::Float64 = get(kwargs,:cutoff,1E-13)

  N = length(sites)

  ValType = determineValType(data(ampo))

  Vs = [Dict{QN,Matrix{ValType}}() for n=1:N+1]
  tempMPO = [QNMatElem{MPOTerm}[] for n=1:N]

  crosses_bond(t::MPOTerm,n::Int) = (ops(t)[1].site <= n <= ops(t)[end].site)

  rightmap = Dict{QN,Vector{OpTerm}}()
  next_rightmap = Dict{QN,Vector{OpTerm}}()

  for n=1:N

    leftbond_coefs = Dict{QN,Vector{MatElem{ValType}}}()

    leftmap = Dict{QN,Vector{OpTerm}}()
    for term in data(ampo)
      crosses_bond(term,n) || continue

      left::OpTerm   = filter(t->(t.site < n),ops(term))
      onsite::OpTerm = filter(t->(t.site == n),ops(term))
      right::OpTerm  = filter(t->(t.site > n),ops(term))

      function calcQN(term::OpTerm)
        q = QN()
        for st in term
          op_tensor = op(sites[site(st)],name(st))
          q -= flux(op_tensor)
        end
        return q
      end
      lqn = calcQN(left)
      sqn = calcQN(onsite)

      q_leftmap = get!(leftmap,lqn,OpTerm[])
      q_rightmap = get!(rightmap,lqn,OpTerm[])

      bond_row = -1
      bond_col = -1
      if !isempty(left)
        bond_row = posInLink!(q_leftmap,left)
        bond_col = posInLink!(q_rightmap,mult(onsite,right))
        bond_coef = convert(ValType,coef(term))
        q_leftbond_coefs = get!(leftbond_coefs,lqn,MatElem{ValType}[])
        push!(q_leftbond_coefs,MatElem(bond_row,bond_col,bond_coef))
      end

      rqn = sqn+lqn
      q_next_rightmap = get!(next_rightmap,rqn,OpTerm[])
      A_row = bond_col
      A_col = posInLink!(q_next_rightmap,right)
      site_coef = 1.0+0.0im
      if A_row == -1
        site_coef = coef(term)
      end
      if isempty(onsite)
        if isfermionic(right,sites)
          push!(onsite,SiteOp("F",n))
        else
          push!(onsite,SiteOp("Id",n))
        end
      end
      el = QNMatElem(lqn,rqn,A_row,A_col,MPOTerm(site_coef,onsite))
      push!(tempMPO[n],el)
    end
    rightmap = next_rightmap
    next_rightmap = Dict{QN,Vector{OpTerm}}()

    remove_dups!(tempMPO[n])

    if n > 1 && !isempty(leftbond_coefs)
      for (q,mat) in leftbond_coefs
        M = toMatrix(mat)
        U,S,V = svd(M)
        P = S.^2
        truncate!(P;maxdim=maxdim,cutoff=cutoff,mindim=mindim)
        tdim = length(P)
        nc = size(M,2)
        Vs[n][q] = Matrix{ValType}(V[1:nc,1:tdim])
      end
    end
  end

  #
  # Make MPO link indices
  #
  d0 = 2
  llinks = Vector{QNIndex}(undef,N+1)
  llinks[1] = Index(QN()=>d0;tags="Link,l=0")
  for n=1:N
    qi = Vector{Pair{QN,Int}}()
    if !haskey(Vs[n+1],QN())
      # Make sure QN=zero is first in list of sectors
      push!(qi,QN()=>d0)
    end
    for (q,Vq) in Vs[n+1]
      cols = size(Vq,2)
      if q==QN()
        # Make sure QN=zero is first in list of sectors
        insert!(qi,1,q=>d0+cols)
      else
        push!(qi,q=>cols)
      end
    end
    llinks[n+1] = Index(qi...;tags="Link,l=$n")
  end

  H = MPO(N)

  # Constants which define MPO start/end scheme
  startState = 2
  endState = 1

  for n=1:N
    finalMPO = Dict{Tuple{QN,OpTerm},Matrix{ValType}}()

    ll = llinks[n]
    rl = llinks[n+1]

    function defaultMat(ll,rl,lqn,rqn)
      ldim = qnblockdim(ll,lqn)
      rdim = qnblockdim(rl,rqn)
      return zeros(ValType,ldim,rdim)
    end

    idTerm = [SiteOp("Id",n)]
    finalMPO[(QN(),idTerm)] = defaultMat(ll,rl,QN(),QN())
    idM = finalMPO[(QN(),idTerm)]
    idM[1,1] = 1.0
    idM[2,2] = 1.0

    for el in tempMPO[n]
      t = el.val
      (abs(coef(t)) > eps()) || continue
      A_row = el.row
      A_col = el.col

      M = get!(finalMPO,(el.rowqn,ops(t)),defaultMat(ll,rl,el.rowqn,el.colqn))

      # rowShift and colShift account for
      # special entries in the zero-QN sector
      # of the MPO
      rowShift = (el.rowqn == QN()) ? 2 : 0
      colShift = (el.colqn == QN()) ? 2 : 0

      ct = convert(ValType,coef(t))
      if A_row==-1 && A_col==-1 #onsite term
        M[startState,endState] += ct
      elseif A_row==-1 #term starting on site n
        VR = Vs[n+1][el.colqn]
        for c=1:size(VR,2)
          z = ct*VR[A_col,c]
          M[startState,colShift+c] += z
        end
      elseif A_col==-1 #term ending on site n
        VL = Vs[n][el.rowqn]
        for r=1:size(VL,2)
          z = ct*conj(VL[A_row,r])
          M[rowShift+r,endState] += z
        end
      else
        VL = Vs[n][el.rowqn]
        VR = Vs[n+1][el.colqn]
        for r=1:size(VL,2),c=1:size(VR,2)
          z = ct*conj(VL[A_row,r])*VR[A_col,c]
          M[rowShift+r,colShift+c] += z
        end
      end
    end

    s = sites[n]
    #H[n] = emptyITensor(dag(s),s',dag(ll),rl)
    H[n] = emptyITensor(Any)
    for (q_op,M) in finalMPO
      op_prod = q_op[2]
      Op = computeSiteProd(sites,op_prod)

      rq = q_op[1]
      sq = flux(Op)
      cq = rq-sq

      rn = qnblocknum(ll,rq)
      cn = qnblocknum(rl,cq)

      #TODO: wrap following 3 lines into a function
      block = (rn,cn)
      T = BlockSparseTensor(ValType,[block],IndexSet(dag(ll),rl))
      blockview(T,block) .= M

      IT = itensor(T)
      H[n] += IT * Op
    end

  end

  L = emptyITensor(llinks[1])
  L[startState] = 1.0

  R = emptyITensor(dag(llinks[N+1]))
  R[endState] = 1.0

  H[1] *= L
  H[N] *= R

  return H
end #qn_svdMPO

function sorteachterm(ampo::AutoMPO, sites)
  ampo = copy(ampo)
  isless_site(o1::SiteOp, o2::SiteOp) = site(o1) < site(o2)
  for t in data(ampo)
    Nt = length(t.ops)

    # Sort operators in t by site order,
    # and keep the permutation used, perm, for analysis below
    perm = Vector{Int}(undef,Nt)
    sortperm!(perm,t.ops, alg=InsertionSort, lt=isless_site)

    t.ops = t.ops[perm]

    # Identify fermionic operators,
    # zeroing perm for bosonic operators,
    # and inserting string "F" operators
    parity = +1
    for n=Nt:-1:1
      fermionic = has_fermion_string(sites[site(t.ops[n])],name(t.ops[n]))
      if parity == -1
        # Put Jordan-Wigner string emanating
        # from fermionic operators to the right
        # (Remaining F operators will be put in by svdMPO)
        t.ops[n] = SiteOp("$(name(t.ops[n]))*F",site(t.ops[n]))
      end
      if fermionic
        parity = -parity
      else
        # Ignore bosonic operators in perm
        # by zeroing corresponding entries
        perm[n] = 0
      end
    end
    if parity == -1
      error("Parity-odd fermionic terms not yet supported by AutoMPO")
    end
    # Keep only fermionic op positions (non-zero entries)
    filter!(!iszero,perm)
    # Account for anti-commuting, fermionic operators
    # during above sort; put resulting sign into coef

    t.coef *= parity_sign(perm)
  end
  return ampo
end

"""
    MPO(ampo::AutoMPO,sites::Vector{<:Index};kwargs...)

Convert an AutoMPO object `ampo` to an
MPO, with indices given by `sites`. The
resulting MPO will have the indices
`sites[1], sites[1]', sites[2], sites[2]'`
etc. The conversion is done by an algorithm
that compresses the MPO resulting from adding
the AutoMPO terms together, often achieving
the minimum possible bond dimension.

# Examples
```julia
ampo = AutoMPO()
ampo += ("Sz",1,"Sz",2)
ampo += ("Sz",2,"Sz",3)
ampo += ("Sz",3,"Sz",4)

sites = siteinds("S=1/2",4)
H = MPO(ampo,sites)
```
"""
function MPO(ampo::AutoMPO,
             sites::Vector{<:Index};
             kwargs...)::MPO
  ampo = sorteachterm(ampo,sites)
  #for t in data(ampo)
  #  @show t
  #end
  if hasqns(sites[1])
    return qn_svdMPO(ampo,sites;kwargs...)
  end
  return svdMPO(ampo,sites;kwargs...)
end

@deprecate toMPO(args...; kwargs...) MPO(args...; kwargs...)


text:
using BinDeps
using Compat

@BinDeps.setup

sqlite3_lib = library_dependency("sqlite3_lib", aliases = ["sqlite3","sqlite3-64","libsqlite3","libsqlite3-0"])

if is_windows()
    using WinRPM
    provides(WinRPM.RPM, "libsqlite3-0",sqlite3_lib, os = :Windows)
elseif is_apple()
    using Homebrew
    provides(Homebrew.HB, "sqlite", sqlite3_lib, os = :Darwin)
elseif is_linux()
    provides(Yum, Dict("sqlite3-devel" => sqlite3_lib))
    provides(AptGet, Dict("sqlite3" => sqlite3_lib))
end

@BinDeps.install Dict(:sqlite3_lib => :sqlite3_lib)


text:
# Test cases for PTI RAW file parser
using Memento, Memento.Test

@testset "Test PTI '.raw' file parser" begin
    @testset "Check PTI exception handling" begin
        setlevel!(getlogger(PowerModels), "warn")

        @test_nowarn PowerModels.parse_pti("../test/data/parser_test_a.raw")
        @test_throws(getlogger(PowerModels),
                    MethodError,
                    @test_warn(getlogger(PowerModels),
                               "This feature is incomplete, and will currently only return RAW data from a PTI file",
                               PowerModels.parse_file("../test/data/frankenstein_00.raw")))
        @test_throws(getlogger(PowerModels),
                     ErrorException,
                     PowerModels.parse_pti("../test/data/parser_test_b.raw"))
        @test_warn(getlogger(PowerModels),
                   "Version 32 of PTI format is unsupported, parser may not function correctly.",
                   PowerModels.parse_pti("../test/data/parser_test_c.raw"))
        @test_throws(getlogger(PowerModels), ErrorException, PowerModels.parse_pti("../test/data/parser_test_d.raw"))
        @test_warn(getlogger(PowerModels), "GNE DEVICE parsing is not supported.", PowerModels.parse_pti("../test/data/parser_test_h.raw"))

        setlevel!(getlogger(PowerModels), "error")
    end

    @testset "4-bus frankenstein file" begin
        data_dict = PowerModels.parse_pti("../test/data/frankenstein_00.raw")
        @test isa(data_dict, Dict)

        @test length(data_dict["CASE IDENTIFICATION"]) == 1
        @test length(data_dict["CASE IDENTIFICATION"][1]) == 8

        @test length(data_dict["BUS"]) == 4
        for item in data_dict["BUS"]
            @test length(item) == 13
        end

        @test length(data_dict["LOAD"]) == 2
        for item in data_dict["LOAD"]
            @test length(item) == 13
        end

        @test length(data_dict["FIXED SHUNT"]) == 1
        @test length(data_dict["FIXED SHUNT"][1]) == 5

        @test length(data_dict["GENERATOR"]) == 3
        for item in data_dict["GENERATOR"]
            @test length(item) == 28
        end

        @test length(data_dict["TRANSFORMER"]) == 1
        @test length(data_dict["TRANSFORMER"][1]) == 43

        @test length(data_dict["AREA INTERCHANGE"]) == 1
        @test length(data_dict["AREA INTERCHANGE"][1]) == 5

        @test length(data_dict["ZONE"]) == 1
        @test length(data_dict["ZONE"][1]) == 2

        @test length(data_dict["OWNER"]) == 1
        @test length(data_dict["OWNER"][1]) == 2
    end

    @testset "20-bus frankenstein file" begin
        data_dict = PowerModels.parse_pti("../test/data/frankenstein_20.raw")
        @test isa(data_dict, Dict)

        @test length(data_dict["BRANCH"]) == 2
        for item in data_dict["BRANCH"]
            @test length(item) == 24
        end

        @test length(data_dict["TRANSFORMER"][2]) == 83

        @test length(data_dict["SWITCHED SHUNT"]) == 2
        @test length(data_dict["SWITCHED SHUNT"][1]) == 16
        @test length(data_dict["SWITCHED SHUNT"][2]) == 12
    end

    @testset "70-bus frankenstein file" begin
        data_dict = PowerModels.parse_pti("../test/data/frankenstein_70.raw")
        @test isa(data_dict, Dict)

        @test length(data_dict["TWO-TERMINAL DC"]) == 1
        @test length(data_dict["TWO-TERMINAL DC"][1]) == 46

        @test length(data_dict["VOLTAGE SOURCE CONVERTER"]) == 1
        @test length(data_dict["VOLTAGE SOURCE CONVERTER"][1]) == 12
        @test length(data_dict["VOLTAGE SOURCE CONVERTER"][1]["CONVERTER BUSES"]) == 2
        for item in data_dict["VOLTAGE SOURCE CONVERTER"][1]["CONVERTER BUSES"]
            @test length(item) == 15
        end

        @test length(data_dict["IMPEDANCE CORRECTION"]) == 2
        @test length(data_dict["IMPEDANCE CORRECTION"][1]) == 23
        @test length(data_dict["IMPEDANCE CORRECTION"][2]) == 15

        @test length(data_dict["MULTI-TERMINAL DC"]) == 1
        @test length(data_dict["MULTI-TERMINAL DC"][1]) == 11

        @test length(data_dict["MULTI-TERMINAL DC"][1]["CONV"]) == 3
        for item in data_dict["MULTI-TERMINAL DC"][1]["CONV"]
            @test length(item) == 16
        end

        @test length(data_dict["MULTI-TERMINAL DC"][1]["DCBS"]) == 3
        for item in data_dict["MULTI-TERMINAL DC"][1]["DCBS"]
            @test length(item) == 8
        end

        @test length(data_dict["MULTI-TERMINAL DC"][1]["DCLN"]) == 3
        for item in data_dict["MULTI-TERMINAL DC"][1]["DCLN"]
            @test length(item) == 6
        end

        @test length(data_dict["FACTS CONTROL DEVICE"]) == 2
        for item in data_dict["FACTS CONTROL DEVICE"]
            @test length(item) == 19
        end

        data_dict = PowerModels.parse_pti("../test/data/parser_test_e.raw")
        @test length(data_dict["MULTI-TERMINAL DC"][1]) == 10

        data_dict = PowerModels.parse_pti("../test/data/parser_test_f.raw")
        @test length(data_dict["MULTI-TERMINAL DC"][1]) == 9

        data_dict = PowerModels.parse_pti("../test/data/parser_test_g.raw")
        @test length(data_dict["MULTI-TERMINAL DC"][1]) == 8

    end

    @testset "0-bus case file" begin
        @test_throws(getlogger(PowerModels), ArgumentError, PowerModels.parse_pti("../test/data/case0.raw"))
    end

    @testset "73-bus case file" begin
        data_dict = PowerModels.parse_pti("../test/data/case73.raw")
        @test isa(data_dict, Dict)

        @test length(data_dict["BUS"]) == 73
        for item in data_dict["BUS"]
            @test length(item) == 13
        end

        @test length(data_dict["LOAD"]) == 51
        for item in data_dict["LOAD"]
            @test length(item) == 13
        end

        @test length(data_dict["GENERATOR"]) == 99
        for item in data_dict["GENERATOR"]
            @test length(item) == 28
        end

        @test length(data_dict["BRANCH"]) == 105
        for item in data_dict["BRANCH"]
            @test length(item) == 24
        end

        @test length(data_dict["TRANSFORMER"]) == 15
        for item in data_dict["TRANSFORMER"]
            @test length(item) == 43
        end

        @test length(data_dict["AREA INTERCHANGE"]) == 3
        @test length(data_dict["ZONE"]) == 3
        @test length(data_dict["OWNER"]) == 1

        @test length(data_dict["SWITCHED SHUNT"]) == 3
        for item in data_dict["SWITCHED SHUNT"]
            @test length(item) == 12
        end
    end
end