text:
function iterate!(job::BasicMCJob, ::Type{SMMALA}, ::Type{Univariate})
  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    job.sstate.tune.proposed += 1
  end

  job.sstate.μ = job.pstate.value+0.5*job.sstate.tune.step*job.sstate.oldfirstterm
  job.sstate.pstate.value = job.sstate.μ+job.sstate.sqrttunestep*job.sstate.cholinvtensor*randn()

  job.parameter.uptotensorlogtarget!(job.sstate.pstate)

  job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget

  job.sstate.ratio += (
    0.5*(
      log(job.sstate.tune.step*job.sstate.oldinvtensor)
      +abs2(job.sstate.pstate.value-job.sstate.μ)*job.pstate.tensorlogtarget/job.sstate.tune.step
    )
  )

  job.sstate.newinvtensor = inv(job.sstate.pstate.tensorlogtarget)

  job.sstate.newfirstterm = job.sstate.newinvtensor*job.sstate.pstate.gradlogtarget

  job.sstate.μ = job.sstate.pstate.value+0.5*job.sstate.tune.step*job.sstate.newfirstterm

  job.sstate.ratio -= (
    0.5*(
      log(job.sstate.tune.step*job.sstate.newinvtensor)
      +abs2(job.pstate.value-job.sstate.μ)*job.sstate.pstate.tensorlogtarget/job.sstate.tune.step
    )
  )

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = job.sstate.pstate.value

    job.pstate.gradlogtarget = job.sstate.pstate.gradlogtarget
    if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
      job.pstate.gradloglikelihood = job.sstate.pstate.gradloglikelihood
    end
    if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
      job.pstate.gradlogprior = job.sstate.pstate.gradlogprior
    end

    job.pstate.tensorlogtarget = job.sstate.pstate.tensorlogtarget
    if in(:tensorloglikelihood, job.outopts[:monitor]) && job.parameter.tensorloglikelihood! != nothing
      job.pstate.tensorloglikelihood = job.sstate.pstate.tensorloglikelihood
    end
    if in(:tensorlogprior, job.outopts[:monitor]) && job.parameter.tensorlogprior! != nothing
      job.pstate.tensorlogprior = job.sstate.pstate.tensorlogprior
    end

    job.sstate.oldinvtensor = job.sstate.newinvtensor
    job.sstate.cholinvtensor = chol(job.sstate.newinvtensor)
    job.sstate.oldfirstterm = job.sstate.newfirstterm

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end

    if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
      rate!(job.sstate.tune)

      if isa(job.tuner, AcceptanceRateMCTuner)
        tune!(job.sstate.tune, job.tuner)
        job.sstate.sqrttunestep = sqrt(job.sstate.tune.step)
      end

      if job.tuner.verbose
        println(
          "Burnin iteration ",
          job.fmt_iter(job.sstate.tune.totproposed),
          " of ",
          job.range.burnin,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate),
          " % acceptance rate"
        )
      end

      reset_burnin!(job.sstate.tune)
    end
  end
end

function iterate!(job::BasicMCJob, ::Type{SMMALA}, ::Type{Multivariate})
  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    job.sstate.tune.proposed += 1
  end

  job.sstate.μ[:] = job.pstate.value+0.5*job.sstate.tune.step*job.sstate.oldfirstterm
  job.sstate.pstate.value[:] = job.sstate.μ+job.sstate.sqrttunestep*job.sstate.cholinvtensor*randn(job.pstate.size)

  job.parameter.uptotensorlogtarget!(job.sstate.pstate)

  if job.sampler.transform != nothing
    job.sstate.pstate.tensorlogtarget[:, :] = job.sampler.transform(job.sstate.pstate.tensorlogtarget)
  end

  job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget

  job.sstate.ratio += (
    0.5*(
      logdet(job.sstate.tune.step*job.sstate.oldinvtensor)
      +dot(
        job.sstate.pstate.value-job.sstate.μ,
        job.pstate.tensorlogtarget*(job.sstate.pstate.value-job.sstate.μ)
      )/job.sstate.tune.step
    )
  )

  job.sstate.newinvtensor[:, :] = inv(job.sstate.pstate.tensorlogtarget)

  job.sstate.newfirstterm[:] = job.sstate.newinvtensor*job.sstate.pstate.gradlogtarget

  job.sstate.μ[:] = job.sstate.pstate.value+0.5*job.sstate.tune.step*job.sstate.newfirstterm

  job.sstate.ratio -= (
    0.5*(
      logdet(job.sstate.tune.step*job.sstate.newinvtensor)
      +dot(
        job.pstate.value-job.sstate.μ,
        job.sstate.pstate.tensorlogtarget*(job.pstate.value-job.sstate.μ)
      )/job.sstate.tune.step
    )
  )

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = copy(job.sstate.pstate.value)

    job.pstate.gradlogtarget = copy(job.sstate.pstate.gradlogtarget)
    if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
      job.pstate.gradloglikelihood = copy(job.sstate.pstate.gradloglikelihood)
    end
    if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
      job.pstate.gradlogprior = copy(job.sstate.pstate.gradlogprior)
    end

    job.pstate.tensorlogtarget = copy(job.sstate.pstate.tensorlogtarget)
    if in(:tensorloglikelihood, job.outopts[:monitor]) && job.parameter.tensorloglikelihood! != nothing
      job.pstate.tensorloglikelihood = copy(job.sstate.pstate.tensorloglikelihood)
    end
    if in(:tensorlogprior, job.outopts[:monitor]) && job.parameter.tensorlogprior! != nothing
      job.pstate.tensorlogprior = copy(job.sstate.pstate.tensorlogprior)
    end

    job.sstate.oldinvtensor = copy(job.sstate.newinvtensor)
    job.sstate.cholinvtensor[:, :] = ctranspose(chol(Hermitian(job.sstate.newinvtensor)))
    job.sstate.oldfirstterm = copy(job.sstate.newfirstterm)

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end

    if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
      rate!(job.sstate.tune)

      if isa(job.tuner, AcceptanceRateMCTuner)
        tune!(job.sstate.tune, job.tuner)
        job.sstate.sqrttunestep = sqrt(job.sstate.tune.step)
      end

      if job.tuner.verbose
        println(
          "Burnin iteration ",
          job.fmt_iter(job.sstate.tune.totproposed),
          " of ",
          job.range.burnin,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate),
          " % acceptance rate"
        )
      end

      reset_burnin!(job.sstate.tune)
    end
  end
end


text:
function start_load_failure!(cur_package::Module, cur_cabinet::FileCabinet; verbose::Bool=false)

  if verbose

    println("\ninvalid files:\n")

    for cur_info in cur_cabinet.file_infos

      println(cur_info.name * " - " * string(cur_info.undef))

    end

    println("-----")

  end

  cur_error = nothing

  for (cur_index, cur_info) in enumerate(cur_cabinet.file_infos)

    try
      load_invalid_file(cur_package, cur_info)
    catch tmp_error
      cur_error = tmp_error

      verbose || break
      ( cur_index == length(cur_cabinet.file_infos) ) && break

      showerror(STDERR, tmp_error)
      println("\n-----")
    end

  end

  rethrow(cur_error)

end

ext: jl


text:
using NBInclude

@testset "Jupyter notebooks" begin
    notebooks = ["Pitching Plate - K07",
                 "Point Source Demo",
                 "Translating Plate - 20°",
                 "Translating Plate - 60°",
                 "Vortex Sheet Roll-up",
                 "Moving bodies",
                 "Airfoil in uniform flow"]

    notebook_dir = joinpath(splitdir(@__FILE__())[1], "../binder/notebooks")

    @testset "$notebook" for notebook in notebooks
        @nbinclude(joinpath(notebook_dir, "$notebook.ipynb"); softscope=true)
        @test true
    end
end


text:
using Unicode

abstract type CallInfo end

# Call could be resolved to a singular MI
struct MICallInfo <: CallInfo
    mi::MethodInstance
    rt
    function MICallInfo(mi::MethodInstance, @nospecialize(rt))
        if isa(rt, LimitedAccuracy)
            return LimitedCallInfo(new(mi, ignorelimited(rt)))
        else
            return new(mi, rt)
        end
    end
end
get_mi(ci::MICallInfo) = ci.mi
get_rt(ci::CallInfo) = ci.rt

abstract type WrappedCallInfo <: CallInfo end

get_wrapped(ci::WrappedCallInfo) = ci.wrapped
ignorewrappers(ci::CallInfo) = ci
ignorewrappers(ci::WrappedCallInfo) = ignorewrappers(get_wrapped(ci))
get_mi(ci::WrappedCallInfo) = get_mi(ignorewrappers(ci))
get_rt(ci::WrappedCallInfo) = get_rt(ignorewrappers(ci))

# only appears when inspecting pre-optimization states
struct LimitedCallInfo <: WrappedCallInfo
    wrapped::CallInfo
end

# uncached callsite, we can't recurse into this call
struct UncachedCallInfo <: WrappedCallInfo
    wrapped::CallInfo
end

struct PureCallInfo <: CallInfo
    argtypes::Vector{Any}
    rt
    PureCallInfo(argtypes::Vector{Any}, @nospecialize(rt)) =
        new(argtypes, rt)
end
get_mi(::PureCallInfo) = nothing

# Failed
struct FailedCallInfo <: CallInfo
    sig
    rt
end
get_mi(ci::FailedCallInfo) = fail(ci)
get_rt(ci::FailedCallInfo) = fail(ci)
function fail(ci::FailedCallInfo)
    @error "MethodInstance extraction failed" ci.sig ci.rt
    return nothing
end

# Generated
struct GeneratedCallInfo <: CallInfo
    sig
    rt
end
function get_mi(genci::GeneratedCallInfo)
    @error "Can't extract MethodInstance from call to generated functions" genci.sig genci.rt
    return nothing
end

struct MultiCallInfo <: CallInfo
    sig
    rt
    callinfos::Vector{CallInfo}
end
# actual code-error
get_mi(ci::MultiCallInfo) = error("Can't extract MethodInstance from multiple call informations")

struct TaskCallInfo <: CallInfo
    ci::CallInfo
end
get_mi(tci::TaskCallInfo) = get_mi(tci.ci)
get_rt(tci::TaskCallInfo) = get_rt(tci.ci)

struct InvokeCallInfo <: CallInfo
    ci::MICallInfo
    InvokeCallInfo(mi::MethodInstance, @nospecialize(rt)) =
        new(MICallInfo(mi, rt))
end
get_mi(ci::InvokeCallInfo) = get_mi(ci.ci)
get_rt(ci::InvokeCallInfo) = get_rt(ci.ci)

# OpaqueClosure CallInfo
struct OCCallInfo <: CallInfo
    ci::MICallInfo
end
OCCallInfo(mi::MethodInstance, rt) = OCCallInfo(MICallInfo(mi, rt))
get_mi(tci::OCCallInfo) = get_mi(tci.ci)
get_rt(tci::OCCallInfo) = get_rt(tci.ci)

# Special handling for ReturnTypeCall
struct ReturnTypeCallInfo <: CallInfo
    vmi::CallInfo # virtualized method call
end
get_mi((; vmi)::ReturnTypeCallInfo) = isa(vmi, FailedCallInfo) ? nothing : get_mi(vmi)
get_rt((; vmi)::ReturnTypeCallInfo) = Type{isa(vmi, FailedCallInfo) ? Union{} : widenconst(get_rt(vmi))}

struct ConstPropCallInfo <: CallInfo
    mi::CallInfo
    result::InferenceResult
end
get_mi(cpci::ConstPropCallInfo) = cpci.result.linfo
get_rt(cpci::ConstPropCallInfo) = get_rt(cpci.mi)

struct ConstEvalCallInfo <: CallInfo
    mi::CallInfo
    argtypes::ArgTypes
end
get_mi(crci::ConstEvalCallInfo) = get_mi(crci.mi)
get_rt(crci::ConstEvalCallInfo) = get_rt(crci.mi)

# CUDA callsite
struct CuCallInfo <: CallInfo
     cumi::MICallInfo
end
get_mi(gci::CuCallInfo) = get_mi(gci.cumi)
get_rt(gci::CuCallInfo) = get_rt(gci.cumi)

struct Callsite
    id::Int # ssa-id
    info::CallInfo
    head::Symbol
end
get_mi(c::Callsite) = get_mi(c.info)

# Callsite printing
mutable struct TextWidthLimiter <: IO
    io::IO
    width::Int
    limit::Int
end
TextWidthLimiter(io::IO, limit) = TextWidthLimiter(io, 0, limit)

# Uncomment this to debug display of a TextWidthLimiter (e.g., triggered by `print(limited, args1, args...)`)
# Base.show(io::IO, twl::TextWidthLimiter) = error("do not display")

has_space(limiter::TextWidthLimiter, width::Int) = limiter.width + width < limiter.limit - 1
has_space(limiter::TextWidthLimiter, s) = has_space(limiter, textwidth(string(s)))
has_space(::IO, s) = true

function Base.print(io::TextWidthLimiter, s::String)
    io.width == io.limit && return
    width = textwidth(s)
    if has_space(io, width)
        print(io.io, s)
        io.width += width
        return
    end
    for c in graphemes(s)
        cwidth = textwidth(c)
        if has_space(io, cwidth)
            print(io.io, c)
            io.width += cwidth
        else
            break
        end
    end
    print(io.io, '…')
    io.width += 1
    return
end

function Base.print(io::TextWidthLimiter, c::Char)
    tw = textwidth(c)
    if has_space(io, tw)
        print(io.io, c)
        io.width += tw
        return
    end
    return
end

function Base.take!(io::TextWidthLimiter)
    io.width = 0
    return take!(io.io)
end

function headstring(@nospecialize(T))
    if isvarargtype(T)
        T = unwrapva(T)
    elseif isa(T, TypeVar)
        return string(T.name)
    end
    T = widenconst(T)
    if T isa Union || T === Union{}
        return string(T)::String
    elseif T isa UnionAll
        return headstring(Base.unwrap_unionall(T))
    else
        return string(T.name.name)::String
    end
end

function __show_limited(limiter, name, tt, @nospecialize(rt))
    vastring(@nospecialize(T)) = (isvarargtype(T) ? headstring(T)*"..." : string(T)::String)

    if !has_space(limiter, name)
        print(limiter, '…')
        return
    end
    print(limiter, string(name))
    pstrings = String[vastring(T) for T in tt]
    headstrings = String[
        T isa DataType && isempty(T.parameters) ? headstring(T) : string(headstring(T), "{…}")
        for T in tt
    ]
    print(limiter, "(")
    if length(pstrings) != 0
        # See if we have space to print all the parameters fully
        if has_space(limiter, sum(textwidth, pstrings) + 3*length(pstrings))
            join(limiter, (string("::", T) for T in pstrings), ",")
        # Alright, see if we at least have enough space for each head
        elseif has_space(limiter, sum(textwidth, headstrings) + 3*length(pstrings))
            join(limiter, (string("::", T) for T in headstrings), ",")
        # Fine, what about just indicating the number of arguments
        elseif has_space(limiter, 2*(length(tt)))
            join(limiter, ("…" for _ in pstrings), ",")
        else
            print(limiter, "…")
        end
    end
    print(limiter, ")")

    # If we have space for the return type, print it
    rts = string(rt)::String
    if has_space(limiter, textwidth(rts)+2)
        print(limiter, string("::", rts))
    elseif has_space(limiter, 3)
        print(limiter, "::…")
    end
    return
end

function show_callinfo(limiter, mici::MICallInfo)
    mi = mici.mi
    tt = (Base.unwrap_unionall(mi.specTypes)::DataType).parameters[2:end]
    name = (mi.def::Method).name
    rt = get_rt(mici)
    __show_limited(limiter, name, tt, rt)
end

function show_callinfo(limiter, ci::Union{MultiCallInfo, FailedCallInfo, GeneratedCallInfo})
    types = (ci.sig::DataType).parameters
    ft, tt = types[1], types[2:end]
    f = Compiler.singleton_type(ft)
    if f !== nothing
        name = "→ $f"
    elseif ft isa Union
        name = "→ (::Union{$(join(String[String(nameof(T)) for T in Base.uniontypes(ft)], ", "))})"
    else
        name = "→ (::$(nameof(ft)))"
    end
    __show_limited(limiter, name::String, tt, get_rt(ci))
end

function show_callinfo(limiter, (; argtypes, rt)::PureCallInfo)
    ft, tt... = argtypes
    f = Compiler.singleton_type(ft)
    name = isnothing(f) ? "unknown" : string(f)
    __show_limited(limiter, name, tt, rt)
end

function show_callinfo(limiter, ci::ConstPropCallInfo)
    # XXX: The first argument could be const-overriden too
    name = ci.result.linfo.def.name
    tt = ci.result.argtypes[2:end]
    __show_limited(limiter, name, tt, get_rt(ignorewrappers(ci.mi)::MICallInfo))
end

function show_callinfo(limiter, ci::ConstEvalCallInfo)
    # XXX: The first argument could be const-overriden too
    name = get_mi(ci).def.name
    tt = ci.argtypes[2:end]
    __show_limited(limiter, name, tt, get_rt(ci))
end

function show_callinfo(limiter, (; vmi)::ReturnTypeCallInfo)
    if isa(vmi, FailedCallInfo)
        ft = Base.tuple_type_head(vmi.sig)
        f = Compiler.singleton_type(ft)
        name = isnothing(f) ? "unknown" : string(f)
        tt = Base.tuple_type_tail(vmi.sig).parameters
        __show_limited(limiter, name, tt, vmi.rt)
    else
        show_callinfo(limiter, vmi)
    end
end

function Base.show(io::IO, c::Callsite)
    limit = get(io, :limit, false)::Bool
    cols = limit ? (displaysize(io)::Tuple{Int,Int})[2] : typemax(Int)
    optimize = get(io, :optimize, true)::Bool
    iswarn = get(io, :iswarn, false)::Bool
    info = c.info
    rt = get_rt(info)
    if iswarn && is_type_unstable(rt)
        printstyled(io, '%'; color=:red)
    else
        print(io, '%')
    end
    limiter = TextWidthLimiter(io, cols)
    limiter.width += 1   # for the '%' character
    print(limiter, string(c.id, ' '))
    if isa(info, MICallInfo)
        optimize ? print(limiter, string(" = ", c.head, ' ')) : print(limiter, " = ")
        show_callinfo(limiter, info)
    elseif isa(info, WrappedCallInfo)
        wrapped_callinfo(limiter, info)
        show_callinfo(limiter, ignorewrappers(info))
    elseif isa(info, PureCallInfo)
        print(limiter, " = < pure > ")
        show_callinfo(limiter, info)
    elseif info isa MultiCallInfo
        print(limiter, " = call ")
        show_callinfo(limiter, info)
    elseif info isa FailedCallInfo ||
           info isa GeneratedCallInfo
        print(limiter, " = call ")
        show_callinfo(limiter, info)
    elseif info isa TaskCallInfo
        print(limiter, " = task < ")
        show_callinfo(limiter, info.ci)
        print(limiter, " >")
    elseif info isa InvokeCallInfo
        print(limiter, " = invoke < ")
        show_callinfo(limiter, info.ci)
        print(limiter, " >")
    elseif isa(info, ReturnTypeCallInfo)
        print(limiter, " = return_type < ")
        show_callinfo(limiter, info)
        print(limiter, " >")
    elseif isa(info, CuCallInfo)
        print(limiter, " = cucall < ")
        show_callinfo(limiter, info.cumi)
        print(limiter, " >")
    elseif isa(info, ConstPropCallInfo)
        print(limiter, " = < constprop > ")
        show_callinfo(limiter, info)
    elseif isa(info, ConstEvalCallInfo)
        print(limiter, " = < consteval > ")
        show_callinfo(limiter, info)
    elseif isa(info, OCCallInfo)
        print(limiter, " = < opaque closure call > ")
        show_callinfo(limiter, info.ci)
    end
    return
end

function wrapped_callinfo(limiter, ci::WrappedCallInfo)
    print(limiter, " = < ")
    _wrapped_callinfo(limiter, ci)
    ci = get_wrapped(ci)
    while isa(ci, WrappedCallInfo)
        print(limiter, ", ")
        _wrapped_callinfo(limiter, ci)
        ci = get_wrapped(ci)
    end
    print(limiter, " > ")
end
_wrapped_callinfo(limiter, ::LimitedCallInfo)  = print(limiter, "limited")
_wrapped_callinfo(limiter, ::UncachedCallInfo) = print(limiter, "uncached")

is_callsite(cs::Callsite, mi::MethodInstance) = is_callsite(cs.info, mi)
is_callsite(info::MICallInfo, mi::MethodInstance) = get_mi(info) === mi
is_callsite(info::WrappedCallInfo, mi::MethodInstance) = is_callsite(get_wrapped(info), mi)
is_callsite(info::ConstPropCallInfo, mi::MethodInstance) = is_callsite(info.mi, mi)
is_callsite(info::TaskCallInfo, mi::MethodInstance) = is_callsite(info.ci, mi)
is_callsite(info::InvokeCallInfo, mi::MethodInstance) = is_callsite(info.ci, mi)
is_callsite(info::ReturnTypeCallInfo, mi::MethodInstance) = is_callsite(info.vmi, mi)
is_callsite(info::CuCallInfo, mi::MethodInstance) = is_callsite(info.cumi, mi)
function is_callsite(info::MultiCallInfo, mi::MethodInstance)
    for csi in info.callinfos
        is_callsite(csi, mi) && return true
    end
    return false
end
is_callsite(::CallInfo, mi::MethodInstance) = false


text:
if basename(pwd()) == "aoc"
    cd("2015/11")
end
include("../../helpers/helpers.jl")

fromdigits(n; base = 10) = foldl((a, b) -> base * a + b, n, init = 0)

alphabet = "abcdefghjkmnpqrstuvwxyz"
function nextpassword(password::AbstractString)
    n = toencoded(password) + 1
    while !isvalidpass(n)
        n += 1
    end
    fromencoded(n)
end

function toencoded(password::AbstractString)::Int
    fromdigits((findfirst(==(c), alphabet::String) - 1 for c ∈ password), base = 23)
end

function fromencoded(passnum::Integer)::String
    join(getindex.(alphabet::String, reverse(digits(passnum, base = 23)) .+ 1))
end

isvalidpass(password::AbstractString) = isvalidpass(toencoded(password))
isvalidpass(passnum::Integer) = isvalidpass(reverse(digits(passnum, base = 23)))
function isvalidpass(passdigits::AbstractVector{<:Integer})::Bool
    length(passdigits) == 8 &&
        any(a + 1 == b && b + 1 == c for (a, b, c) ∈ eachcons(passdigits, 3)) &&
        count(a == b for (a, b) ∈ Iterators.unique(eachcons(passdigits, 2))) > 1
end

# @time nextpassword("cqjxjnds")
#   0.064901 seconds (306.64 k allocations: 37.737 MiB)
# "cqjxxyzz"

# @time nextpassword(ans)
#   0.275939 seconds (1.19 M allocations: 146.090 MiB, 10.42% gc time)
# "cqkaabcc"

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# custom help for byrow()
function Docs.getdoc(x::typeof(byrow), y)
    if y == Union{}
        return _get_doc_byrow("default")
    elseif y == Tuple{typeof(sum)}
        return _get_doc_byrow("sum")
    elseif y == Tuple{typeof(mean)}
        return _get_doc_byrow("mean")
    elseif y == Tuple{typeof(all)}
        return _get_doc_byrow("all")
    elseif y == Tuple{typeof(any)}
        return _get_doc_byrow("any")
    elseif y == Tuple{typeof(count)}
        return _get_doc_byrow("count")
    elseif y == Tuple{typeof(prod)}
        return _get_doc_byrow("prod")
    elseif y == Tuple{typeof(isequal)}
        return _get_doc_byrow("isequal")
    elseif y == Tuple{typeof(isless)}
        return _get_doc_byrow("isless")
    elseif y == Tuple{typeof(in)}
        return _get_doc_byrow("in")
    end
end

function _get_doc_byrow(fun; text = byrow_docs_text)
    split_text = split(text, "@@@@")
    loc = findfirst(==(fun), split_text)
    Markdown.parse(split_text[loc+1])
end

byrow_docs_text = """
@@@@default@@@@
    byrow(ds, fun, cols; ...)

Perform a row-wise operation specified by `fun` on selected columns `cols`. Generally,
`fun` can be any function that returns a scalar value for each row.

`byrow` is fine tuned for the following operations, to get extra help for each of them search help for `byrow(fun)`, e.g. `?byrow(sum)`;

# Reduction operations

- `all`
- `any`
- `argmax`
- `argmin`
- `coalesce`
- `count`
- `findfirst`
- `findlast`
- `hash`
- `in`
- `isequal`
- `isless`
- `issorted`
- `join`
- `mapreduce`
- `maximum`
- `mean`
- `minimum`
- `nunique`
- `prod`
- `select`
- `std`
- `sum`
- `var`

# Special operations

- `cummax`
- `cummax!`
- `cummin`
- `cummin!`
- `cumprod`
- `cumprod!`
- `cumsum`
- `cumsum!`
- `fill`
- `fill!`
- `sort`
- `sort!`
- `stdze`
- `stdze!`
@@@@sum@@@@
    byrow(ds, sum, cols = names(ds, Number); [by = identity, threads])

Sum results of calling function `by` on each element of each row of `ds`. If `cols` is not specified, `byrow`
computes sum for all numeric columns in `ds`.

Passing `threads = false` disables multitrheaded computations.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

## Example

```jldoctest
julia> ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, sum, :)
3-element Vector{Float64}:
 3.0
 3.5
 7.0
```
@@@@mean@@@@
    byrow(ds, mean, cols = names(ds, Number); [by = identity, threads])

Compute mean of the results of calling function `by` on each element of each row of `ds`. If `cols` is not specified, `byrow`
computes mean for all numeric columns in `ds`.

Passing `threads = false` disables multitrheaded computations.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

## Example

```jldoctest
julia> ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, mean, :)
3-element Vector{Float64}:
 1.5
 1.75
 3.5
```
@@@@all@@@@
    byrow(ds, all, cols = :; [by = isequal(true), threads, mapformats = false])

Test whether all elements in each row in selected columns are `true`, when `by` is passed, determine whether predicate `by` returns `true` for all elements in the row.

By default, `byrow` uses the actual values for test the elements, however, passing `mapformats = true`
change this to the formatted values.

Each columns in `cols` may have its own `by`. This may be achieved by passing a vector of predicates to `by`.

Passing `threads = false` disables multitrheaded computations.

See [`filter`](@ref) or [`filter!`](@ref)

## Example

```jldoctest
julia> ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, all, :, by = [==(2), >(1)])
3-element Vector{Bool}:
 0
 1
 0
```
@@@@any@@@@
    byrow(ds, any, cols = :; [by = isequal(true), threads, mapformats = false])

Test whether any elements in each row in selected columns is `true`, when `by` is passed, determine whether predicate `by` returns `true` for any elements in the row.

By default, `byrow` uses the actual values for test the elements, however, passing `mapformats = true`
change this to the formatted values.

Each columns in `cols` may have its own `by`. This may be achieved by passing a vector of predicates to `by`.

Passing `threads = false` disables multitrheaded computations.

See [`filter`](@ref) or [`filter!`](@ref)

## Example

```jldoctest
julia> ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, any, :, by = [==(2), >(1)])
3-element Vector{Bool}:
 1
 1
 1
```
@@@@count@@@@
    byrow(ds, count, cols = :; [by = isequal(true), threads])

Count the number of elements in each row for selected columns which the function `by` returns `true`.

Passing `threads = false` disables multitrheaded computations.

## Example

```jldoctest
julia> julia> ds = Dataset(x = [1,2,3], y = [2, 6, 5])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         2
   2 │        2         6
   3 │        3         5

julia> byrow(ds, count, :, by = isodd)
3-element Vector{Int32}:
 1
 0
 2
```
@@@@prod@@@@
    byrow(ds, prod, cols = names(ds, Number); [by = identity, threads])

Return the product of the results of calling function `by` on each element of each row of `ds`. If `cols` is not specified, `byrow`
computes product for all numeric columns in `ds`.

Passing `threads = false` disables multitrheaded computations.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

## Example

```jldoctest
julia>  ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, prod, :)
3-element Vector{Union{Missing, Float64}}:
  2.0
  3.0
 12.0
```
@@@@isequal@@@@
    byrow(ds, isequal, cols; [with = nothing, threads])

Returns a boolean vector which is `true` if all values in the corresponding row are equal (using `isequal`).
 Optionally, a vector of values can be passed view the `with` keyword argument to compare values in selected
 columns with the passed vector.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(isless)`](@ref) and [`byrow(in)`](@ref)

## Examples

```jldoctest
julia> ds = Dataset(x1 = [1,2,3,1,2,3], x2 = [1,2,1,2,1,2])
6×2 Dataset
 Row │ x1        x2
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         1
   2 │        2         2
   3 │        3         1
   4 │        1         2
   5 │        2         1
   6 │        3         2

julia> byrow(ds, isequal, [1,2])
6-element Vector{Bool}:
 1
 1
 0
 0
 0
 0

julia> byrow(ds, isequal, [1,2], with = [2,2,2,3,3,3])
6-element Vector{Bool}:
 0
 1
 0
 0
 0
 0
```
@@@@isless@@@@
    byrow(ds, isless, cols, [with, threads, rev = false, lt = isless])

Return a boolean vector which is true if all values in corresponding row for selected `cols` are less than value given by the `with` keyword argument. A vector, or a column name can be passed via `with`.

Passing `rev = true` returns true if all values are greater than passed values via `with`.

By default, the comparison is done via `isless` function, however, user may change it by passing a function via the `lt` keyword argument. The function passed to `lt` must accept two arguments where it takes its first argument from `cols` and its second argument from `with`. However, if `rev = true` the function passed as `lt` will take its first argument from `with` and its second argument from `cols`. The function passed as `lt` must return `true` or `false`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(isequal)`](@ref) and [`byrow(in)`](@ref)

## Examples

```jldoctest
julia> ds = Dataset(x1 = [1,2,3,1,2,3], x2 = [1,2,1,2,1,2], x3 = 6:-1:1)
6×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1         1         6
   2 │        2         2         5
   3 │        3         1         4
   4 │        1         2         3
   5 │        2         1         2
   6 │        3         2         1

julia> byrow(ds, isless, [1,2], with = :x3)
6-element Vector{Bool}:
 1
 1
 1
 1
 0
 0

julia> byrow(ds, isless, 1:2, with = :x3, lt = (x,y) -> isless(x^2, y))
6-element Vector{Bool}:
 1
 1
 0
 0
 0
 0

julia> ds = Dataset(x1 = [1,2,3,1,2,3], x2 = [1,2,1,2,1,2],
                    x3 = [(1,2), (2,3), (1,2), (4,5), (1,2), (1,2)])
6×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Tuple…?
─────┼──────────────────────────────
   1 │        1         1  (1, 2)
   2 │        2         2  (2, 3)
   3 │        3         1  (1, 2)
   4 │        1         2  (4, 5)
   5 │        2         1  (1, 2)
   6 │        3         2  (1, 2)

julia> byrow(ds, isless, 1:2, with = :x3, lt = in)
6-element Vector{Bool}:
 1
 1
 0
 0
 1
 0
```
@@@@in@@@@
    byrow(ds, in, cols; [item, threads, eq = isequal])

Return a boolean vector which its elements are true if in a row the value of `item` is equal to any values from `cols`. The equality is checked via the function passed as `eq`. User can pass a vector of values or a column name to `item`.

The function passed as `eq` must accept two arguments where it takes its first argument from `item` and its second argument from `cols`. The function passed as `eq` must return `true` or `false`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(isequal)`](@ref) and [`byrow(isless)`](@ref)

## Examples

```jldoctest
julia> ds = Dataset(x1 = [1,2,3,1,2,3], x2 = [1,2,1,2,1,2], x3 = 6:-1:1)
6×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1         1         6
   2 │        2         2         5
   3 │        3         1         4
   4 │        1         2         3
   5 │        2         1         2
   6 │        3         2         1

julia> byrow(ds, in, r"x", item = [1,2,3,4,5,6])
6-element Vector{Bool}:
 1
 1
 1
 0
 0
 0

julia> byrow(ds, in, [2,3], item = :x1, eq = isless)
6-element Vector{Bool}:
 1
 1
 1
 1
 0
 0
julia> byrow(ds, in, r"x", item = [5,4,5,4,5,4], eq = (x,y) -> x+y == 11)
6-element Vector{Bool}:
 1
 0
 0
 0
 0
 0
```
"""


text:
using Test: @testset, @test, @inferred
using HankelTensors: Hankel, hankel

@testset "hankel" begin
    channel_size = (3,2)
    input_size = (5,5)
    kernel_size = (3,2)
    batch_size = (2,)
    v = randn(channel_size..., input_size..., batch_size...)
    A = @inferred Hankel(v, channel_size, kernel_size)
    @test (@inferred size(A)) == (channel_size..., kernel_size..., (input_size .- kernel_size .+ 1)..., batch_size...)
    @test A == @inferred hankel(v, channel_size, kernel_size)
    @test ndims(A) == ndims(v) + length(kernel_size)
    for c in CartesianIndices(channel_size), j in CartesianIndices(kernel_size), k in CartesianIndices(input_size .- kernel_size .+ 1), b in CartesianIndices(batch_size)
        i = CartesianIndex(Tuple(j) .+ Tuple(k) .- 1)
        @test A[c,j,k,b] == v[c,i,b]
        #@inferred A[c,j,k,b]  # https://github.com/JuliaLang/julia/issues/44059
    end
end


text:
function ODMMetaData(arg0::ODMFile)
    return ODMMetaData((ODMFile,), arg0)
end

function get_center_body(obj::ODMMetaData)
    return jcall(obj, "getCenterBody", CelestialBody, ())
end

function get_center_name(obj::ODMMetaData)
    return jcall(obj, "getCenterName", JString, ())
end

function get_comment(obj::ODMMetaData)
    return jcall(obj, "getComment", List, ())
end

function get_frame(obj::ODMMetaData)
    return jcall(obj, "getFrame", Frame, ())
end

function get_frame_epoch(obj::ODMMetaData)
    return jcall(obj, "getFrameEpoch", AbsoluteDate, ())
end

function get_frame_epoch_string(obj::ODMMetaData)
    return jcall(obj, "getFrameEpochString", JString, ())
end

function get_frame_string(obj::ODMMetaData)
    return jcall(obj, "getFrameString", JString, ())
end

function get_has_creatable_body(obj::ODMMetaData)
    return jcall(obj, "getHasCreatableBody", jboolean, ())
end

function get_launch_number(obj::ODMMetaData)
    return jcall(obj, "getLaunchNumber", jint, ())
end

function get_launch_piece(obj::ODMMetaData)
    return jcall(obj, "getLaunchPiece", JString, ())
end

function get_launch_year(obj::ODMMetaData)
    return jcall(obj, "getLaunchYear", jint, ())
end

function get_object_id(obj::ODMMetaData)
    return jcall(obj, "getObjectID", JString, ())
end

function get_object_name(obj::ODMMetaData)
    return jcall(obj, "getObjectName", JString, ())
end

function get_odm_file(obj::ODMMetaData)
    return jcall(obj, "getODMFile", ODMFile, ())
end

function get_ref_frame(obj::ODMMetaData)
    return jcall(obj, "getRefFrame", Frame, ())
end

function get_time_scale(obj::ODMMetaData)
    return jcall(obj, "getTimeScale", TimeScale, ())
end

function get_time_system(obj::ODMMetaData)
    return jcall(obj, "getTimeSystem", CcsdsTimeScale, ())
end

function set_center_body(obj::ODMMetaData, arg0::CelestialBody)
    return jcall(obj, "setCenterBody", void, (CelestialBody,), arg0)
end

function set_center_name(obj::ODMMetaData, arg0::JString)
    return jcall(obj, "setCenterName", void, (JString,), arg0)
end

function set_comment(obj::ODMMetaData, arg0::List)
    return jcall(obj, "setComment", void, (List,), arg0)
end

function set_frame_epoch(obj::ODMMetaData, arg0::AbsoluteDate)
    return jcall(obj, "setFrameEpoch", void, (AbsoluteDate,), arg0)
end

function set_frame_epoch_string(obj::ODMMetaData, arg0::JString)
    return jcall(obj, "setFrameEpochString", void, (JString,), arg0)
end

function set_frame_string(obj::ODMMetaData, arg0::JString)
    return jcall(obj, "setFrameString", void, (JString,), arg0)
end

function set_has_creatable_body(obj::ODMMetaData, arg0::jboolean)
    return jcall(obj, "setHasCreatableBody", void, (jboolean,), arg0)
end

function set_launch_number(obj::ODMMetaData, arg0::jint)
    return jcall(obj, "setLaunchNumber", void, (jint,), arg0)
end

function set_launch_piece(obj::ODMMetaData, arg0::JString)
    return jcall(obj, "setLaunchPiece", void, (JString,), arg0)
end

function set_launch_year(obj::ODMMetaData, arg0::jint)
    return jcall(obj, "setLaunchYear", void, (jint,), arg0)
end

function set_object_id(obj::ODMMetaData, arg0::JString)
    return jcall(obj, "setObjectID", void, (JString,), arg0)
end

function set_object_name(obj::ODMMetaData, arg0::JString)
    return jcall(obj, "setObjectName", void, (JString,), arg0)
end

function set_ref_frame(obj::ODMMetaData, arg0::Frame)
    return jcall(obj, "setRefFrame", void, (Frame,), arg0)
end

function set_time_system(obj::ODMMetaData, arg0::CcsdsTimeScale)
    return jcall(obj, "setTimeSystem", void, (CcsdsTimeScale,), arg0)
end


text:
using YaoBase

export ShiftGate, shift

"""
    ShiftGate <: PrimitiveBlock

Phase shift gate.

# Definition

```math
\\begin{pmatrix}
1 & 0\\
0 & e^(im θ)
\\end{pmatrix}
```
"""
mutable struct ShiftGate{T} <: PrimitiveBlock{2}
    theta::T
end

ShiftGate(theta::Integer) = ShiftGate(Float64(theta))
nqudits(blk::ShiftGate) = 1

"""
    shift(θ)

Create a [`ShiftGate`](@ref) with phase `θ`.

# Example

```jldoctest; setup=:(using YaoBlocks)
julia> shift(0.1)
shift(0.1)
```
"""
shift(θ) = ShiftGate(θ)
mat(::Type{T}, gate::ShiftGate) where {T} = Diagonal(T[1.0, exp(im * gate.theta)])

cache_key(gate::ShiftGate) = gate.theta

# parametric interface
niparams(::Type{<:ShiftGate}) = 1
getiparams(x::ShiftGate) = x.theta
setiparams!(r::ShiftGate, param::Number) = (r.theta = param; r)
setiparams(r::ShiftGate, param::Number) = ShiftGate(param)

Base.adjoint(blk::ShiftGate) = ShiftGate(-blk.theta)
Base.copy(block::ShiftGate{T}) where {T} = ShiftGate{T}(block.theta)
Base.:(==)(lhs::ShiftGate, rhs::ShiftGate) = lhs.theta == rhs.theta

# fallback to matrix method if it is not real
YaoBase.isunitary(r::ShiftGate{<:Real}) = true

function YaoBase.isunitary(r::ShiftGate)
    isreal(r.theta) && return true
    @warn "θ in ShiftGate is not real, got θ=$(r.theta), fallback to matrix-based method"
    return isunitary(mat(r))
end

iparams_range(::ShiftGate{T}) where {T} = ((zero(T), T(2.0 * pi)),)


text:
"""
    CouplingVariableKey

Key to map oupling variables. `block_id` identifies the problem block containing
the variable, and `coupling_id` identifies a set of variables whose values are equal.
"""
struct CouplingVariableKey
    block_id::Int
    coupling_id
end

"""
    CouplingVariableRef

Coupling variable reference with `key::CouplingVariableKey`.
"""
struct CouplingVariableRef
    key::CouplingVariableKey
    ref::JuMP.VariableRef
    function CouplingVariableRef(key::CouplingVariableKey, ref::JuMP.VariableRef)
        return new(key, ref)
    end
end

CouplingVariableRef(block_id::Int, coupling_id, ref::JuMP.VariableRef) = CouplingVariableRef(CouplingVariableKey(block_id, coupling_id), ref)

"""
    BlockModel

Block model struture contrains a set of `JuMP.Model` objects, each of which
represents a sub-block model with the information of how to couple these block
models.
"""
abstract type AbstractBlockModel end

mutable struct BlockModel <: AbstractBlockModel
    model::Dict{Int,JuMP.Model} # Dictionary of block models
    coupling_variables::Vector{CouplingVariableRef} # array of variables that couple block models
    variables_by_couple::Dict{Any,Vector{CouplingVariableKey}} # maps `couple_id` to `CouplingVariableKey`

    dual_bound::Float64
    dual_solution::Vector{Float64}

    # TODO: These may be available with heuristics.
    primal_bound::Float64
    primal_solution::Dict{Int, Float64} #coupling_id : value
    combined_weights::Dict{Int, Float64} # block_id : value
    record::Dict{Any, Any}

    function BlockModel()
        return new(
            Dict(),
            [],
            Dict(),
            0.0,
            [],
            +Inf,
            Dict(),
            Dict(),
            Dict())
    end
end

"""
    add_block_model!

Add block model `model` to `block_model::AbstractBlockModel` with `block_id`.
"""
function add_block_model!(block_model::AbstractBlockModel, block_id::Integer, model::JuMP.Model)
    block_model.model[block_id] = model
end

"""
    num_blocks

Number of blocks in `block_model::AbstractBlockModel`
"""
num_blocks(block_model::AbstractBlockModel) = length(block_model.model)

"""
    block_model

This returns a dictionary of `JuMP.Model` objects.
"""
block_model(block_model::AbstractBlockModel) = block_model.model

"""
    block_model

This returns a `JuMP.Model` object for a given `block_id`.
"""
block_model(block_model::AbstractBlockModel, block_id::Integer) = block_model.model[block_id]

"""
    has_block_model

This returns true if `block_model::AbstractBlockModel` has key `block_id::Integer`; false otherwise.
"""
has_block_model(block_model::AbstractBlockModel, block_id::Integer) = haskey(block_model.model, block_id)

"""
    num_coupling_variables

This returns the number of coupling variables in `block_model::AbstractBlockModel`.
"""
num_coupling_variables(block_model::AbstractBlockModel) = length(block_model.coupling_variables)

"""
    coupling_variables

This returns the array of coupling variables in `block_model::AbstractBlockModel`.
"""
coupling_variables(block_model::AbstractBlockModel) = block_model.coupling_variables

"""
    set_coupling_variables!

This sets coupling variables `variables` to `block_model::AbstractBlockModel`.
"""
function set_coupling_variables!(block_model::AbstractBlockModel, variables::Vector{CouplingVariableRef})
    block_model.coupling_variables = variables
end

"""
    set_variables_by_couple!

This sets `BlockModel.variables_by_couple`.
"""
function set_variables_by_couple!(block_model::AbstractBlockModel, variables::Vector{CouplingVariableKey})
    block_model.variables_by_couple = Dict(v.coupling_id => [] for v in variables)
    for v in variables
        push!(block_model.variables_by_couple[v.coupling_id], v)
    end
end

"""
    set_block_weights!
    set the weights of the blocks that will be used by the primal heuristics
"""
function set_block_weights!(block_model::AbstractBlockModel, weights::Dict{Int, Float64})
    block_model.combined_weights = weights
end

"""
    dual_objective_value

This returns the dual objective value obtained from a method.
"""
dual_objective_value(block_model::AbstractBlockModel) = block_model.dual_bound

"""
    dual_solution

This returns a vector of dual solution obtained from a method.
"""
dual_solution(block_model::AbstractBlockModel) = block_model.dual_solution

"""
primal_objective_value

This returns the best primal objective value obtained from a method.
"""
primal_objective_value(block_model::AbstractBlockModel) = block_model.primal_bound

"""
primal_solution

This returns a vector of the best primal solution obtained from a method.
"""
primal_solution(block_model::AbstractBlockModel) = block_model.primal_solution


text:
doc =
"""
usage:
simulation.jl <outputdir> <builddir> <refine> <scheme> <metric_tensor> <run> <vtk>

     builddir path      Directory to write discretization.
     outputdir path     Directory to write output data to.
     refine int         Level of grid refinement to use. refine = 0 (no grid
                                                                  refinement)
     scheme string      Choose between 'staggered' or 'collocated'
     metric_tensor string  Set to 'modified' to use the modified metric tensor
                           discretization.
     run     int        1 (Run simulation in Julia) 0 (Otherwise)
     vtk     int        How often to write vtk files

This example code solves the acoustic wave equation in covariant form on a
curvilinear staggered grid

The governing equations are discretized in space and the zero pressure boundary
condition is applied on all boundaries. The resulting discretization is stored
in a sparse matrix `A`. The solution is advanced in time by solving the linear
system

dq/dt = A*q

using a low-storage Runge-Kutta time integration scheme.
"""

"""
Parse command line arguments
"""

if length(ARGS) != 7
        println(doc)
        exit(-1)
end

builddir=string(ARGS[1], "/")
outputdir=string(ARGS[2], "/")
refine=parse(Int32, ARGS[3])
scheme=ARGS[4]
metric_tensor=ARGS[5]
run=parse(Int32, ARGS[6])
vtk_stride=parse(Int32, ARGS[7])

"""
User settings

nx, ny: Number of grid points in each direction (nodes)
dt: Time step size
num_steps: Number of time steps to take
info_stride: How often to display information
outputdir: Save vtk snapshots to outputdir/{p, v1, v2}_%d.vtk
r1_s, r2_s: Source location in parameter space
r1_r, r2_r: Receiver location in parameter space
num_moment: Number of moment conditions used in source discretization
num_smoothness Number of smoothness conditions used in source discretization
run: Run simulation in Julia
export: Export the discretization for running with the C++ code
"""

ref=2^refine
mx=64*ref
my=32*ref
nx=mx + 1
ny=my + 1
dt=2/mx
num_steps=250*ref
info_stride=100*ref
vtk_stride=vtk_stride*ref
num_moment=4
num_smoothness=4
r1_s=0.5
r2_s=0.5
r1_r=0.4
r2_r=0.4
export_all=true

"""
Define the mapping function x = X(r^1, r^2), y = Y(r^1, r^2) that describes the
geometry of the computational domain.
"""
function mapping(r1, r2)
        x = 10.0 .* r1
        sigma = 10/1.05
        y = (5 .+ 1 .* exp.(-sigma^2*(r1 .- 0.5).^2)) .* r2
        return x, y
end

"""
Define the initial conditions for p, v1, v2. Make sure that p is of the same
size as the grid coordinate vectors (xp, yp), and similarly for v1, v2.
"""
function initial_conditions(xp, yp, x1, y1, x2, y2)
        p = zeros(size(r1_p))
        v1 = zeros(size(r1_1))
        v2 = zeros(size(r1_2))
        return p, v1, v2
end

"""
Define the source time function that introduces a time dependent forcing in the
pressure equation via a singular source term.
"""
function source_time_function(t)
        return exp.( -5e2 .* (t  .- 2.0).^2)
end

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

using Printf
using SparseArrays
import sbp
using sbp.Grid: grid_xp, grid_2d_x, grid_2d_y
using sbp.VTK: vtk_write
using sbp.LSRK4:lsrk4

if scheme == "staggered"
        println("Using staggered discretization")
        ac = sbp.StaggeredAcoustic
else
        println("Using collocated discretization")
        ac = sbp.CollocatedAcoustic
end

# Number of grid points and grids
if scheme == "staggered"
        nxp, nyp = ac.grid_points("p", nx, ny)
        nx1, ny1 = ac.grid_points("v1", nx, ny)
        nx2, ny2 = ac.grid_points("v2", nx, ny)
        r1_p, r2_p = ac.grids("p", nx, ny)
        r1_1, r2_1 = ac.grids("v1", nx, ny)
        r1_2, r2_2 = ac.grids("v2", nx, ny)
        r1, r2 = ac.grids("node", nx, ny)
else
        nxp = nx
        nyp = ny
        nx1 = nx
        ny1 = ny
        nx2 = nx
        ny2 = ny
        r1_p, r2_p = ac.grid(nx, ny)
        r1_1, r2_1 = ac.grid(nx, ny)
        r1_2, r2_2 = ac.grid(nx, ny)
        r1, r2 = ac.grid(nx, ny)
end

x, y = mapping(r1, r2)
xp, yp = mapping(r1_p, r2_p)
x1, y1 = mapping(r1_1, r2_1)
x2, y2 = mapping(r1_2, r2_2)

# Pack initial conditions into `q = [p, v1, v2]`
p, v1, v2 = initial_conditions(xp, yp, x1, y1, x2, y2)
q = vcat(p, v1, v2)
sizes = [0, size(p,1), size(v1,1), size(v2,1)]
ranges = accumulate(+, sizes)

"""

Discretize in space using either the staggered or collocated scheme

"""

@time begin
println("Discretizing in space...")
@printf("refine = %d mx = %d my = %d dt = %g \n", refine, mx, my, dt)
if scheme == "staggered"
        ops = ac.init_operators(nx, ny, sbp.OP2019.build_operators, x, y)
else
        st = sbp.Strand1994
        desc_D, desc_H = sbp.Strand1994.load_description(order=4)
        ops = ac.init_operators(nx, ny, st, desc_D, desc_H, x, y)
end
if metric_tensor == "modified" && scheme == "staggered"
        println("using modified metric tensor discretization")
        Gp = ac.modified_contravariant_metric_tensor(ops)
else
        println("using positivity preserving metric tensor discretization")
        Gp = ac.contravariant_metric_tensor(ops)
end

Ap = ac.divergence(ops)
Av = ac.gradient(ops, Gp)
S = ac.pressure_sat(ops, Gp)
A = -ac.spatial_discretization(Ap, Av - S)
d = ac.pressure_source(ops, r1_s, r2_s, num_moment, num_smoothness)
r = ac.pressure_receiver(ops, r1_r, r2_r, num_moment, num_smoothness)
xs, ys = mapping(r1_s, r2_s)
xr, yr = mapping(r1_r, r2_r)
@printf("Source location: (%f, %f) \n", xs, ys)
@printf("Receiver location: (%f, %f) \n", xr, yr)

end

A = dropzeros(A)

println("Number of non-zeros in A: ", length(A.nzval))

"""
Export to C++ code

"""

if export_all == true
        cfg = Dict("dt" => dt,
                   "nxp" => nxp,
                   "nyp" => nyp,
                   "ninfo" => 20*ref,
                   "ncsv" => ref,
                   "nvtk" => vtk_stride,
                   "ndump" => num_steps,
                   "nrecv" => 1,
                   "nt" => num_steps,
                   "outputdir" => outputdir,
                   )

        println("Writing discretization to ", builddir)
        sbp.lmv.write_config(string(builddir,  "config.txt"), cfg,
                                                               skip_types=true)
        sbp.lmv.write_vector(string(builddir, "q.bin"), q)
        Aout = convert(SparseMatrixCSC{Float64, Int32}, A)
        sbp.lmv.write_csr_matrix(string(builddir, "A.bin"), Aout)
        r = Array(r)[:]
        sbp.lmv.write_vector(string(builddir, "r.bin"), r)
        d = Array(d)[:]
        sbp.lmv.write_vector(string(builddir, "d.bin"), d)
        sbp.lmv.write_vector(string(builddir, "xp.bin"), xp)
        sbp.lmv.write_vector(string(builddir, "yp.bin"), yp)

end

if run == false
        exit()
end

"""
Advance solution in time

"""

scheme = (u, t) -> A * u + d * source_time_function(t)

dq = zeros(size(q))
curr_vtk = 0
elapsed = 0
@printf("step \t time \t time per step \t elapsed time \n")
for step in 1:num_steps

        tk = (step - 1) * dt

        if step % vtk_stride == 0
                outfile=field -> @sprintf("%s/%s_%d.vtk", outputdir, field,
                                          curr_vtk)
                p = q[1+ranges[1]:ranges[2]]
                v1 = q[1+ranges[2]:ranges[3]]
                v2 = q[1+ranges[3]:ranges[4]]
                vtk_elapsed = @elapsed begin
                        vtk_write(outfile("p"), xp, yp, p, nxp, nyp)
                        vtk_write(outfile("v1"), x1, y1, v1, nx1, ny1)
                        vtk_write(outfile("v2"), x2, y2, v2, nx2, ny2)
                end
                @printf("Wrote: %s in %.3f s \n", outfile("{p, v1, v2}"),
                        vtk_elapsed)
                global curr_vtk += 1
        end

        global elapsed += @elapsed global q = lsrk4(scheme, q, dq, tk, dt)
        if step % info_stride == 0
                @printf("%05d \t %.3f \t %.3f s \t %.3f s \n", step, tk, elapsed
                        / step, elapsed)
        end

end

max_stars_repo_name: ooreilly/sbp

max_issues_repo_name: ooreilly/sbp

max_forks_repo_name: ooreilly/sbp


text:
@testset "Misc/MPolyAbsFact" begin

  f = Hecke.MPolyFact.example(wildanger_field(3, 13)[1], 3, 5)
  @test length(factor_absolute(f)) >= 2

  f = Hecke.MPolyFact.example(cyclotomic_field(4)[1], 3, 5)
  @test length(factor_absolute(f)) >= 2
end


text:
module PipsNlpInterface

#Import Packages
using SparseArrays
using LinearAlgebra
using DataStructures
import MPI

import JuMP
import MathOptInterface
const MOI = MathOptInterface
using Plasmo

include("PipsNlpCInterface.jl")
using .PipsNlpSolver

include("pips_utils.jl")
include("nlp_evaluator.jl")

export pipsnlp_solve

mutable struct PIPSNLPData
    d::Union{JuMP.NLPEvaluator,Nothing}
    n::Int
    m::Int
    local_m::Int                    #rank row dimension
    jacnnz::Int                     #Jacobian num nonzeros
    hessnnz::Int                    #Hessian num nonzeros
    firstIeq::Vector{Int}           #first stage equality constraint jacobian nonzero rows
    firstJeq::Vector{Int}           #first stage equality constraint jacobian nonzero cols
    firstVeq::Vector{Float64}       #first stage equality constraint jacobian values
    secondIeq::Vector{Int}          #second stage equality constraint jacobian nonzero rows
    secondJeq::Vector{Int}          #second stage equality constraint jacobian nonzero cols
    secondVeq::Vector{Float64}      #second stage equality constraint jacobian values
    firstIineq::Vector{Int}         #first stage inequality constraint jacobian nonzero rows
    firstJineq::Vector{Int}
    firstVineq::Vector{Float64}
    secondIineq::Vector{Int}
    secondJineq::Vector{Int}
    secondVineq::Vector{Float64}
    num_eqconnect::Int              #num equalities connecting first to second stage
    num_ineqconnect::Int            #num inequalities connecting first to second stage
    eqconnect_lb::Vector{Float64}   #connection equality lower bounds
    eqconnect_ub::Vector{Float64}   #connection equality upper bounds
    ineqconnect_lb::Vector{Float64} #connection inequality lower bounds
    ineqconnect_ub::Vector{Float64} #connection inequality upper bounds
    eq_idx::Vector{Int}             #equality connect indices
    ineq_idx::Vector{Int}           #inequality connect indices
    firstJeqmat
    secondJeqmat
    firstJineqmat
    secondJineqmat
    linkIeq::Vector{Int}            #equality link jacobian rows
    linkJeq::Vector{Int}            #equality link jacobian cols
    linkVeq::Vector{Float64}        #equality link jacobian values
    linkIineq::Vector{Int}
    linkJineq::Vector{Int}
    linkVineq::Vector{Float64}
    x_sol::Vector{Float64}          #solution
    coreid::Int                     #rank
    loaded::Bool
    local_unsym_hessnnz::Int        #local hessian num nonzeros
end
PIPSNLPData() = PIPSNLPData(nothing,0,0,0,0,0,Int[],Int[], Float64[], Int[], Int[], Float64[],Int[],Int[],Float64[], Int[], Int[],
                            Float64[], 0, 0, Float64[],Float64[],Float64[],Float64[], Int[], Int[],
                            nothing, nothing, nothing, nothing, Int[],Int[], Float64[], Int[], Int[], Float64[],Float64[], 0, false, 0)

function _get_pips_data(node::OptiNode)
    return node.ext[:pips_data]
end

function pipsnlp_solve(graph::OptiGraph) #Assume graph variables and constraints are first stage
    #Check structure: 3 possibilities
    #1.) optigraph without subgraphs (only link constraints, no first stage -> second stage constraints)
    #2.) optigraph with single subgraph (first stage and linkconstraints between optinodes)
    #3.) optigraph with multiple subgraphs (first stage and linkconstraints between subgraphs) Not yet supported

    if !(has_subgraphs(graph))
        first_stage = OptiNode() #empty first stage
        submodels = all_nodes(graph)
    elseif length(graph.subgraphs) == 1
        num_nodes(graph) == 1 || error("Incompatible optigraph structure.  Optigraph contains a single subgraph, but has multiple optinodes in highest level graph.")
        first_stage = getnode(graph,1)
        submodels = all_nodes(graph.subgraphs[1])
    end

    comm = MPI.COMM_WORLD
    if MPI.Comm_rank(comm) == 0
        println("Building Model for PIPS-NLP")
    end

    n_sub_models = length(submodels)

    model_list = [first_stage; submodels]

    for (idx,node) in enumerate(model_list)
        node.ext[:pips_data] = PIPSNLPData()
        constraint_data = get_constraint_data(node)
        node.ext[:constraint_data] = constraint_data
    end
    first_stage_data = _get_pips_data(first_stage)

    #Assume every worker has every linkconstraint that connects across nodes
    #TODO: MPI.ALLGATHER! to setup linkconstraints when they are different across workers
    #TODO: Check whether workers can use zeros for bounds they don't need.
    #TODO: Improve distribute performance.

    #IDENTIFY LINK CONSTRAINTS
    link_connect_eq,link_connect_ineq,linkeqconstraints,linkineqconstraints  = _identify_linkconstraints(graph)

    #FIRST STAGE -> SECOND STAGE LINK CONSTRAINTS
    for link in link_connect_eq
        node = _get_subnode(first_stage,link)::OptiNode
        local_data = _get_pips_data(node)
        firstIeq = local_data.firstIeq
        firstJeq = local_data.firstJeq
        firstVeq = local_data.firstVeq
        secondIeq = local_data.secondIeq
        secondJeq = local_data.secondJeq
        secondVeq = local_data.secondVeq
        push!(local_data.eqconnect_lb, link.set.value)
        push!(local_data.eqconnect_ub, link.set.value)
        local_data.num_eqconnect += 1
        row = local_data.num_eqconnect
        for (var,coeff) in link.func.terms
            if getnode(var) == first_stage
                push!(firstIeq, row)
                push!(firstJeq, var.index.value)
                push!(firstVeq, coeff)
            else
                @assert getnode(var) == node
                push!(secondIeq, row)
                push!(secondJeq, var.index.value)
                push!(secondVeq, coeff)
            end
        end
    end

    for link in link_connect_ineq
        node = _get_subnode(first_stage,link)
        firstIineq = local_data.firstIineq
        firstJineq = local_data.firstJineq
        firstVineq = local_data.firstVineq
        secondIineq = local_data.secondIineq
        secondJineq = local_data.secondJineq
        secondVineq = local_data.secondVineq
        if isa(link.set,MOI.LessThan)
            push!(local_data.ineqconnect_lb, -Inf)
            push!(local_data.ineqconnect_ub, link.set.upper)
        elseif isa(link.set,MOI.GreaterThan)
            push!(local_data.ineqconnect_lb, link.set.lower)
            push!(local_data.ineqconnect_ub, Inf)
        elseif isa(link.set,MOI.Interval)
            push!(local_data.ineqconnect_lb, link.set.lower)
            push!(local_data.ineqconnect_ub, link.set.upper)
        end
        local_data.num_ineqconnect += 1
        row = local_data.num_ineqconnect
        for (var,coeff) in link.func.terms
            if getnode(var) == first_stage
                push!(firstIineq, row)
                push!(firstJineq, var.index.value)
                push!(firstVineq, ind)
            else
                @assert getnode(var) == node
                push!(secondIineq, row)
                push!(secondJineq, var.index.value)
                push!(secondVineq, ind)
            end
        end
    end

    #SECOND STAGE LINK CONSTRAINTS
    nlinkeq = length(linkeqconstraints)
    nlinkineq = length(linkineqconstraints)
    ineqlink_lb = zeros(nlinkineq)
    ineqlink_ub = zeros(nlinkineq)
    eqlink_lb = zeros(nlinkeq)
    eqlink_ub = zeros(nlinkeq)

    #INEQUALITY CONSTRAINTS
    #Inequality bounds
    #NOTE: linkconstraint bounds need to be on every rank.
    for (idx,link) in enumerate(linkineqconstraints)
        row = idx
        if isa(link.set,MOI.LessThan)
            ineqlink_lb[row] = -Inf
            ineqlink_ub[row] = link.set.upper
        elseif isa(link.set,MOI.GreaterThan)
            ineqlink_lb[row] = link.set.lower
            ineqlink_ub[row] = Inf
        elseif isa(link.set,MOI.Interval)
            ineqlink_lb[row] = link.set.lower
            ineqlink_ub[row] = link.set.upper
        end
        #Inequality values
        for (var,coeff) in link.func.terms
            node = getnode(var)
            local_data = _get_pips_data(node)
            linkIineq = local_data.linkIineq
            linkJineq = local_data.linkJineq
            linkVineq = local_data.linkVineq
            push!(linkIineq, row)
            push!(linkJineq, var.index.value)
            push!(linkVineq, coeff)
        end
    end

    #EQUALITY CONSTRAINTS
    for (idx,link) in enumerate(linkeqconstraints)
        row = idx
        #Equality bounds
        eqlink_lb[row] = link.set.value
        eqlink_ub[row] = link.set.value

        #Equality values
        for (var,coeff) in link.func.terms
            node = getnode(var)
            local_data = _get_pips_data(node)
            linkIeq = local_data.linkIeq
            linkJeq = local_data.linkJeq
            linkVeq = local_data.linkVeq
            push!(linkIeq, row)              #the variable row
            push!(linkJeq, var.index.value)  #the variable column
            push!(linkVeq, coeff)            #the coefficient
        end
    end

    #Add constraint data for each node to help with function evaluations

    #Initialize
    function str_init_x0(nodeid, x0)
        node = model_list[nodeid+1]
        jump_initval = copy(JuMP.start_value.(JuMP.all_variables(node)))
        local_initval = Vector{Float64}(undef,JuMP.num_variables(node))

        nothing_values = isa.(jump_initval,Nothing)
        float_values = .!(nothing_values)
        local_initval[nothing_values] .= 0  #set to 1 by default
        local_initval[float_values] .= jump_initval[float_values]
        original_copy(local_initval,x0)
    end

    #Problem Info
    function str_prob_info(nodeid,flag,mode,col_lb,col_ub,row_lb,row_ub)
        if flag != 1
            node = model_list[nodeid+1]
            local_data = _get_pips_data(node)
            #Do a warm start on the node
            if !(local_data.loaded)
                local_data.loaded = true
                if (nodeid > 0)
                    if haskey(node.ext, :warmStart)
    	                if node.ext[:warmStart] == true
              	      	   solve(node)
                   	   	end
                    end
            	end

    			#nlp_lb, nlp_ub = JuMP.constraintbounds(node)  #This is every constraint in the model
                nlp_lb, nlp_ub = constraintbounds(node.ext[:constraint_data])       #This is every constraint in the model
         		local_data.local_m  = length(nlp_lb)          #number of local constraints (rows)

    			newRowId = Array{Int}(undef,local_data.local_m)
    			eqId = 1
    			ineqId = 1

                #Go through all constraints, check if they are equality of inequality
                for c in 1:local_data.local_m
                    if nlp_lb[c] == nlp_ub[c]
                        push!(local_data.eq_idx, c)
                        newRowId[c] = eqId
                        eqId +=  1
                    else
        				push!(local_data.ineq_idx, c)
        				newRowId[c] = ineqId
        				ineqId += 1
                    end
         		end

                #Local node data
                local_data.m  = local_data.local_m + local_data.num_eqconnect + local_data.num_ineqconnect
                local_data.n = JuMP.num_variables(node)
                local_data.x_sol = zeros(Float64,local_data.n)
                local_data.firstJeqmat = sparse(local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_data.num_eqconnect, first_stage_data.n)
                local_data.secondJeqmat = sparse(local_data.secondIeq, local_data.secondJeq, local_data.secondVeq, local_data.num_eqconnect, local_data.n)
                local_data.firstJineqmat = sparse(local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_data.num_ineqconnect, first_stage_data.n)
                local_data.secondJineqmat = sparse(local_data.secondIineq, local_data.secondJineq, local_data.secondVineq, local_data.num_ineqconnect, local_data.n)

                if node.nlp_data == nothing
                    JuMP._init_NLP(node)
                end
                local_data.d = JuMP.NLPEvaluator(node)
                MOI.initialize(local_data.d, [:Grad,:Jac, :Hess])
                #Ijac, Jjac = jac_structure(local_data.d)
                jac_structure = pips_jacobian_structure(local_data.d)
                Ijac = [jac_structure[i][1] for i = 1:length(jac_structure)]
                Jjac = [jac_structure[j][2] for j = 1:length(jac_structure)]
                Ijaceq = Int[]
                Jjaceq = Int[]
                Ijacineq = Int[]
                Jjacineq = Int[]
                jac_eq_index = Int[]
                jac_ineq_index = Int[]
                for i in 1:length(Ijac)
                    c = Ijac[i]
                    if nlp_lb[c] == nlp_ub[c]
                        modifiedrow = newRowId[c]
                        push!(Ijaceq, modifiedrow)
                        push!(Jjaceq, Jjac[i])
                        push!(jac_eq_index, i)
                    else
                        modifiedrow = newRowId[c]
                        push!(Ijacineq, modifiedrow)
                        push!(Jjacineq, Jjac[i])
                        push!(jac_ineq_index,i)
                    end
                end

        		node.ext[:Ijaceq] = Ijaceq
         		node.ext[:Jjaceq] = Jjaceq
                node.ext[:Ijacineq] = Ijacineq
         		node.ext[:Jjacineq] = Jjacineq
         		node.ext[:jac_eq_index] = jac_eq_index
         		node.ext[:jac_ineq_index] = jac_ineq_index
         		#Ihess, Jhess = hesslag_structure(local_data.d)
                hess_structure = pips_hessian_lagrangian_structure(local_data.d)
                Ihess = [hess_structure[i][1] for i = 1:length(hess_structure)]
                Jhess = [hess_structure[j][2] for j = 1:length(hess_structure)]

         		Hmap = Bool[]
         		node_Hrows = Int[]
         		node_Hcols = Int[]
         		for i in 1:length(Ihess)
            	    if Jhess[i] <= Ihess[i]
               	        push!(node_Hrows, Ihess[i])
                   		push!(node_Hcols, Jhess[i])
                   		push!(Hmap, true)
            	    else
                        push!(Hmap, false)
                    end
                end
         		val = ones(Float64,length(node_Hrows))
         		mat = sparse(node_Hrows,node_Hcols,val, local_data.n, local_data.n)
         		node.ext[:Hrows] = node_Hrows
         		node.ext[:Hcols] = node_Hcols
         		node.ext[:Hmap] = Hmap
         		local_hessnnz = length(mat.rowval)
         		local_data.local_unsym_hessnnz = length(Ihess)
         		local_data.hessnnz = local_hessnnz
		    end

	 	    if mode == :Values
               	#nlp_lb, nlp_ub = JuMP.constraintbounds(node)
                nlp_lb, nlp_ub = constraintbounds(node.ext[:constraint_data])
    			eq_lb=Float64[]
    			eq_ub=Float64[]
    			ineq_lb=Float64[]
    			ineq_ub=Float64[]
    			for i in 1: length(nlp_lb)
    			    if nlp_lb[i] == nlp_ub[i]
    			       push!(eq_lb, nlp_lb[i])
    			       push!(eq_ub, nlp_ub[i])
    			    else
    			       push!(ineq_lb, nlp_lb[i])
    			       push!(ineq_ub, nlp_ub[i])
    			    end
    			end

    			if nodeid !=  0 #0 is the master
    			   eq_lb = [eq_lb; local_data.eqconnect_lb]
    			   eq_ub = [eq_ub; local_data.eqconnect_ub]
    			   ineq_lb = [ineq_lb; local_data.ineqconnect_lb]
    			   ineq_ub = [ineq_ub; local_data.ineqconnect_ub]
    			end

    			original_copy([eq_lb;ineq_lb], row_lb)
    			original_copy([eq_ub;ineq_ub], row_ub)
                colLower = variablelowerbounds(node)
                colUpper = variableupperbounds(node)
                original_copy(colLower, col_lb)
    			original_copy(colUpper, col_ub)
            end
		    return (local_data.n,local_data.m)
		else
		    if mode == :Values
                # println("row_lb: ",row_lb)
                # println("row_ub: ",row_ub)
                # println(eqlink_lb)
                # println(eqlink_ub)
                # println(ineqlink_lb)
                # println(ineqlink_ub)
                original_copy([eqlink_lb; ineqlink_lb], row_lb)
                original_copy([eqlink_ub; ineqlink_ub], row_ub)
            end
            return (0, nlinkeq + nlinkineq)
		end
    end

    #x0 is first stage variable values, x1 is local values
    function str_eval_f(nodeid,x0,x1)
    	node = model_list[nodeid+1] #Julia doesn't start index at 0
        local_data = _get_pips_data(node)
        local_d = _get_pips_data(node).d
        if nodeid ==  0
            local_x = x0
        else
            local_x = x1
        end
        #check objective sign
        #local_scl = (node.objSense == :Min) ? 1.0 : -1.0
        local_scl = (JuMP.objective_sense(node) == MOI.MAX_SENSE) ? -1.0 : 1.0
        f = local_scl*pips_eval_objective(local_d,local_x)
        return f
    end

    function str_eval_g(nodeid,x0,x1,new_eq_g, new_inq_g)

        node = model_list[nodeid+1]
        local_data = _get_pips_data(node)
        local_d = _get_pips_data(node).d
        if nodeid ==  0
            local_x = x0
        else
            local_x = x1
        end
        local_g = Array{Float64}(undef,local_data.local_m)
        pips_eval_constraint(local_d, local_g, local_x)

        #evaluates local constraints and connecting constraints with first stage
        new_eq_g[1:end] = [local_g[local_data.eq_idx]; local_data.firstJeqmat*x0+local_data.secondJeqmat*x1]
        new_inq_g[1:end] = [local_g[local_data.ineq_idx]; local_data.firstJineqmat*x0+local_data.secondJineqmat*x1]
	    return Int32(1)
    end

    function str_write_solution(id::Integer, x::Vector{Float64}, y_eq::Vector{Float64}, y_ieq::Vector{Float64})
        node = model_list[id+1]
        local_data = _get_pips_data(node)
        local_data.x_sol = copy(x)

        #TODO: Grab dual values

        #NOTE: This won't work anymore
        # if id == 0
        #     #Add objective value to node
        #     node.ext[:objective] = str_eval_f(id, x, nothing)
        # else
        #     #node.objVal = str_eval_f(id, nothing, x)
        #     node.ext[:objective] = str_eval_f(id, nothing, x)
        # end
        rank = MPI.Comm_rank(comm)
        local_data.coreid = rank
    end

    function str_eval_grad_f(rowid,colid,x0,x1,new_grad_f)
        #println("eval_grad_f")
        node = model_list[rowid+1]
        if rowid == colid
            local_data = _get_pips_data(node)
            local_d = _get_pips_data(node).d
            if colid ==  0
                local_x = x0
            else
                local_x = x1
            end
            local_grad_f = Array{Float64}(undef,local_data.n)
            pips_eval_objective_gradient(local_d, local_grad_f, local_x)
            #eval_grad_f(local_d, local_grad_f, local_x)
            #local_scl = (node.objSense == :Min) ? 1.0 : -1.0
            local_scl = (JuMP.objective_sense(node) == MOI.MAX_SENSE) ? -1.0 : 1.0
            #scale!(local_grad_f,local_scl)
            rmul!(local_grad_f,local_scl)
            original_copy(local_grad_f, new_grad_f)
        elseif colid == 0
            new_grad_f[1:end] .= 0
        else
            @assert(false)
        end
        return Int32(1)
    end

    #NOTE: Why does it subtract 1?
    function array_copy(src,dest)
        @assert(length(src)==length(dest))
        for i in 1:length(src)
            dest[i] = src[i]-1  #NOTE: Is this copying indices to a C vector?
        end
    end

    function original_copy(src,dest)
        @assert(length(src)==length(dest))
        for i in 1:length(src)
            dest[i] = src[i]
        end
    end

    function str_eval_jac_g(rowid,colid,flag, x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
        if flag != 1 #populate parent child structure
            node = model_list[rowid+1]
            local_data = _get_pips_data(node)
            local_m_eq = length(local_data.eq_idx)
            local_m_ineq = length(local_data.ineq_idx)
            if mode == :Structure
                if (rowid == colid)
                    Ieq=[node.ext[:Ijaceq];local_data.secondIeq .+ local_m_eq]
                    Jeq=[node.ext[:Jjaceq];local_data.secondJeq]
                    Veq= ones(Float64, length(Ieq))
                    Iineq=[node.ext[:Ijacineq];local_data.secondIineq .+ local_m_ineq]
                    Jineq=[node.ext[:Jjacineq];local_data.secondJineq]
                    Vineq=ones(Float64, length(Iineq))
                    eqmat = sparse(Ieq, Jeq, Veq, local_m_eq + local_data.num_eqconnect, local_data.n)
                    ineqmat = sparse(Iineq, Jineq, Vineq, local_m_ineq + local_data.num_ineqconnect, local_data.n)
                else
                    eqmat = sparse(local_m_eq .+ local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_m_eq .+ local_data.num_eqconnect, first_stage_data.n)
                    ineqmat = sparse(local_m_ineq .+ local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_m_ineq .+ local_data.num_ineqconnect, first_stage_data.n)
                end
                return(length(eqmat.rowval), length(ineqmat.rowval))
            else #mode = :Values
                if rowid == colid    #evaluate Jacobian for block
                    if colid ==  0   #if it's the root node
                        local_x = x0
                    else
                        local_x = x1
                    end
                    local_values = Array{Float64}(undef,length(node.ext[:Ijaceq])+length(node.ext[:Ijacineq]))
                    #eval_jac_g(local_data.d, local_values, local_x)
                    pips_eval_constraint_jacobian(local_data.d, local_values, local_x)
                    jac_eq_index = node.ext[:jac_eq_index]
                    jac_ineq_index = node.ext[:jac_ineq_index]
                    Ieq=[node.ext[:Ijaceq];local_data.secondIeq .+ local_m_eq]
                    Jeq=[node.ext[:Jjaceq];local_data.secondJeq]
                    Veq=[local_values[jac_eq_index];local_data.secondVeq]
                    Iineq=[node.ext[:Ijacineq];local_data.secondIineq .+ local_m_ineq]
                    Jineq=[node.ext[:Jjacineq];local_data.secondJineq]
                    Vineq=[local_values[jac_ineq_index];local_data.secondVineq]
                    eqmat = sparseKeepZero(Ieq, Jeq, Veq, local_m_eq + local_data.num_eqconnect, local_data.n)
                    ineqmat = sparseKeepZero(Iineq, Jineq, Vineq, local_m_ineq + local_data.num_ineqconnect, local_data.n)
                else #evaluate
                    eqmat = sparseKeepZero(local_m_eq .+ local_data.firstIeq, local_data.firstJeq, local_data.firstVeq, local_m_eq .+ local_data.num_eqconnect, first_stage_data.n)
                    ineqmat = sparseKeepZero(local_m_ineq .+ local_data.firstIineq, local_data.firstJineq, local_data.firstVineq, local_m_ineq .+ local_data.num_ineqconnect, first_stage_data.n)
                end
                if length(eqmat.nzval) > 0
                    array_copy(eqmat.rowval,e_rowidx)
                    array_copy(eqmat.colptr,e_colptr)
                    original_copy(eqmat.nzval,e_values)
                end
                if length(ineqmat.nzval) > 0
                    array_copy(ineqmat.rowval,i_rowidx)
                    array_copy(ineqmat.colptr,i_colptr)
                    original_copy(ineqmat.nzval, i_values)
                end
            end
        else #populate linkconstraint structure
            node = model_list[rowid+1]
            local_data = _get_pips_data(node)
            linkIeq = local_data.linkIeq
            linkJeq = local_data.linkJeq
            linkVeq = local_data.linkVeq
            linkIineq = local_data.linkIineq
            linkJineq = local_data.linkJineq
            linkVineq = local_data.linkVineq

            if mode == :Structure
                return(length(linkVeq), length(linkVineq))
            else
                eqmat = sparse(linkIeq, linkJeq, linkVeq, nlinkeq, local_data.n)
                ineqmat = sparse(linkIineq, linkJineq, linkVineq, nlinkineq, local_data.n)
                if length(eqmat.nzval) > 0
                    array_copy(eqmat.rowval,e_rowidx)
                    array_copy(eqmat.colptr,e_colptr)
                    original_copy(eqmat.nzval,e_values)
                end
                if length(ineqmat.nzval) > 0
                    array_copy(ineqmat.rowval,i_rowidx)
                    array_copy(ineqmat.colptr,i_colptr)
                    original_copy(ineqmat.nzval, i_values)
                end
            end
        end
        return Int32(1)
    end

    function str_eval_h(rowid,colid,x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
        #println("eval_h")
        node = model_list[colid+1]
        local_data = _get_pips_data(node)
        if mode == :Structure
            if rowid == colid
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_val = ones(Float64,length(node_Hrows))
                mat = sparseKeepZero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                return length(mat.rowval)
            elseif colid == 0
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_Hrows,node_Hcols = exchange(node_Hrows,node_Hcols)
                node_val = ones(Float64,length(node_Hrows))
                mat = sparseKeepZero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                return length(mat.rowval)
            else
                return 0
            end
        else #fill in values
            if rowid == colid
                node_Hmap = node.ext[:Hmap]
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                if colid ==  0
                    local_x = x0
                else
                    local_x = x1
                end
                local_unsym_values = Array{Float64}(undef,local_data.local_unsym_hessnnz)
                node_val = ones(Float64,length(node_Hrows))
                #local_scl = (node.objSense == :Min) ? 1.0 : -1.0
                local_scl = (JuMP.objective_sense(node) == MOI.MAX_SENSE) ? -1.0 : 1.0
                local_m_eq = length(local_data.eq_idx)
                local_m_ineq = length(local_data.ineq_idx)
                local_lambda = zeros(Float64, local_data.local_m)
                for i in 1:local_m_eq
                    local_lambda[local_data.eq_idx[i]] = lambda[i]
                end
                for i in 1:local_m_ineq
                    local_lambda[local_data.ineq_idx[i]] = lambda[i+local_m_eq]
                end

                #eval_hesslag(local_data.d, local_unsym_values, local_x, obj_factor*local_scl, local_lambda)
                pips_eval_hessian_lagrangian(local_data.d, local_unsym_values, local_x, obj_factor*local_scl, local_lambda)
                local_sym_index=1
                for i in 1:local_data.local_unsym_hessnnz
                    if node_Hmap[i]
                        node_val[local_sym_index] = local_unsym_values[i]
                        local_sym_index +=1
                    end
                end
                node_Hrows,node_Hcols = exchange(node_Hrows,node_Hcols)
                mat = sparseKeepZero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                array_copy(mat.rowval,rowidx)
                array_copy(mat.colptr,colptr)
                original_copy(mat.nzval,values)
            elseif colid ==	0
                node_Hrows  = node.ext[:Hrows]
                node_Hcols = node.ext[:Hcols]
                node_val = zeros(Float64,length(node_Hrows))
                node_Hrows,node_Hcols = exchange(node_Hrows,node_Hcols)
                mat = sparseKeepZero(node_Hrows,node_Hcols,node_val, local_data.n, local_data.n)
                array_copy(mat.rowval,rowidx)
                array_copy(mat.colptr,colptr)
                original_copy(mat.nzval,values)
            else
            end
        end
        return Int32(1)
    end

    if !MPI.Initialized()
        MPI.Init()
    end

    comm = MPI.COMM_WORLD
    if (MPI.Comm_rank(comm) == 0)
        t1 = time()
    end

    root = 0
    rank = MPI.Comm_rank(comm)

    #Create PipsModel (The PIPS interface model) and pass all the functions it requires
    model = PipsModel(:Min,0, n_sub_models,str_init_x0, str_prob_info, str_eval_f, str_eval_g, str_eval_grad_f, str_eval_jac_g, str_eval_h,str_write_solution)

    ###################################
    #PIPS-NLP FUNCTIONS
    prob = createProblemStruct(comm, model, true)
    if rank == 0
        println("Created PIPS-NLP Problem Struct")
    end

    MPI.Barrier(comm)
    ret = solveProblemStruct(prob)
    ###################################

    if rank == 0
        println("Solved PIPS-NLP Problem Struct")
    end

    if rank == 0
        println("Timing Results:")
        println("init_x0  ",prob.t_jl_init_x0)
        println("str_init_x0  ", prob.t_jl_str_prob_info)
        println("eval_f  ", prob.t_jl_eval_f)
        println("eval_g0  ",prob.t_jl_eval_g)
        println("eval_grad_f  ", prob.t_jl_eval_grad_f)
        println("eval_jac  ", prob.t_jl_eval_jac_g)
        println("str_eval_jac  ", prob.t_jl_str_eval_jac_g)
        println("eval_h  ",  prob.t_jl_eval_h)
        println("str_eval_h ",  prob.t_jl_str_eval_h)
        println("eval_write_solution  ",  prob.t_jl_write_solution)
        println("PIPS-NLP time:   ",  time() - t1, " (s)")
    end

    #TODO.  Put solution onto optinodes
    for (idx,node) in enumerate(model_list)  #set solution values for each model
        local_data = _get_pips_data(node)
        if idx != 1 #all cores have the first stage
            coreid = zeros(Int, 1)
            sc = MPI.Reduce(local_data.coreid, MPI.SUM, root, comm)

            if rank == root
                coreid[1] = sc
            end

            MPI.Bcast!(coreid, length(coreid), root, comm)

            n = zeros(Int, 1)
            n[1] = local_data.n
            MPI.Bcast!(n,length(n),coreid[1], comm)

            if rank != coreid[1]
                local_data.n = n[1]
                local_data.x_sol = zeros(Float64,local_data.n)
            end
            MPI.Bcast!(local_data.x_sol, local_data.n, coreid[1], comm)
        end
        #node.ext[:colVal] = local_data.x_sol
        vars = MOI.get(node,MOI.ListOfVariableIndices())
        primals = OrderedDict(zip(vars,local_data.x_sol))

        id = graph.id
        src = JuMP.backend(node)
        src.primals[id] = primals
        src.last_solution_id = id

        #TODO set duals

        #TODO: Load results into a graph backend with termination status etc...
    end

    #TODO: Better return status.  Hook up a simple MOI interface that drops the model results into a graph backend CachingOptimizer
    status = :Unknown
    if ret == 0
        status = :Optimal
    elseif ret == 1
        status = :Not_Finished
    elseif	ret == 2
        status = :Maximum_Iterations_Exceeded
    elseif	ret == 3
        status = :Infeasible_Problem_Detected
    elseif ret == 4
        status = :Restoration_needed
    else
        status = :UnknownStatus
    end

    return status
end  #end pips_nlp_solve

end #end module


text:
using Base.Dates
using Base.Test

using MarketData

using TimeSeries

@testset "split" begin

@testset "find methods" begin
    @testset "find returns correct row numbers array" begin
        @test cl[find(cl .> op)].timestamp[1] == Date(2000,1,3)
        @test length(find(cl .> op))          == 244
    end

    @testset "findwhen returns correct Dates array" begin
       @test findwhen(cl .> op)[2]      == Date(2000,1,5)
       @test length(findwhen(cl .> op)) == 244
    end
end

@testset "split date operations" begin
    @testset "from and to correctly subset non-zero and zero-length time arrays" begin
        @test length(from(cl, Date(2001,12,28))) == 2
        @test length(from(cl, Date(2002,1,1)))   == 0
        @test length(from(from(cl, Date(2002,1,1)), Date(2012,1,1))) == 0

        @test length(to(cl, Date(2000,1,4)))     == 2
        @test length(to(cl, Date(1999,1,4)))     == 0
        @test length(to(to(cl, Date(1999,1,4)), Date(1912,1,1))) == 0
    end

    @testset "when method correctly subset" begin
        @test when(cl, day, 4).timestamp[1]              == Date(2000,1,4)
        @test when(cl, dayname, "Friday").timestamp[1]   == Date(2000,1,7)
        @test when(cl, week, 5).timestamp[1]             == Date(2000,1,31)
        @test when(cl, month, 5).timestamp[1]            == Date(2000,5,1)
        @test when(cl, monthname, "June").timestamp[1]   == Date(2000,6,1)
        @test when(cl, year, 2001).timestamp[1]          == Date(2001,1,2)
        @test when(cl, dayofweek, 1).timestamp[1]        == Date(2000,1,3)
        # all the days in the nth week of each month
        @test when(cl, dayofweekofmonth, 5).timestamp[1] == Date(2000,1,31)
        @test when(cl, dayofyear, 365).timestamp[1]      == Date(2001,12,31)
        @test when(cl, quarterofyear, 4).timestamp[1]    == Date(2000,10,2)
        @test when(cl, dayofquarter, 1).timestamp[1]     == Date(2001,10,1)
    end
end

@testset "element wrappers" begin
    @testset "type element wrappers isolate elements" begin
        @test isa(timestamp(cl), Array{Date,1})       == true
        @test isa(values(cl), Array{Float64,1})       == true
        @test isa(values(ohlc), Array{Float64,2})     == true
        @test isa(colnames(cl), Array{String, 1}) == true
    end
end

@testset "head and tail methods" begin
    @testset "head and tail methods work with default n value on single column TimeArray" begin
        @test length(head(cl))      == 1
        @test length(head(cl,1))    == 1
        @test head(cl).timestamp[1] == Date(2000,1,3)
        @test head(cl).values[1]    == 111.94
        @test head(cl).meta         == "AAPL"

        @test length(tail(cl))      == 1
        @test length(tail(cl,1))    == 1
        @test tail(cl).timestamp[1] == Date(2001,12,31)
        @test tail(cl).values[1]    == 21.9
        @test tail(cl).meta         == "AAPL"
    end

    @testset "head and tail methods work with default n value on multi column TimeArray" begin
        @test length(head(ohlc))      == 1
        @test length(head(ohlc,1))    == 1
        @test head(ohlc).timestamp[1] == Date(2000,1,3)
        @test head(ohlc).values       == [104.88  112.5  101.69  111.94]
        @test head(ohlc).meta         == "AAPL"

        @test length(tail(ohlc))      == 1
        @test length(tail(ohlc,1))    == 1
        @test tail(ohlc).timestamp[1] == Date(2001,12,31)
        @test tail(ohlc).values       == [22.51  22.66  21.83  21.9]
        @test tail(ohlc).meta         == "AAPL"
    end

    @testset "head and tail methods work with custom periods on single column TimeArray" begin
        @test length(head(cl, 2))   == 2
        @test length(head(cl, 500)) == length(cl)
        @test length(tail(cl, 2))   == 2
        @test length(tail(cl, 500)) == length(cl)
    end

    @testset "head and tail methods work with custom periods on multi column TimeArray" begin
        @test length(head(ohlc, 2))   == 2
        @test length(head(ohlc, 500)) == length(ohlc)
        @test length(tail(ohlc, 2))   == 2
        @test length(tail(ohlc, 500)) == length(ohlc)
    end
end

end  # @testset "split"


text:
#simple.jl

include("kernels_naive.jl")
include("kernels_pure.jl")

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

#Helper routines that initialize, execute, time, and validate the kernels

function runhpl(n)
    #Initialize
    hpl(rand(1, 1), rand(1)) #Precompile
    A = randn(n, n)
    b = randn(n)
    A′= copy(A)

    #Run
    t = @elapsed x=hpl(A, b)

    #Validate
    r = b - A′*x
    ϵ = eps()
    r₀= norm(r,Inf)
    nrmA1 = norm(A′, 1)
    r₁= r₀ / (ϵ*nrmA1*n)
    r₂= r₀ / (ϵ*nrmA1*norm(x,1))
    r₃= r₀ / (ϵ*norm(A′,Inf)*norm(x,Inf)*n)
    err = max(r₁, r₂, r₃)
    if err ≥ 16
        warn("Error", err, "exceeds allowed value of 16")
    end
    return t
end

function runhpl2(n)
    #Initialize
    hpl2(rand(1, 1), rand(1)) #Precompile
    A = randn(n, n)
    b = randn(n)
    A′= copy(A)

    #Run
    t = @elapsed x=hpl2(A, b)

    #Validate
    r = b - A′*x
    ϵ = eps()
    r₀= norm(r,Inf)
    nrmA1 = norm(A′, 1)
    r₁= r₀ / (ϵ*nrmA1*n)
    r₂= r₀ / (ϵ*nrmA1*norm(x,1))
    r₃= r₀ / (ϵ*norm(A′,Inf)*norm(x,Inf)*n)
    err = max(r₁, r₂, r₃)
    if err ≥ 16
        warn("Error", err, "exceeds allowed value of 16")
    end
    return t
end

function runrandomupdate(m, nupdate=4m, seed=1)
    #Initialize
    randomupdate!(UInt64[1], 1) #Precompile
    srand(seed)
    T = Array{UInt64}(m)
    for i = 1:m
        T[i] = i
    end

    #Run
    t = @elapsed randomupdate!(T, nupdate)

    #Validate
    srand(seed)
    randomupdate!(T, nupdate)
    err = 0
    for i = 1:m
        err += (T[i] != i)
    end
    err /= m
    if err > 0.01
        warn("Error rate = ", 100err, "%, exceeding 1% allowance")
    end

    return t
end

function runstreamtriad(m, ntrials=10)
    #Initialize
    streamtriad!([0.0], [1.0], 2.0, [3.0])
    a = zeros(m)
    b = randn(m)
    c = randn(m)
    α = randn()

    #Run
    t = Inf
    for i in 1:ntrials
        t = min(t, @elapsed streamtriad!(a, b, α, c))
    end
    return t
end

function runfft(m)
    #Initialize
    FFT([1.0+0.0im])
    z = Vector{Complex128}(m)
    for i in 1:m
        z[i] = randn() + im*randn()
    end
    z′ = copy(z)

    #Run
    t = @elapsed FFT(z)

    #Validate
    ifft!(z)
    err = maxabs(z′ - z)
    errtol = 16*eps()*log2(m)
    if err > errtol
        warn("Error = $err exceeds tolerance of $errtol")
    end

    return t
end

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

#Run the tests

println("Serial HPL")
#        1234567890123456789012345678901234567890123456
println("Problem size (n) | Run time (sec.) | Gigaflops")
for e in 3:0.2:3.8
    n = round(Int, 10^e)
    @printf("%16d | ", n)
    t = runhpl(n)
    @printf("%15.6f | %15.9f\n", t, 1e-9*(2n^3/3 + 3n^2/2)/t)
end

println("\n\n\nSerial HPL v2 - pure Julia kernel")
#        1234567890123456789012345678901234567890123456
println("Problem size (n) | Run time (sec.) | Gigaflops")
for e in 3:0.2:3.4
    n = round(Int, 10^e)
    @printf("%16d | ", n)
    t = runhpl2(n)
    @printf("%15.6f | %15.9f\n", t, 1e-9*(2n^3/3 + 3n^2/2)/t)
end

println("\n\n\nSerial random update")
#        12345678901234567890123456789012345678901234567890123
println("Problem size (m) | Run time (sec.) | Gigaupdates/sec.")
for e in 18:24
    m = 2^e
    @printf("%16d | ", m)
    t = runrandomupdate(m)
    @printf("%15.6f | %15.9f\n", t, m*1e-9/t)
end

println("\n\n\nSerial STREAM triad")
#        123456789012345678901234567890123456789012345678901
println("Problem size (m) | Run time (sec.) | Gigabytes/sec.")
for e in 6.5:0.5:8.0
    m = round(Int, 10^e)
    @printf("%16d | ", m)
    t = runstreamtriad(m)
    @printf("%15.6f | %15.9f\n", t, 1e-9*24*10*m/t)
end

println("\n\n\nSerial FFT")
#        123456789012345678901234567890123456789012345678901
println("Problem size (m) | Run time (sec.) | Gigaflops")
for e in 6:0.5:7
    m = round(Int, 10^e)
    @printf("%16d | ", m)
    t = runfft(m)
    @printf("%15.6f | %15.9f\n", t, 5e-9*m*log2(m)/t)
end


text:
using Documenter

const pages = [
    "Home" => "index.md",
    "File" => "file.md",
    "Subdir" => "subdir/index.md",
    "Subfile" => "subdir/file.md",
]

@info "Building builds/default"
makedocs(sitename="Test", pages = pages, build="builds/default")

@info "Building builds/absolute"
mkdir(joinpath(@__DIR__, "builds/absolute-workdir"))
makedocs(sitename="Test", pages = pages, build="builds/absolute", workdir=joinpath(@__DIR__, "builds/absolute-workdir"))

@info "Building builds/relative"
mkdir(joinpath(@__DIR__, "builds/relative-workdir"))
makedocs(sitename="Test", pages = pages, build="builds/relative", workdir="builds/relative-workdir")


text:
# This file was generated, do not modify it. # hide
m = MLJ.machine(tm, features, targets)

MLJ.fit!(m)


text:
T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-4,
    cost_tolerance_intermediate=1.0e-3,
    constraint_tolerance=max_con_viol,
    penalty_scaling=10.,
    penalty_initial=1.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3,
    opts_al=opts_al,
    opts_pn=opts_pn)

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

prob = copy(Problems.pendulum_static)
U0 = deepcopy(controls(prob))
ilqr = iLQRSolver(prob)
initial_controls!(ilqr, U0)
Z = ilqr.Z
@btime TO.cost!($ilqr.obj, $Z)
ilqr.obj[prob.N-1].Q
solve!(ilqr)

@btime begin
    initial_controls!($ilqr, $U0)
    solve!($ilqr)
end

# Solve with AL-iLQR
prob = copy(Problems.pendulum_static)
U0 = deepcopy(controls(prob))
al = AugmentedLagrangianSolver(prob, opts_al)
initial_controls!(al, U0)
solve!(al)
max_violation(al)

@btime begin
    initial_controls!($al, $U0)
    solve!($al)
end

# Solve with ALTRO
prob = copy(Problems.pendulum_static)
U0 = deepcopy(controls(prob))
altro = ALTROSolver(prob, opts_altro)
initial_controls!(altro, U0)
solve!(altro)
max_violation(altro)

@btime begin
    initial_controls!($altro, $U0)
    solve!($altro)
end
max_violation(altro)

# Ipopt
prob_ipopt = copy(Problems.parallel_park_static)
prob_ipopt = TO.change_integration(prob_ipopt, HermiteSimpson)
rollout!(prob_ipopt)
ipopt = DIRCOLSolver(prob_ipopt, opts_ipopt)
ipopt.opts.verbose = false
solve!(ipopt)
max_violation(ipopt)

@btime solve!($ipopt)


text:
using Test, Workbooks

###  W O R K B O O K S
wb = Workbook()
S = wb["Sheet1"]

@testset "references       " begin
    @test Workbooks.LocalPosition(1,1) == Workbooks.LocalPosition("A1")
    @test Workbooks.LocalPosition(1,1) == Workbooks.LocalPosition("A",1)
    @test Workbooks.LocalPosition(1,1) == Workbooks.LocalRef("A1")

    @test Workbooks.GlobalPosition("Sheet1!A1") == Workbooks.GlobalPosition("Sheet1","A1")

    @test_throws AssertionError Workbooks.LocalPosition("A0")
end

@testset "basic formulas   " begin
    # test constant assignment
    S["A1"] = "3"
    @test S["A1"].value == 3
    # testing constant formulas
    S["A2"] = "=5"
    @test S["A2"].value == 5
    # the next value will be rewritten
    S["B1"] = 7
    # test a formula
    S["B2"] = "=1+A1*(B1-A2)" # now = 1 + 3*(7-5) = 7
    @test S["B2"].value == 7
    # testing formulas without inputs
    S["D1"] = "=6*6"
    @test S["D1"].value == 36
    # testing formulas with missing inpuits
    S["C3"] = "=B2/B3"
    @test ismissing(S["C3"].value)
    # supplying the missing input
    S["B3"] = "2"
    @test S["C3"].value == 7/2    # C3 should now be 7/2 = 3.5
    # refreshing an existing input
    S["B1"] = "11"
    @test S["B2"].value == 19     # B2 = 1 + 3*(11-5) = 19
    @test S["C3"].value == 19/2   # C3 = 19/2 = 9.5
    # test repeated references
    S["A3"] = "=1/A1 + 2/A2 + 3/A1"
    @test S["A3"].value == 1/3 + 2/5 + 3/3
    # test simplest reference
    S["C1"] = "=A3"
    @test S["C1"].value == S["A3"].value
    # test a circular reference
    @test_throws ErrorException S["A1"] = "=A1"
    @test_throws ErrorException S["A3"] = "=C1"
end

@testset "evaluating ranges" begin
    # reference to a range
    S["D2"] = "=A1:A3"
    @test all(S["D2"].value .== Real[S["A1"].value, S["A2"].value, S["A3"].value])
    # summing the range
    S["D3"] = "=sum(D2)"
    @test S["D3"].value == sum(getfield.(S["A1:A3"],:value))
    @test S["D3"].value == (9 + 2/5 + 1/3)
    # refreshing with a range reference
    S["A3"] = "7"
    @test S["D3"].value == 3 + 5 + 7
end

@testset "cross-sheet refs " begin
    global S2 = wb["Sheet2"]
    S2["A1"] = "1"
    S2["A2"] = "=Sheet1!A2"
    @test S2["A2"].value == S["A2"].value
end

@testset "sheets           " begin
    @test size(wb["Sheet1"]) == (3,4)
    @test size(S2) == (2,1)
end

### J e W e L s
@testset "JWLs             " begin
    jwl = JWL("test1")
    jwl.wb = wb
    Workbooks.updatefiles(jwl)

    jwl.name = "test2"
    Workbooks.write(jwl)

    jwl_copy = Workbooks.read("test2.jwl")
    wb_copy = jwl_copy.wb
    for (sname,S) in wb_copy.sheets
        @test Workbooks.hassheet(wb,sname)
        if Workbooks.hassheet(wb,sname)
            a = axes(S)
            @test size(S) == size(wb[sname])
            for i in a[1]
                for j in a[2]
                    if ismissing(S[i,j])
                        @test ismissing(wb[sname][i,j])
                    else
                        @test Workbooks.formula(S[i,j]) == Workbooks.formula(wb[sname][i,j])
                    end
                end
            end
            @test all(skipmissing(Workbooks.array(S)) .== skipmissing(Workbooks.array(wb[sname])))
        end
    end

    T_jwl = Workbooks.readabs("test.jwl")
    @test hasmethod("DILaum9y", Tuple{})
    @test DILaum9y() == "JWLs rock"
end