text:
# Data loader and preprocessing funs for dataframes:
#
# (c) A. Dominik, 2020

"""
    dataframe_read(fname)

Read a data table from an CSV-file with one sample
per row and return a DataFrame with the data.
(ODS-support is removed because of PyCall compatibility issues
of the OdsIO package).
"""
function dataframe_read(fname)

    if occursin(r".*\.ods$", fname)
        println("Reading ODS-files is no longer supported!")
        return nothing
        # return readODS(fname)
    elseif occursin(r".*\.csv$", fname)
        return readCSV(fname)
    else
        println("Error: unknown file format!")
        println("Please support csv or ods file")
        return nothing
    end
end

# function readODS(fname)
#
#     printf("Reading data from ODS: $fname")
#     return OdsIO.ods_read(fname, retType="DataFrame")
# end

function readCSV(fname)

    println("Reading data from CSV: $fname")
    return DataFrames.DataFrame(CSV.File(fname, header=true))
end

"""
    dataframe_minibatches(data::DataFrames.DataFrame; size=256, ignore=[], teaching="y", o...)

Make Knet-conform minibatches of type `Knet.data` from a dataframe
with one sample per row.

### Arguments:
+ `ignore`: defines a list of column names to be ignored
+ `teaching="y"`: defines the column name with teaching input. Default is "y".
                `teaching` is handled differently, depending on its type:
                If `Int`, the teaching input is interpreted as
                class ids and directly used for training (this assumes that
                the values range from 1..n). If type is a String, values are
                interpreted as class labels and convertet to numeric class IDs
                by calling `mk_class_ids()`. The list of valid lables and their
                order can be created by calling `mk_class_ids(data.y)[2]`.
                If teaching is a scalar value, regression context is assumed,
                and the value is used unchanged for training.
+ other keyword arguments: all keyword arguments accepted by
                `Knet.minibatch()` may be used.
"""
function dataframe_minibatches(data; size=16, ignore=[], teaching="y", o...)

    if !(ignore isa(AbstractArray))
        ignore = [ignore]
    end
    if teaching != nothing
        push!(ignore, teaching)
    end
    cols = filter(c->!(c in ignore), names(data))
    x = convert2KnetArray(data[!,cols])
    x = permutedims(x)

    if teaching == nothing
        return Knet.minibatch(x, size; o...)
    else
        # care for type of teaching column:
        #
        t_type = eltype(data[!,teaching])
        if t_type <: AbstractString       # make class_ids from labels
            teach = mk_class_ids(data[!,teaching])[1]
            y = permutedims(Array{UInt8}(teach))

        elseif t_type <: Int              # take values as class_ids
            y = permutedims(Array{UInt8}(data[!,teaching]))

        elseif t_type <: Real             # use as is -> probably regression
            y = permutedims(convert2KnetArray(data[!,teaching]))

        else
            println("Don't know how to handle teaching input of type $t_type")
            return nothing
        end
        return Knet.minibatch(x, y, size; o...)
    end
end

"""
    function mk_class_ids(labels)

Take a list with n class labels for n instances and return a list of
n class-IDs (of type Int) and an array of lables with the array index
of each label corresponds its ID.

### Arguments:
+ `labels`: List of labels (typically Strings)

### Result values:
+ array of class-IDs in the same order as the input
+ array of unique class-IDs ordered by their ID.

### Examples:
```
julia> labels = ["blue", "red", "red", "red", "green", "blue", "blue"]
7-element Array{String,1}:
 "blue"
 "red"
 "red"
 "red"
 "green"
 "blue"
 "blue"

julia> mk_class_ids(labels)[1]
7-element Array{Int64,1}:
 1
 3
 3
 3
 2
 1
 1

 julia> mk_class_ids(labels)[2]
3-element Array{String,1}:
 "blue"
 "green"
 "red"
```
"""
function mk_class_ids(labels)

    l = sort(unique(labels))
    ids = [findfirst(x->x==c, l) for c in labels]
    return ids, l
end

"""
    function dataframe_split(df::DataFrames.DataFrame;
                             teaching="y", fr=0.2, balanced=true)

Split data, organised row-wise in a DataFrame into train and valid sets.

### Arguments:
+ `df`: data
+ `teaching="y"`: name or index of column with teaching input (y)
+ `fr=0.2`: fraction of data to be used for validation
+ `shuffle=true`: shuffle the rows of the dataframe.
+ `balanced=true`: if `true`, result datasets will be balanced by oversampling.
              Returned datasets will be bigger as expected
              but include the same numbers of samples for each class.
"""
function dataframe_split(df::DataFrames.DataFrame; teaching="y",
                         fr=0.2, shuffle=true, balanced=false)

    if shuffle
        df .= df[Random.randperm(DataFrames.nrow(df)),:]
    end
    ((trn,ytrn), (vld,yvld)) = do_split(df, df[:,teaching], at=fr)

    if balanced
        (trn,ytrn) = do_balance(trn, ytrn)
        (vld,yvld) = do_balance(vld, yvld)
    end

    return trn, vld
end


text:
##############################################################################
#
# Specialization of AbstractMPCCModel based on automatic differentiation
# Based on AbastractNLPModel:
# https://github.com/JuliaSmoothOptimizers/NLPModels.jl
#
# TODO:
# - export
# - doc
# - list of functions not adapted: jth_hprod, jth_hprod!, ghjvprod, ghjvprod!, jth_con, jth_congrad, jth_congrad!, jth_sparse_congrad,
# - nnzj for G and H
#
##############################################################################

# Methods to be overridden in other packages.
"""
    f = obj(nlp, x)
Evaluate ``f(x)``, the objective function of `mpcc` at `x`.
"""
obj(::AbstractMPCCModel, ::AbstractVector) =
  throw(NotImplementedError("obj"))

"""
    g = grad(nlp, x)
Evaluate ``∇f(x)``, the gradient of the objective function at `x`.
"""
function grad(nlp::AbstractMPCCModel, x::AbstractVector)
  g = similar(x)
  return grad!(nlp, x, g)
end

"""
    g = grad!(nlp, x, g)
Evaluate ``∇f(x)``, the gradient of the objective function at `x` in place.
"""
grad!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector) =
  throw(NotImplementedError("grad!"))

"""
    c = cons(nlp, x)
Evaluate ``c(x)``, the constraints at `x`.
"""
function cons(nlp::AbstractMPCCModel, x::AbstractVector)
  c = similar(x, nlp.meta.ncon)
  return cons!(nlp, x, c)
end

"""
    c = cons!(nlp, x, c)
    Return the vector of the constraints
    lc <= c(x) <= uc,
    lccG <= G(x),
    lccH <= H(x)
"""
function cons!(nlp :: AbstractMPCCModel, x :: AbstractVector, c :: AbstractVector)
  if nlp.meta.ncc > 0
   G, H = consG(nlp, x), consH(nlp, x)
  else
   G, H = Float64[], Float64[]
  end

  if nlp.meta.ncon > 0
   vnl = cons_nl(nlp, x)
  else
   vnl = Float64[]
  end

  c = vcat(vnl, G, H)
  return c
end

"""
    c = viol(nlp, x)
Evaluate ``c(x)``, the constraints at `x`.
"""
function viol(nlp::AbstractMPCCModel, x::AbstractVector)
  c = similar(x, nlp.meta.nvar+nlp.meta.ncon+3*nlp.meta.ncc)
  return viol!(nlp, x, c)
end

"""
    c = viol!(nlp, x, c)
    Return the vector of the constraints
    lx <= x <= ux
    lc <= c(x) <= uc,
    lccG <= G(x),
    lccH <= H(x),
    G(x) .* H(x) <= 0
"""
function viol!(mod :: AbstractMPCCModel, x :: AbstractVector, c :: AbstractVector)
  n, ncon, ncc = mod.meta.nvar, mod.meta.ncon, mod.meta.ncc
  if ncc > 0
   cG, cH = consG(mod, x), consH(mod, x)

   c[n+ncon+1:n+ncon+ncc]     = max.(mod.meta.lccG - cG, 0.0)
   c[n+ncon+ncc+1:n+ncon+2*ncc] = max.(mod.meta.lccH - cH, 0.0)
   c[n+ncon+2*ncc+1:n+ncon+3*ncc] = max.(cH.*cG, 0.0)
  end

  if mod.meta.ncon > 0
   cx = cons_nl(mod, x)
   c[n+1:n+ncon] = max.(max.(mod.meta.lcon - cx, 0.0),max.(cx - mod.meta.ucon, 0.0))
  end

  c[1:n] = max.(max.(mod.meta.lvar - x, 0.0),max.(x - mod.meta.uvar, 0.0))
  return c
end

"""
Evaluate ``c(x)``, the constraints at `x`.
"""
function cons_nl(nlp :: AbstractMPCCModel, x :: AbstractVector)
 c = similar(x, nlp.meta.ncon)
 return cons_nl!(nlp, x, c)
end

cons_nl!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector) =
  throw(NotImplementedError("cons_nl!"))

"""
Evaluate ``G(x)``, the constraints at `x`.
"""
function consG(nlp :: AbstractMPCCModel, x :: AbstractVector)
 c = similar(x, nlp.meta.ncc)
 return consG!(nlp, x, c)
end

consG!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector) =
       throw(NotImplementedError("consG!"))

"""
Evaluate ``G(x)``, the constraints at `x`.
"""
function consH(nlp :: AbstractMPCCModel, x :: AbstractVector)
 c = similar(x, nlp.meta.ncc)
 return consH!(nlp, x, c)
end

consH!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector) =
              throw(NotImplementedError("consH!"))

jth_con(::AbstractMPCCModel, ::AbstractVector, ::Integer) =
  throw(NotImplementedError("jth_con"))

function jth_congrad(nlp::AbstractMPCCModel, x::AbstractVector, j::Integer)
  g = Vector{eltype(x)}(undef, nlp.meta.nvar)
  return jth_congrad!(nlp, x, j, g)
end

jth_congrad!(::AbstractMPCCModel, ::AbstractVector, ::Integer, ::AbstractVector) =
  throw(NotImplementedError("jth_congrad!"))

jth_sparse_congrad(::AbstractMPCCModel, ::AbstractVector, ::Integer) =
  throw(NotImplementedError("jth_sparse_congrad"))

"""
    f, c = objcons(nlp, x)
Evaluate ``f(x)`` and ``c(x)`` at `x`.
"""
function objcons(nlp :: AbstractMPCCModel, x :: AbstractVector)
  c = similar(x, nlp.meta.ncon)
  return objcons!(nlp, x, c)
end

"""
    f = objcons!(nlp, x, c)
Evaluate ``f(x)`` and ``c(x)`` at `x`. `c` is overwritten with the value of ``c(x)``.
"""
function objcons!(nlp :: AbstractMPCCModel, x :: AbstractVector, c :: AbstractVector)
  f = obj(nlp, x)
  nlp.meta.ncon > 0 && cons!(nlp, x, c)
  return f, c
end

"""
    f, g = objgrad(nlp, x)
Evaluate ``f(x)`` and ``∇f(x)`` at `x`.
"""
function objgrad(nlp :: AbstractMPCCModel, x :: AbstractVector)
  g = similar(x)
  return objgrad!(nlp, x, g)
end

"""
    f, g = objgrad!(nlp, x, g)
Evaluate ``f(x)`` and ``∇f(x)`` at `x`. `g` is overwritten with the
value of ``∇f(x)``.
"""
function objgrad!(nlp :: AbstractMPCCModel, x :: AbstractVector, g :: AbstractVector)
  f = obj(nlp, x)
  grad!(nlp, x, g)
  return f, g
end

"""
    (rows,cols) = jac_structure(nlp)
Return the structure of the constraint's Jacobian in sparse coordinate format.
"""
function jac_structure(nlp :: AbstractMPCCModel)
  rows = Vector{Int}(undef, nlp.meta.nnzj+2 * nlp.meta.nvar * nlp.meta.ncc)
  cols = Vector{Int}(undef, nlp.meta.nnzj+2 * nlp.meta.nvar * nlp.meta.ncc)
  return jac_structure!(nlp, rows, cols)
end

"""
    jac_structure!(nlp, rows, cols)
Return the structure of the constraint's Jacobian in sparse coordinate format in place.
"""
function jac_structure!(nlp :: AbstractMPCCModel, rows :: AbstractVector{<:Integer}, cols :: AbstractVector{<:Integer})
 m, n = nlp.meta.ncon+2*nlp.meta.ncc, nlp.meta.nvar
 I = ((i,j) for i = 1:m, j = 1:n)
 rows[1 : n*m] .= getindex.(I, 1)[:]
 cols[1 : n*m] .= getindex.(I, 2)[:]
 return rows, cols
end
#I = findall(!iszero, Hx)
#return (getindex.(I, 1), getindex.(I, 2), Hx[I])

"""
    vals = jac_coord!(nlp, x, vals)
Evaluate ``∇c(x)``, the constraint's Jacobian at `x` in sparse coordinate format,
rewriting `vals`.
"""
function jac_coord!(nlp :: AbstractMPCCModel, x :: AbstractVector, vals ::AbstractVector)
 m, n = nlp.meta.ncon+2*nlp.meta.ncc, nlp.meta.nvar
 Jx = vcat(jac_nl(nlp, x),jacG(nlp, x),jacH(nlp, x))
 vals[1 : n*m] .= Jx[:]
 return vals
end
"""
    vals = jac_coord(nlp, x)
Evaluate ``∇c(x)``, the constraint's Jacobian at `x` in sparse coordinate format.
"""
function jac_coord(nlp :: AbstractMPCCModel, x :: AbstractVector)
  vals = Vector{eltype(x)}(undef, nlp.meta.nnzj+2*nlp.meta.nvar * nlp.meta.ncc)
  return jac_coord!(nlp, x, vals)
end

"""
    Jx = jac(nlp, x)
Evaluate ``[∇c(x), -∇G(x), -∇H(x)]v``, the constraint's Jacobian at `x` as a sparse matrix.
"""
function jac(nlp :: AbstractMPCCModel, x :: AbstractVector)
  rows, cols = jac_structure(nlp)
  vals = jac_coord(nlp, x)
  return sparse(rows, cols, vals, nlp.meta.ncon+2*nlp.meta.ncc, nlp.meta.nvar)
end

"""
    (rows,cols) = jac_nl_structure(nlp)
Return the structure of the constraint's Jacobian in sparse coordinate format.
"""
function jac_nl_structure(nlp :: AbstractMPCCModel)
  rows = Vector{Int}(undef, nlp.meta.nnzj)
  cols = Vector{Int}(undef, nlp.meta.nnzj)
  return jac_nl_structure!(nlp, rows, cols)
end

"""
    jac_nl_structure!(nlp, rows, cols)
Return the structure of the constraint's Jacobian in sparse coordinate format in place.
"""
jac_nl_structure!(:: AbstractMPCCModel, :: AbstractVector{<:Integer}, :: AbstractVector{<:Integer}) = throw(NotImplementedError("jac_nl_structure!"))

"""
    vals = jac_nl_coord!(nlp, x, vals)
Evaluate ``∇c(x)``, the constraint's Jacobian at `x` in sparse coordinate format,
rewriting `vals`.
"""
jac_nl_coord!(:: AbstractMPCCModel, :: AbstractVector, ::AbstractVector) = throw(NotImplementedError("jac_nl_coord!"))

"""
    vals = jac_nl_coord(nlp, x)
Evaluate ``∇c(x)``, the constraint's Jacobian at `x` in sparse coordinate format.
"""
function jac_nl_coord(nlp :: AbstractMPCCModel, x :: AbstractVector)
  vals = Vector{eltype(x)}(undef, nlp.meta.nnzj)
  return jac_nl_coord!(nlp, x, vals)
end

"""
    Jx = jac_nl(nlp, x)
Evaluate ``∇c(x)``, the constraint's Jacobian at `x` as a sparse matrix.
"""
function jac_nl(nlp :: AbstractMPCCModel, x :: AbstractVector)
  rows, cols = jac_nl_structure(nlp)
  vals = jac_nl_coord(nlp, x)
  return sparse(rows, cols, vals, nlp.meta.ncon, nlp.meta.nvar)
end

"""
    (rows,cols) = jacG_structure(nlp)
Return the structure of the constraint's Jacobian in sparse coordinate format.
"""
function jacG_structure(nlp :: AbstractMPCCModel)
  rows = Vector{Int}(undef, nlp.meta.nvar * nlp.meta.ncc)
  cols = Vector{Int}(undef, nlp.meta.nvar * nlp.meta.ncc)
  return jacG_structure!(nlp, rows, cols)
end

"""
    jacG_structure!(nlp, rows, cols)
Return the structure of the constraint's Jacobian in sparse coordinate format in place.
"""
jacG_structure!(:: AbstractMPCCModel, :: AbstractVector{<:Integer}, :: AbstractVector{<:Integer}) = throw(NotImplementedError("jacG_structure!"))

"""
    vals = jacG_coord!(nlp, x, vals)
Evaluate ``∇G(x)``, the constraint's Jacobian at `x` in sparse coordinate format,
rewriting `vals`.
"""
jacG_coord!(:: AbstractMPCCModel, :: AbstractVector, ::AbstractVector) = throw(NotImplementedError("jacG_coord!"))

"""
    vals = jacG_coord(nlp, x)
Evaluate ``∇G(x)``, the constraint's Jacobian at `x` in sparse coordinate format.
"""
function jacG_coord(nlp :: AbstractMPCCModel, x :: AbstractVector)
  vals = Vector{eltype(x)}(undef, nlp.meta.nvar * nlp.meta.ncc)
  return jacG_coord!(nlp, x, vals)
end
"""
    Jx = jacG(nlp, x)
Evaluate ``∇G(x)``, the constraint's Jacobian at `x` as a sparse matrix.
"""
function jacG(nlp :: AbstractMPCCModel, x :: AbstractVector)
  rows, cols = jacG_structure(nlp)
  vals = jacG_coord(nlp, x)
  return sparse(rows, cols, vals, nlp.meta.ncc, nlp.meta.nvar)
end

"""
    (rows,cols) = jacH_structure(nlp)
Return the structure of the constraint's Jacobian in sparse coordinate format.
"""
function jacH_structure(nlp :: AbstractMPCCModel)
  rows = Vector{Int}(undef, nlp.meta.nvar * nlp.meta.ncc)
  cols = Vector{Int}(undef, nlp.meta.nvar * nlp.meta.ncc)
  return jacH_structure!(nlp, rows, cols)
end

"""
    jacH_structure!(nlp, rows, cols)
Return the structure of the constraint's Jacobian in sparse coordinate format in place.
"""
jacH_structure!(:: AbstractMPCCModel, :: AbstractVector{<:Integer}, :: AbstractVector{<:Integer}) = throw(NotImplementedError("jacH_structure!"))

"""
    vals = jacH_coord!(nlp, x, vals)
Evaluate ``∇H(x)``, the constraint's Jacobian at `x` in sparse coordinate format,
rewriting `vals`.
"""
jacH_coord!(:: AbstractMPCCModel, :: AbstractVector, ::AbstractVector) = throw(NotImplementedError("jacH_coord!"))

"""
    vals = jacH_coord(nlp, x)
Evaluate ``∇H(x)``, the constraint's Jacobian at `x` in sparse coordinate format.
"""
function jacH_coord(nlp :: AbstractMPCCModel, x :: AbstractVector)
  vals = Vector{eltype(x)}(undef, nlp.meta.nvar * nlp.meta.ncc)
  return jacH_coord!(nlp, x, vals)
end
"""
    Jx = jacH(nlp, x)
Evaluate ``∇H(x)``, the constraint's Jacobian at `x` as a sparse matrix.
"""
function jacH(nlp :: AbstractMPCCModel, x :: AbstractVector)
  rows, cols = jacH_structure(nlp)
  vals = jacH_coord(nlp, x)
  return sparse(rows, cols, vals, nlp.meta.ncc, nlp.meta.nvar)
end

"""
    Jv = jprod(nlp, x, v)
Evaluate ``[∇c(x), -∇G(x), -∇H(x)]v``, the Jacobian-vector product at `x`.
"""
function jprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector)
  Jv = similar(v, nlp.meta.ncon+2*nlp.meta.ncc)
  return jprod!(nlp, x, v, Jv)
end

"""
    Jv = jnlprod(nlp, x, v)
Evaluate ``∇c(x)v``, the Jacobian-vector product at `x`.
"""
function jnlprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector)
  Jv = similar(v, nlp.meta.ncon)
  return jnlprod!(nlp, x, v, Jv)
end

"""
    JGv = jGprod(nlp, x, v)
Evaluate ``∇G(x)v``, the Jacobian-vector product at `x`.
"""
function jGprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector)
  Jv = similar(v, nlp.meta.ncc)
  return jGprod!(nlp, x, v, Jv)
end

"""
    JHv = jprod(nlp, x, v)
Evaluate ``∇H(x)v``, the Jacobian-vector product at `x`.
"""
function jHprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector)
  Jv = similar(v, nlp.meta.ncc)
  return jHprod!(nlp, x, v, Jv)
end

"""
    Jv = jprod!(nlp, x, v, Jv)
Evaluate ``[∇c(x), -∇G(x), -∇H(x)]v``, the Jacobian-vector product at `x` in place.
"""
function jprod!(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector)
  ncon, ncc = nlp.meta.ncon, nlp.meta.ncc
  if ncon > 0
      Jv[1:ncon] = jnlprod(nlp, x, v)
  end
  if ncc > 0
      Jv[ncon+1:ncon+ncc]       = - jGprod(nlp, x, v)
      Jv[ncon+ncc+1:ncon+2*ncc] = - jHprod(nlp, x, v)
  end
  return Jv
end

"""
 Jv = jnlprod!(nlp, x, v, Jv)
Evaluate ``∇c(x)v``, the Jacobian-vector product at `x` in place.
"""
jnlprod!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
    throw(NotImplementedError("jnlprod!"))

"""
  JGv = jGprod!(nlp, x, v, Jv)
Evaluate ``∇G(x)v``, the Jacobian-vector product at `x` in place.
"""
jGprod!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
    throw(NotImplementedError("jGprod!"))

"""
  JHv = jHprod!(nlp, x, v, Jv)
Evaluate ``∇H(x)v``, the Jacobian-vector product at `x` in place.
"""
jHprod!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
      throw(NotImplementedError("jHprod!"))

"""
    Jtv = jtprod(nlp, x, v, Jtv)
Evaluate ``∇c(x)^Tv``, the transposed-Jacobian-vector product at `x`.
"""
function jtprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector)
  Jtv = similar(x)
  return jtprod!(nlp, x, v, Jtv)
end
"""
    Jtv = jnltprod(nlp, x, v, Jtv)
Evaluate ``[∇c(x), -∇G(x), -∇H(x)]v``, the transposed-Jacobian-vector product at `x`.
"""
function jnltprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector)
  Jtv = similar(x)
  return jnltprod!(nlp, x, v, Jtv)
end
"""
    JGtv = jGtprod(nlp, x, v, Jtv)
Evaluate ``∇G(x)^Tv``, the transposed-Jacobian-vector product at `x`.
"""
function jGtprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector)
  Jtv = similar(x)
  return jGtprod!(nlp, x, v, Jtv)
end
"""
    JHtv = jtprod(nlp, x, v, Jtv)
Evaluate ``∇H(x)^Tv``, the transposed-Jacobian-vector product at `x`.
"""
function jHtprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector)
  Jtv = similar(x)
  return jHtprod!(nlp, x, v, Jtv)
end

"""
    Jtv = jtprod!(nlp, x, v, Jtv)
Evaluate ``[∇c(x), -∇G(x), -∇H(x)]^Tv``, the transposed-Jacobian-vector product at `x` in place.
"""
function jtprod!(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector)
  ncon, ncc = nlp.meta.ncon, nlp.meta.ncc
  Jv = zeros(nlp.meta.nvar)
  if ncon > 0
     Jv += jnltprod(nlp, x, v[1:ncon])
  end
  if ncc > 0
      Jv += - jGtprod(nlp, x, v[ncon+1:ncon+ncc])
      Jv += - jHtprod(nlp, x, v[ncon+ncc+1:ncon+2*ncc])
  end
  return Jv
end

"""
 Jtv = jnltprod!(nlp, x, v, Jtv)
Evaluate ``∇c(x)^Tv``, the transposed-Jacobian-vector product at `x` in place.
"""
jnltprod!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
    throw(NotImplementedError("jnltprod!"))
"""
 JGtv = jGtprod!(nlp, x, v, Jtv)
 Evaluate ``∇G(x)^Tv``, the transposed-Jacobian-vector product at `x` in place.
"""
jGtprod!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
    throw(NotImplementedError("jGtprod!"))
"""
JHtv = jHtprod!(nlp, x, v, Jtv)
Evaluate ``∇H(x)^Tv``, the transposed-Jacobian-vector product at `x` in place.
"""
jHtprod!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
      throw(NotImplementedError("jHtprod!"))

"""
    J = jac_op(nlp, x)
Return the Jacobian at `x` as a linear operator.
The resulting object may be used as if it were a matrix, e.g., `J * v` or
`J' * v`.
"""
function jac_op(nlp :: AbstractMPCCModel, x :: AbstractVector)
  prod = @closure v -> jprod(nlp, x, v)
  ctprod = @closure v -> jtprod(nlp, x, v)
  return LinearOperator{eltype(x)}(nlp.meta.ncon, nlp.meta.nvar,
                                   false, false, prod, ctprod, ctprod)
end

"""
    J = jac_op!(nlp, x, Jv, Jtv)
Return the Jacobian at `x` as a linear operator.
The resulting object may be used as if it were a matrix, e.g., `J * v` or
`J' * v`. The values `Jv` and `Jtv` are used as preallocated storage for the
operations.
"""
function jac_op!(nlp :: AbstractMPCCModel, x :: AbstractVector,
                 Jv :: AbstractVector, Jtv :: AbstractVector)
  prod = @closure v -> jprod!(nlp, x, v, Jv)
  ctprod = @closure v -> jtprod!(nlp, x, v, Jtv)
  return LinearOperator{eltype(x)}(nlp.meta.ncon, nlp.meta.nvar,
                                   false, false, prod, ctprod, ctprod)
end

"""
    J = jacG_op(nlp, x)
Return the Jacobian at `x` as a linear operator.
The resulting object may be used as if it were a matrix, e.g., `J * v` or
`J' * v`.
"""
function jacG_op(nlp :: AbstractMPCCModel, x :: AbstractVector)
  prod = @closure v -> jGprod(nlp, x, v)
  ctprod = @closure v -> jGtprod(nlp, x, v)
  return LinearOperator{eltype(x)}(nlp.meta.ncon, nlp.meta.nvar,
                                   false, false, prod, ctprod, ctprod)
end

"""
    J = jacG_op!(nlp, x, Jv, Jtv)
Return the Jacobian at `x` as a linear operator.
The resulting object may be used as if it were a matrix, e.g., `J * v` or
`J' * v`. The values `Jv` and `Jtv` are used as preallocated storage for the
operations.
"""
function jacG_op!(nlp :: AbstractMPCCModel, x :: AbstractVector,
                 Jv :: AbstractVector, Jtv :: AbstractVector)
  prod = @closure v -> jGprod!(nlp, x, v, Jv)
  ctprod = @closure v -> jGtprod!(nlp, x, v, Jtv)
  return LinearOperator{eltype(x)}(nlp.meta.ncon, nlp.meta.nvar,
                                   false, false, prod, ctprod, ctprod)
end

"""
    J = jacH_op(nlp, x)
Return the Jacobian at `x` as a linear operator.
The resulting object may be used as if it were a matrix, e.g., `J * v` or
`J' * v`.
"""
function jacH_op(nlp :: AbstractMPCCModel, x :: AbstractVector)
  prod = @closure v -> jHprod(nlp, x, v)
  ctprod = @closure v -> jHtprod(nlp, x, v)
  return LinearOperator{eltype(x)}(nlp.meta.ncon, nlp.meta.nvar,
                                   false, false, prod, ctprod, ctprod)
end

"""
    J = jacH_op!(nlp, x, Jv, Jtv)
Return the Jacobian at `x` as a linear operator.
The resulting object may be used as if it were a matrix, e.g., `J * v` or
`J' * v`. The values `Jv` and `Jtv` are used as preallocated storage for the
operations.
"""
function jacH_op!(nlp :: AbstractMPCCModel, x :: AbstractVector,
                 Jv :: AbstractVector, Jtv :: AbstractVector)
  prod   = @closure v -> jHprod!(nlp, x, v, Jv)
  ctprod = @closure v -> jHtprod!(nlp, x, v, Jtv)
  return LinearOperator{eltype(x)}(nlp.meta.ncon, nlp.meta.nvar,
                                   false, false, prod, ctprod, ctprod)
end

function jth_hprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector, j::Integer)
  hv = Vector{eltype(x)}(undef, nlp.meta.nvar)
  return jth_hprod!(nlp, x, v, j, hv)
end

jth_hprod!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::Integer, ::AbstractVector) =
  throw(NotImplementedError("jth_hprod!"))

function ghjvprod(nlp::AbstractMPCCModel, x::AbstractVector, g::AbstractVector, v::AbstractVector)
  gHv = Vector{eltype(x)}(undef, nlp.meta.ncon)
  return ghjvprod!(nlp, x, g, v, gHv)
end

ghjvprod!(::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
  throw(NotImplementedError("ghjvprod!"))

"""
    (rows,cols) = hess_structure(nlp)
Return the structure of the Lagrangian Hessian in sparse coordinate format.
"""
function hess_structure(nlp :: AbstractMPCCModel)
  rows = Vector{Int}(undef, nlp.meta.nnzh)
  cols = Vector{Int}(undef, nlp.meta.nnzh)
  return hess_structure!(nlp, rows, cols)
end

"""
    hess_structure!(nlp, rows, cols)
Return the structure of the Lagrangian Hessian in sparse coordinate format in place.
"""
hess_structure!(:: AbstractMPCCModel, ::AbstractVector{<: Integer}, ::AbstractVector{<: Integer}) = throw(NotImplementedError("hess_structure!"))

"""
    vals = hess_coord!(nlp, x, vals; obj_weight=1.0)
Evaluate the objective Hessian at `x` in sparse coordinate format,
with objective function scaled by `obj_weight`, i.e.,
OBJECTIVE_HESSIAN, rewriting `vals`.
Only the lower triangle is returned.
"""
function hess_coord!(nlp :: AbstractMPCCModel, x :: AbstractVector, vals :: AbstractVector; obj_weight :: Real=1.0)
  return hess_coord!(nlp, x, zeros(nlp.meta.ncon), vals, obj_weight=obj_weight)
end

"""
    vals = hess_coord!(nlp, x, y, vals; obj_weight=1.0)
Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN, rewriting `vals`.
Only the lower triangle is returned.
"""
hess_coord!(nlp:: AbstractMPCCModel, :: AbstractVector, :: AbstractVector, ::AbstractVector; obj_weight :: Real=1.0) = throw(NotImplementedError("hess_coord!"))

"""
    vals = hess_coord(nlp, x; obj_weight=1.0)
Evaluate the objective Hessian at `x` in sparse coordinate format,
with objective function scaled by `obj_weight`, i.e.,
OBJECTIVE_HESSIAN.
Only the lower triangle is returned.
"""
function hess_coord(nlp :: AbstractMPCCModel, x :: AbstractVector; obj_weight::Real=one(eltype(x)))
  vals = Vector{eltype(x)}(undef, nlp.meta.nnzh)
  return hess_coord!(nlp, x, vals; obj_weight=obj_weight)
end

"""
    vals = hess_coord(nlp, x, y; obj_weight=1.0)
Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN.
Only the lower triangle is returned.
"""
function hess_coord(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector; obj_weight::Real=one(eltype(x)))
  vals = Vector{eltype(x)}(undef, nlp.meta.nnzh)
  return hess_coord!(nlp, x, y, vals; obj_weight=obj_weight)
end

"""
    Hx = hess(nlp, x; obj_weight=1.0)
Evaluate the objective Hessian at `x` as a sparse matrix,
with objective function scaled by `obj_weight`, i.e.,
OBJECTIVE_HESSIAN.
Only the lower triangle is returned.
"""
function hess(nlp::AbstractMPCCModel, x::AbstractVector; obj_weight::Real=one(eltype(x)))
  rows, cols = hess_structure(nlp)
  vals = hess_coord(nlp, x, obj_weight=obj_weight)
  return tril(sparse(rows, cols, vals, nlp.meta.nvar, nlp.meta.nvar))
end

"""
    Hx = hess(nlp, x, y; obj_weight=1.0)
Evaluate the Lagrangian Hessian at `(x,y)` as a sparse matrix,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN.
Only the lower triangle is returned.
"""
function hess(nlp::AbstractMPCCModel, x::AbstractVector, y::AbstractVector; obj_weight::Real=one(eltype(x)))
  rows, cols = hess_structure(nlp)
  vals = hess_coord(nlp, x, y, obj_weight=obj_weight)
  return tril(sparse(rows, cols, vals, nlp.meta.nvar, nlp.meta.nvar))
end

"""
    (rows,cols) = hessG_structure(nlp)
Return the structure of the Lagrangian Hessian in sparse coordinate format.
"""
function hessG_structure(nlp :: AbstractMPCCModel)
  rows = Vector{Int}(undef, nlp.meta.nnzh)
  cols = Vector{Int}(undef, nlp.meta.nnzh)
  return hessG_structure!(nlp, rows, cols)
end

"""
    hessG_structure!(nlp, rows, cols)
Return the structure of the Lagrangian Hessian in sparse coordinate format in place.
"""
hessG_structure!(:: AbstractMPCCModel, ::AbstractVector{<: Integer}, ::AbstractVector{<: Integer}) = throw(NotImplementedError("hessG_structure!"))

"""
    vals = hess_coord!(nlp, x, y, vals)
Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN, rewriting `vals`.
Only the lower triangle is returned.
"""
hessG_coord!(nlp:: AbstractMPCCModel, :: AbstractVector, :: AbstractVector, ::AbstractVector) = throw(NotImplementedError("hessG_coord!"))

"""
    vals = hessG_coord(nlp, x, y)
Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN.
Only the lower triangle is returned.
"""
function hessG_coord(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector)
  vals = Vector{eltype(x)}(undef, nlp.meta.nnzh)
  return hessG_coord!(nlp, x, y, vals)
end

"""
    Hx = hessG(nlp, x, y)
Evaluate the Lagrangian Hessian at `(x,y)` as a sparse matrix,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN.
Only the lower triangle is returned.
"""
function hessG(nlp::AbstractMPCCModel, x::AbstractVector, y::AbstractVector)
  rows, cols = hessG_structure(nlp)
  vals = hessG_coord(nlp, x, y)
  return tril(sparse(rows, cols, vals, nlp.meta.nvar, nlp.meta.nvar))
end

"""
    (rows,cols) = hessG_structure(nlp)
Return the structure of the Lagrangian Hessian in sparse coordinate format.
"""
function hessH_structure(nlp :: AbstractMPCCModel)
  rows = Vector{Int}(undef, nlp.meta.nnzh)
  cols = Vector{Int}(undef, nlp.meta.nnzh)
  return hessH_structure!(nlp, rows, cols)
end

"""
    hessH_structure!(nlp, rows, cols)
Return the structure of the Lagrangian Hessian in sparse coordinate format in place.
"""
hessH_structure!(:: AbstractMPCCModel, ::AbstractVector{<: Integer}, ::AbstractVector{<: Integer}) = throw(NotImplementedError("hessH_structure!"))

"""
    vals = hessH_coord!(nlp, x, y, vals)
Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN, rewriting `vals`.
Only the lower triangle is returned.
"""
hessH_coord!(nlp:: AbstractMPCCModel, :: AbstractVector, :: AbstractVector, ::AbstractVector) = throw(NotImplementedError("hessH_coord!"))

"""
    vals = hessG_coord(nlp, x, y)
Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN.
Only the lower triangle is returned.
"""
function hessH_coord(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector)
  vals = Vector{eltype(x)}(undef, nlp.meta.nnzh)
  return hessH_coord!(nlp, x, y, vals)
end

"""
    Hx = hessH(nlp, x, y)
Evaluate the Lagrangian Hessian at `(x,y)` as a sparse matrix,
with objective function scaled by `obj_weight`, i.e.,
LAGRANGIAN_HESSIAN.
Only the lower triangle is returned.
"""
function hessH(nlp::AbstractMPCCModel, x::AbstractVector, y::AbstractVector)
  rows, cols = hessH_structure(nlp)
  vals = hessH_coord(nlp, x, y)
  return tril(sparse(rows, cols, vals, nlp.meta.nvar, nlp.meta.nvar))
end

"""
    Hv = hprod(nlp, x, v; obj_weight=1.0)
Evaluate the product of the objective Hessian at `x` with the vector `v`,
with objective function scaled by `obj_weight`, where the objective Hessian is
OBJECTIVE_HESSIAN.
"""
function hprod(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector; obj_weight::Real=one(eltype(x)))
  Hv = similar(x)
  return hprod!(nlp, x, v, Hv; obj_weight=obj_weight)
end

"""
    Hv = hprod(nlp, x, y, v; obj_weight=1.0)
Evaluate the product of the Lagrangian Hessian at `(x,y)` with the vector `v`,
with objective function scaled by `obj_weight`, where the Lagrangian Hessian is
LAGRANGIAN_HESSIAN.
"""
function hprod(nlp::AbstractMPCCModel, x::AbstractVector, y::AbstractVector, v::AbstractVector; obj_weight::Real=one(eltype(x)))
  Hv = similar(x)
  return hprod!(nlp, x, y, v, Hv; obj_weight=obj_weight)
end

"""
    Hv = hGprod(nlp, x, y, v)
Evaluate the product of the G Hessian at `(x,y)` with the vector `v`.
"""
function hGprod(nlp::AbstractMPCCModel, x::AbstractVector, y::AbstractVector, v::AbstractVector)
  Hv = similar(x)
  return hGprod!(nlp, x, y, v, Hv)
end

"""
    Hv = hHprod(nlp, x, y, v)
Evaluate the product of the G Hessian at `(x,y)` with the vector `v`.
"""
function hHprod(nlp::AbstractMPCCModel, x::AbstractVector, y::AbstractVector, v::AbstractVector)
  Hv = similar(x)
  return hHprod!(nlp, x, y, v, Hv)
end

"""
    Hv = hprod!(nlp, x, v, Hv; obj_weight=1.0)
Evaluate the product of the objective Hessian at `x` with the vector `v` in
place, with objective function scaled by `obj_weight`, where the objective Hessian is
OBJECTIVE_HESSIAN.
"""
function hprod!(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector, Hv::AbstractVector; obj_weight :: Real=1.0)
  return hprod!(nlp, x, zeros(nlp.meta.ncon), v, Hv, obj_weight=obj_weight)
end

"""
    Hv = hprod!(nlp, x, y, v, Hv; obj_weight=1.0)
Evaluate the product of the Lagrangian Hessian at `(x,y)` with the vector `v` in
place, with objective function scaled by `obj_weight`, where the Lagrangian Hessian is
LAGRANGIAN_HESSIAN.
"""
hprod!(nlp::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector; obj_weight :: Real=1.0) =
  throw(NotImplementedError("hprod!"))

"""
 Hv = hGprod!(nlp, x, v, Hv; obj_weight=1.0)
Evaluate the product of the objective Hessian at `x` with the vector `v` in
place, with objective function scaled by `obj_weight`, where the objective Hessian is
OBJECTIVE_HESSIAN.
"""
function hGprod!(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector, Hv::AbstractVector)
 return hGprod!(nlp, x, zeros(nlp.meta.ncon), v, Hv)
end

"""
Hv = hGprod!(nlp, x, y, v, Hv; obj_weight=1.0)
Evaluate the product of the Lagrangian Hessian at `(x,y)` with the vector `v` in
place, with objective function scaled by `obj_weight`, where the Lagrangian Hessian is
LAGRANGIAN_HESSIAN.
"""
hGprod!(nlp::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
throw(NotImplementedError("hGprod!"))

"""
 Hv = hHprod!(nlp, x, v, Hv; obj_weight=1.0)
Evaluate the product of the objective Hessian at `x` with the vector `v` in
place, with objective function scaled by `obj_weight`, where the objective Hessian is
OBJECTIVE_HESSIAN.
"""
function hHprod!(nlp::AbstractMPCCModel, x::AbstractVector, v::AbstractVector, Hv::AbstractVector)
 return hHprod!(nlp, x, zeros(nlp.meta.ncon), v, Hv)
end

"""
Hv = hHprod!(nlp, x, y, v, Hv; obj_weight=1.0)
Evaluate the product of the Lagrangian Hessian at `(x,y)` with the vector `v` in
place, with objective function scaled by `obj_weight`, where the Lagrangian Hessian is
LAGRANGIAN_HESSIAN.
"""
hHprod!(nlp::AbstractMPCCModel, ::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector) =
throw(NotImplementedError("hHprod!"))

"""
    H = hess_op(nlp, x; obj_weight=1.0)
Return the objective Hessian at `x` with objective function scaled by
`obj_weight` as a linear operator. The resulting object may be used as if it were a
matrix, e.g., `H * v`. The linear operator H represents
OBJECTIVE_HESSIAN.
"""
function hess_op(nlp :: AbstractMPCCModel, x :: AbstractVector; obj_weight::Real=one(eltype(x)))
  prod = @closure v -> hprod(nlp, x, v; obj_weight=obj_weight)
  return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar,
                                   true, true, prod, prod, prod)
end

"""
    H = hess_op(nlp, x, y; obj_weight=1.0)
Return the Lagrangian Hessian at `(x,y)` with objective function scaled by
`obj_weight` as a linear operator. The resulting object may be used as if it were a
matrix, e.g., `H * v`. The linear operator H represents
LAGRANGIAN_HESSIAN.
"""
function hess_op(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector; obj_weight::Real=one(eltype(x)))
  prod = @closure v -> hprod(nlp, x, y, v; obj_weight=obj_weight)
  return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar,
                                   true, true, prod, prod, prod)
end

"""
    H = hessG_op(nlp, x, y; obj_weight=1.0)
Return the Lagrangian Hessian at `(x,y)` with objective function scaled by
`obj_weight` as a linear operator. The resulting object may be used as if it were a
matrix, e.g., `H * v`. The linear operator H represents
LAGRANGIAN_HESSIAN.
"""
function hessG_op(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector)
  prod = @closure v -> hGprod(nlp, x, y, v)
  return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar,
                                   true, true, prod, prod, prod)
end

"""
    H = hessH_op(nlp, x, y; obj_weight=1.0)
Return the Lagrangian Hessian at `(x,y)` with objective function scaled by
`obj_weight` as a linear operator. The resulting object may be used as if it were a
matrix, e.g., `H * v`. The linear operator H represents
LAGRANGIAN_HESSIAN.
"""
function hessH_op(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector)
  prod = @closure v -> hHprod(nlp, x, y, v)
  return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar,
                                   true, true, prod, prod, prod)
end

"""
    H = hess_op!(nlp, x, Hv; obj_weight=1.0)
Return the objective Hessian at `x` with objective function scaled by
`obj_weight` as a linear operator, and storing the result on `Hv`. The resulting
object may be used as if it were a matrix, e.g., `w = H * v`. The vector `Hv` is
used as preallocated storage for the operation.  The linear operator H
represents
OBJECTIVE_HESSIAN.
"""
function hess_op!(nlp :: AbstractMPCCModel, x :: AbstractVector, Hv :: AbstractVector; obj_weight::Real=one(eltype(x)))
  prod = @closure v -> hprod!(nlp, x, v, Hv; obj_weight=obj_weight)
  return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar,
                                   true, true, prod, prod, prod)
end

"""
    H = hess_op!(nlp, x, y, Hv; obj_weight=1.0)
Return the Lagrangian Hessian at `(x,y)` with objective function scaled by
`obj_weight` as a linear operator, and storing the result on `Hv`. The resulting
object may be used as if it were a matrix, e.g., `w = H * v`. The vector `Hv` is
used as preallocated storage for the operation.  The linear operator H
represents
LAGRANGIAN_HESSIAN.
"""
function hess_op!(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector, Hv :: AbstractVector; obj_weight::Real=one(eltype(x)))
  prod = @closure v -> hprod!(nlp, x, y, v, Hv; obj_weight=obj_weight)
  return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar,
                                   true, true, prod, prod, prod)
end

"""
    H = hessG_op!(nlp, x, y, Hv; obj_weight=1.0)
Return the Lagrangian Hessian at `(x,y)` with objective function scaled by
`obj_weight` as a linear operator, and storing the result on `Hv`. The resulting
object may be used as if it were a matrix, e.g., `w = H * v`. The vector `Hv` is
used as preallocated storage for the operation.  The linear operator H
represents
LAGRANGIAN_HESSIAN.
"""
function hessG_op!(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector, Hv :: AbstractVector)
  prod = @closure v -> hGprod!(nlp, x, y, v, Hv)
  return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar,
                                   true, true, prod, prod, prod)
end

"""
    H = hessH_op!(nlp, x, y, Hv; obj_weight=1.0)
Return the Lagrangian Hessian at `(x,y)` with objective function scaled by
`obj_weight` as a linear operator, and storing the result on `Hv`. The resulting
object may be used as if it were a matrix, e.g., `w = H * v`. The vector `Hv` is
used as preallocated storage for the operation.  The linear operator H
represents
LAGRANGIAN_HESSIAN.
"""
function hessH_op!(nlp :: AbstractMPCCModel, x :: AbstractVector, y :: AbstractVector, Hv :: AbstractVector)
  prod = @closure v -> hHprod!(nlp, x, y, v, Hv)
  return LinearOperator{eltype(x)}(nlp.meta.nvar, nlp.meta.nvar,
                                   true, true, prod, prod, prod)
end

ext: jl


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
SUITE["Platooning"] = BenchmarkGroup()

include("Platooning.jl")

PLAD01_BND42, options_PLAD01_BND42 =
    platooning(; deterministic_switching=true, time_horizon=20.,
                 allowed_distance=42.)
PLAD01_BND30, options_PLAD01_BND30 =
    platooning(; deterministic_switching=true, time_horizon=20.,
                 allowed_distance=30.)
PLAN01_UNB50, options_PLAN01_UNB50 =
    platooning(; deterministic_switching=true, time_horizon=Inf,
                 allowed_distance=50.)

# ==============================================================================
# Decomposition-based approach
# ==============================================================================

# algorithm-specific options
options_PLAD01_BND42[:mode] = "check"
options_PLAD01_BND30[:mode] = "check"
options_PLAN01_UNB50[:mode] = "check"

𝑂_common = Options(:partition => [1:10])
𝑂_dense_options_PLAD01_BND42 = merge(𝑂_common, Options(:δ => 0.01))
𝑂_dense_options_PLAD01_BND30 = merge(𝑂_common, Options(:δ => 0.00001))
𝑂_dense_options_PLAN01_UNB50 = merge(𝑂_common, Options(:δ => 0.03))
𝑂_discrete = merge(𝑂_common, Options(:discretization => "nobloating", :δ => 0.1))

opC_dense_PLAD01_BND42 = BFFPSV18(𝑂_dense_options_PLAD01_BND42)
opC_dense_PLAD01_BND30 = BFFPSV18(𝑂_dense_options_PLAD01_BND30)
opC_dense_PLAN01_UNB50 = BFFPSV18(𝑂_dense_options_PLAN01_UNB50)
opC_discrete = BFFPSV18(𝑂_discrete)
opD = LazyDiscretePost(:lazy_R⋂I => true, :lazy_R⋂G => false)

# single run to verify that specification holds
res = solve(PLAD01_BND42, options_PLAD01_BND42, opC_dense_PLAD01_BND42, opD)
@assert res.satisfied
res = solve(PLAD01_BND42, options_PLAD01_BND42, opC_discrete, opD)
@assert res.satisfied
# res = solve(PLAD01_BND30, options_PLAD01_BND30, opC_dense_PLAD01_BND30, opD)
# @assert res.satisfied
# res = solve(PLAD01_BND30, options_PLAD01_BND30, opC_discrete, opD)
# @assert res.satisfied
res = solve(PLAN01_UNB50, options_PLAN01_UNB50, opC_dense_PLAN01_UNB50, opD)
@assert res.satisfied
res = solve(PLAN01_UNB50, options_PLAN01_UNB50, opC_discrete, opD)
@assert res.satisfied

# benchmark
SUITE["Platooning"]["PLAD01_BND42", "dense"] =
    @benchmarkable solve($PLAD01_BND42, $options_PLAD01_BND42, $opC_dense_PLAD01_BND42, $opD)
SUITE["Platooning"]["PLAD01_BND42", "discrete"] =
    @benchmarkable solve($PLAD01_BND42, $options_PLAD01_BND42, $opC_discrete, $opD)
# SUITE["Platooning"]["PLAD01_BND30", "dense"] =
#     @benchmarkable solve($PLAD01_BND30, $options_PLAD01_BND30, $opC_dense_PLAD01_BND30, $opD)
# SUITE["Platooning"]["PLAD01_BND30", "discrete"] =
#     @benchmarkable solve($PLAD01_BND30, $options_PLAD01_BND30, $opC_discrete, $opD)
SUITE["Platooning"]["PLAN01_UNB50", "dense"] =
    @benchmarkable solve($PLAN01_UNB50, $options_PLAN01_UNB50, $opC_dense_PLAN01_UNB50, $opD)
SUITE["Platooning"]["PLAN01_UNB50", "discrete"] =
    @benchmarkable solve($PLAN01_UNB50, $options_PLAN01_UNB50, $opC_discrete, $opD)

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=false)

# return the sample with the smallest time value in each test
println("minimum: ", minimum(results))

# return the median for each test
println("median: ", median(results))

# ==============================================================================
# Create plots
# ==============================================================================

function plot_dashed_line!(v, time_horizon=20.)
    plot!(x->x, x->v, 0.0, time_horizon, line=2, color="red", linestyle=:dash,
          legend=nothing)
end

options_PLAD01_BND42[:mode] = "reach"
options_PLAD01_BND42[:plot_vars] = [0, 1]
options_PLAD01_BND42[:project_reachset] = true
res = solve(PLAD01_BND42, options_PLAD01_BND42, opC_dense_PLAD01_BND42, opD)
plot(res,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{-0.5mm}{\textcolor{white}{.}}",
     ylab=L"x_1\raisebox{2mm}{\textcolor{white}{.}}",
     xtick=[0., 10., 20.], ytick=[-40., -20., 0., 20.],
     bottom_margin=6mm, left_margin=2mm,
     size=(1000, 1000))
plot_dashed_line!(-42.)
savefig(@relpath "PLAD01_BND42.png")

options_PLAN01_UNB50[:mode] = "reach"
options_PLAN01_UNB50[:plot_vars] = [0, 1]
options_PLAN01_UNB50[:project_reachset] = true
res = solve(PLAN01_UNB50, options_PLAN01_UNB50, opC_dense_PLAN01_UNB50, opD)
time_horizon = res.Xk[end].t_start
plot(res,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{-0.5mm}{\textcolor{white}{.}}",
     ylab=L"x_1\raisebox{2mm}{\textcolor{white}{.}}",
     xtick=[0., 100., 200.], ytick=[-40., -20., 0., 20.],
     bottom_margin=6mm, left_margin=2mm,
     size=(1000, 1000))
plot_dashed_line!(-50., time_horizon)
savefig(@relpath "PLAN01_UNB50.png")

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
module SpSkipGramD
# share sparse representations

using CUDArt: HostArray, CudaArray, synchronize, null_stream
using JLD: jldopen
using CycUtils: invertable_unique, get, set!
using Veclib: gamma_kernel!, gamma_kernel_loss!, normalize_columns!

export SparseSkipGramD, train, init, close_device

include("config.jl")
include("array_handler.jl")
include("samples.jl")
include("utils.jl")

type SparseSkipGramD{Tensor} # Tensor could be HostArray{Float32, 2}
    α::Tensor
    β::Tensor
    A::Tensor

    # the following fields are training parameters.
    m::Int
    k::Int
    lr::Float32
    lr0::Float32
    λ::Float32
    ɛ::Float32

    # maps
    idx2word::Array{String, 1}
    word2idx::Dict{String, UInt32}

    # the following fields are named memory for training
    ∇αβ::Tensor
    ᾱ::Tensor
    β̄::Tensor
    a⃗::Tensor
    b⃗::Tensor
    γ::Tensor
    σ::Tensor
    ∇α::Tensor
    M::Tensor
    ∇β::Tensor
    ∇A::Tensor
    idxes::Tensor
    sign::Tensor

    function SparseSkipGramD(
        num_base::Integer,     # number of atom words
        num_vocab::Integer,    # number of vocabulary
        mini_batch::Integer,     # mini-batch size
        negative::Integer,        # number of negative samples per positive sample
        lr::AbstractFloat=0.01,       # learning rate
        λ::AbstractFloat=0.01,         # sparseness weight
        ɛ::AbstractFloat=0.01,
        init_model::String=""   # dumped init_model.
    )
        @assert num_base ≤ num_vocab "The number of atom words must be less than the vocabulary size"

        α_ = randn(Float32, num_base, num_vocab) ./ sqrt(num_base)

        i2w = Array{String, 1}()
        w2i = Dict{String, UInt32}()

        if init_model != ""
            e, c, i2w = jldopen(init_model, "r") do f
                (read(f, "embeddings"), read(f, "contexts"),
                 read(f, "idx2word"))
            end
            @assert num_vocab ≤ length(i2w) "Initial model must contains as many vocabularies as the model to be trained."

            d = size(e, 1)
            B_ = randn(Float32, num_base, num_base) ./ sqrt(num_base)
            C_ = randn(Float32, num_base, num_base) ./ sqrt(num_base)

            i2w = i2w[1:num_vocab]
            for (i, w) in enumerate(i2w)
                w2i[w] = i
            end

            α_[1:num_base, 1:num_base] = eye(num_base)
            copy!(B_, @view e[:, 1:num_base])
            copy!(C_, @view c[:, 1:num_base])
            A_ = B_' * C_
        else
            A_ = abs(randn(Float32, (num_base, num_base)) ./ sqrt(num_base))
        end

        α = Tensor(α_)
        β = α
        A = Tensor(A_)

        ∇αβ = Tensor(Float32, num_base*(negative+2)*mini_batch)
        ∇α, ∇β = splitarray(∇αβ, (num_base, mini_batch),
                                   (num_base, negative+1, mini_batch))

        ᾱ = Tensor(Float32, num_base, mini_batch)
        β̄ = Tensor(Float32, num_base, negative+1, mini_batch)
        a⃗ = Tensor(Float32, num_base, mini_batch)
        b⃗ = Tensor(Float32, num_base, negative+1, mini_batch)
        γ = Tensor(Float32, negative+1, mini_batch)
        σ = Tensor(Float32, negative+1, mini_batch)
        M = Tensor(Float32, num_base, mini_batch)
        ∇A = Tensor(Float32, num_base, num_base)

        idxes = Tensor(UInt32, (negative+2)*mini_batch)

        fill!(∇αβ, 0.0)
        fill!(ᾱ, 0.0)
        fill!(β̄, 0.0)
        fill!(a⃗, 0.0)
        fill!(b⃗, 0.0)
        fill!(γ, 0.0)
        fill!(σ, 0.0)
        fill!(M, 0.0)
        fill!(∇A, 0.0)
        fill!(idxes, 0)

        tmp = -ones(Float32, (negative+1, mini_batch))
        tmp[1, :] = 1.0f0
        s = Tensor(tmp)

        new{Tensor}(α, β, A,                          # model parameters
            mini_batch, negative, lr, lr, λ, ɛ,       # training parameters
            i2w, w2i,                                 # maps
            ∇αβ, ᾱ, β̄, a⃗, b⃗, γ, σ, ∇α, M, ∇β, ∇A, idxes, s
        )
    end
end

function save(model::SparseSkipGramD{Array}, fname::String; save_maps=false)
    jldopen(fname, "w") do f
        write(f, "α", model.α)
        write(f, "A", model.A)

        # the following fields are training parameters.
        write(f, "m", model.m)
        write(f, "k", model.k)
        write(f, "lr", model.lr0)
        write(f, "λ", model.λ)
        write(f, "ɛ", model.ɛ)

        # maps
        if save_maps == true
            write(f, "idx2word", model.idx2word)
            write(f, "word2idx", model.word2idx)
        end
    end
end

num_vocab(model::SparseSkipGramD) = size(model.α, 2)
num_base(model::SparseSkipGramD) = size(model.α, 1)
# num_vecdim(model::SparseSkipGramD) = size(model.B, 1)
num_minibatch(model::SparseSkipGramD) = size(model.a⃗, 2)
num_negative(model::SparseSkipGramD) = size(model.b⃗, 2) - 1

function to_device{M <: Union{SparseSkipGramD{Array},
                    SparseSkipGramD{HostArray}}}(model::M)
    num_base_ = num_base(model)
    num_vocab_ = num_vocab(model)
    mini_batch = num_minibatch(model)
    negative = num_negative(model)

    d_model = SparseSkipGramD{CudaArray}(num_base_, num_vocab_, mini_batch,
                                        negative, model.lr, model.λ, model.ɛ)
    for fn in fieldnames(typeof(d_model))
        if fn == :β
            continue
        end
        dest = getfield(d_model, fn)
        if isa(dest, CudaArray)
             src = getfield(model, fn)
             copy!(dest, src)
        end
    end
    return d_model
end

@generated function get_training_model(model::SparseSkipGramD)
    if BACKEND == "GPU"
        return :(to_device(model))
    else
        return :(model)
    end
end

function compute_grads_and_updateA!(model, batch)
    const m = model.m
    const k = model.k
    const A = model.A
    const M = model.M
    w, c = batch
    const nb = num_base(model)
    const ᾱ = takelast!(model.ᾱ, model.α, w)
    const β̄ = takelast!(model.β̄, model.β, c)
    const a⃗ = model.a⃗
    const b⃗ = model.b⃗
    const γ = model.γ
    const ɛ = model.ɛ
    const ∇α = model.∇α
    const ∇β = model.∇β
    const ∇A = model.∇A

    flat_b⃗ = reshape(b⃗, (nb, (k+1)*m))
    gemm!('T', 'N', 1.0f0, A, ᾱ, 0.0f0, a⃗)
    gemm!('N', 'N', 1.0f0, A, reshape(β̄, (nb, (k+1)*m)), 0.0f0, flat_b⃗)

    α_ = reshape(ᾱ, (nb, 1, m))
    γ_ = reshape(γ, (k+1, 1, m))
    gemm_strided_batched!('T', 'N', 1.0f0, b⃗, α_, 0.0f0, γ_)
    gamma_kernel_loss!(γ, model.σ, model.sign)

    #∇α
    ∇α_ = reshape(∇α, (nb, 1, m))
    gemm_strided_batched!('N', 'N', 1.0f0, b⃗, γ_, 0.0f0, ∇α_)

    # ∇β
    a⃗_ = reshape(a⃗, (nb, 1, m))
    gemm_strided_batched!('N', 'T', 1.0f0, a⃗_, γ_, 0.0f0, ∇β)

    # ∇A
    M_ = reshape(M, (nb, 1, m))
    gemm_strided_batched!('N', 'N', 1.0f0, β̄, γ_, 0.0f0, M_)
    # gemm!('N', 'T', 1.0f0, ᾱ, M, 0.0f0, ∇A)
    # axpy!(ɛ, vec(A), vec(∇A))

    # update A
    gemm!('N', 'T', -model.lr, ᾱ, M, 1.0f0-model.lr*model.ɛ, A)
    nothing
end

function fetch_grad_sparse!(grads, model)
    copy!(grads, model.∇αβ)
    return grads
end

function gpu_sparse_update!(arr, grad, idxes, lr, lambda)
    # make sure the transferring completes.
    CUDArt.synchronize(CUDArt.null_stream)
    ccall((:sparse_update, vec_lib), Void, (Ptr{Float32}, Ptr{Float32},
        Cint, Cint, Ptr{Cuint}, Float32, Float32, Ptr{Void}, Cint),
        arr.ptr, grad.ptr, size(grad, 1), size(grad, 2),
        idxes.ptr, lr, lambda, null_stream.inner.handle, 256)
end

# function update_A!(A, gradA, lr, e, num_thread::Integer=256,
#     stream::Ptr{Void}=null_stream.inner.handle)
#
#     ccall((:update_A_neg, vec_lib), Void,
#           (Ptr{Cfloat}, Ptr{Cfloat}, Cfloat, Cfloat, Cint,
#            Ptr{Void}, Cint),
#            vec(A), vec(gradA), lr, e, length(A), stream, num_thread)
# end

@generated function proximal_update!(model, ∇αβ, idxes)
    if BACKEND == "GPU"
        quote
            const lr_::Float32 = model.lr
            dev = Int(DEVICE)
            siz_αβ = (num_base(model), length(idxes))
            d_∇αβ = CudaArray(model.∇αβ.ptr, siz_αβ, dev)
            copy!(d_∇αβ, ∇αβ)

            ptr_idxes = model.idxes.ptr
            d_idxes = CudaArray(ptr_idxes, (length(idxes),), dev)
            copy!(d_idxes, idxes)
            gpu_sparse_update!(model.α, d_∇αβ, d_idxes, lr_, model.λ)
        end
    else
        quote
            error("Not implemented")
            # lr_::Float32 = model.lr
            # cpu_sparse_update!(model.α, ∇αβ, idxes, lr_, model.λ)
        end
    end
end

function update_sparse_u1!(model, ∇αβ, grad_αβ, uniq_batch)
    idx, inv_idxes = uniq_batch

    ∇αβ = @view ∇αβ[:, 1:length(idx)]
    fill!(∇αβ, 0.0f0)

    # make sure the transferring completes.
    CUDArt.synchronize(CUDArt.null_stream)

    for (i, inv_idx) in enumerate(inv_idxes)
        ∇αβ[:, inv_idx] += grad_αβ[:, i]
    end

    proximal_update!(model, ∇αβ, idx)
    nothing
end

function update_sparse_u2!(model, ∇αβ, grad_αβ, uniq_batch)
    idx, inv_idxes = uniq_batch

    ∇αβ = @view ∇αβ[:, 1:length(idx)]
    fill!(∇αβ, 0.0f0)

    # make sure the transferring completes.
    CUDArt.synchronize(CUDArt.null_stream)

    for (i, inv_idx) in enumerate(inv_idxes)
        ∇αβ[:, inv_idx] += grad_αβ[:, i]
    end
    normalize_columns!(∇αβ)

    proximal_update!(model, ∇αβ, idx)
    nothing
end

@generated function sync_params!(cpu_model, model)
    if BACKEND == "CPU"
        return :(cpu_model)
    else
        quote
            for fn in (:α, :A)  #fieldnames(typeof(model))
                src = getfield(model, fn)
                dest = getfield(cpu_model, fn)
                copy!(dest, src)
            end
            CUDArt.synchronize(CUDArt.null_stream)
            cpu_model
        end
    end
end

const _start_time_loginfo = Ref{Float64}(0.0)
const _next_report_time = Ref{Float64}(0.0)
const _sigmas_loginfo = Ref{HostArray{Float32, 2}}()

function log_progress(model, current, progress=progress)
    global _start_time_loginfo
    global _next_report_time
    global _sigmas_loginfo
    global sampleFactory
    factory = get(sampleFactory)

    if current >= _next_report_time[]
        σ = _sigmas_loginfo[]
        copy!(σ, model.σ)
        CUDArt.synchronize(CUDArt.null_stream)
        loss = 0.0f0
        for x in σ
            loss += -log(x+eps(typeof(x)))
        end
        ppl = exp(loss/(model.m*(model.k+1)))

        _next_report_time[] += 3.0
        st = _start_time_loginfo[]
        factory = get(sampleFactory)
        prog = progress(factory) * 100.0
        elapsed = current - st
        etc = (100.0 - prog) / prog * elapsed
        hours = div(etc, 3600)
        rem = etc-hours*3600
        mins = div(rem, 60)
        secs = rem - mins * 60

        mean_speed = factory.consumed[]/elapsed
        info(@sprintf("ETA: %02d:%02d:%02d - loss: %.4f - ppl: %.3f - progress: %5.2f%% - speed: %.1f %s/s - lr: %g", hours, mins, secs, loss, ppl, prog, mean_speed, factory.corpus_type == String ? "lines" : "bytes", model.lr))
    end
    nothing
end

function init_logger(model, start_time)
    global _start_time_loginfo
    global _next_report_time
    global _sigmas_loginfo
    _start_time_loginfo[] = start_time
    _next_report_time[] = start_time
    _sigmas_loginfo[] = HostArray(Float32, (model.k+1), model.m)
end

function train(model; min_lr=1f-6, normalize_sp=false,
               saveratio=0.05, part=1.0, save_basename="D")

    global sampleFactory
    cpu_model = model
    model = to_device(cpu_model)
    beta = model.lr0 / min_lr - 1.0f0

    info(@sprintf("Train with lr=%g, normalize_sp=%s. Save after every %.2f%% of training. Use %.2f%% of the total training data. Save the trained model to \"%s*.jld\"", model.lr, string(normalize_sp), saveratio*100, part*100, save_basename))

    grads = HostArray(Float32, (num_base(model), (model.k+2)*model.m))
    grad_cache = HostArray(Float32, (num_base(model), (model.k+2)*model.m))

    factory = get(sampleFactory)
    const producer = factory.producer
    const progress_ = factory -> progress(factory)/part
    batch_set = Task(() -> producer(factory, model.m, model.k))

    const update_sparse! = normalize_sp ? update_sparse_u2! : update_sparse_u1!

    st = time()
    init_logger(model, st)

    old_batch = consume(batch_set)
    next_save_point = saveratio
    compute_grads_and_updateA!(model, old_batch)
    for batch in batch_set
        uniq_batch = invertable_unique(old_batch[1], old_batch[2])
        # synchronize(CUDArt.null_stream)
        # TODO: maybe another stream to handle the data transferring
        fetch_grad_sparse!(grads, model)
        update_sparse!(model, grad_cache, grads, uniq_batch)
        old_batch = batch
        compute_grads_and_updateA!(model, old_batch)
        prog = progress_(factory)
        if prog >= next_save_point
            if prog >= 1.0
                break
            end
            sync_params!(cpu_model, model)
            sn = joinpath(project_root,
                @sprintf("%s-checkpoint-%03.0f.jld", save_basename, prog*100.0f0))
            info("Save point reached. saveing to $sn")
            save(cpu_model, sn)
            next_save_point += saveratio
        end
        log_progress(model, time(), progress_)
        model.lr = model.lr0 / (1.0f0 + beta*prog)
    end

    sync_params!(cpu_model, model)
    sn = joinpath(project_root,
        @sprintf("%s-checkpoint-100.jld", save_basename))
    info("Training finished. Saving to $sn")
    save(cpu_model, sn, save_maps=true)
end

end # end of module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Kokkos_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Kokkos")
JLLWrappers.@generate_main_file("Kokkos", UUID("c1216c3d-6bb3-5a2b-bbbf-529b35eba709"))
end  # module Kokkos_jll


text:
using BeastUtils.RTrees

newick = "(((A:0.5,B:0.2):0.1,C:0.3):1.0,(D:0.3,E:0.4):0.7);"
tree = RTrees.parse_newick(newick)

taxa1 = ["A", "B", "C", "D", "E"]

Ψ1 = RTrees.tree_variance(tree, taxa1)

Ψ1_true = [
    1.6 1.1 1.0 0.0 0.0;
    1.1 1.3 1.0 0.0 0.0;
    1.0 1.0 1.3 0.0 0.0;
    0.0 0.0 0.0 1.0 0.7;
    0.0 0.0 0.0 0.7 1.1
    ]

tol = 1e-14

diff1 = maximum(abs.(Ψ1 - Ψ1_true))

@test diff1 < tol

perm = [5, 2, 3, 1, 4]
taxa2 = taxa1[perm]
Ψ2_true = Ψ1_true[perm, perm]
Ψ2 = RTrees.tree_variance(tree, taxa2)

diff2 = maximum(abs.(Ψ2 - Ψ2_true))

@test diff2 < tol


text:
is_call(x) = false
is_call(x::Expr) = x.head == :call

isunless(x::Expr) = x.head == :gotoifnot
isunless(x) = false
islabelnode(x::LabelNode) = true
islabelnode(x) = false
isgoto(x::GotoNode) = true
isgoto(x) = false

function is_unless_label(label, hist, histpos = 1)
    islabelnode(label) || return false
    unless = hist[histpos][1]
    unless_label = unless.args[2]
    unless_label == label.label
end
function is_goto_label(label, hist, histpos)
    islabelnode(label) || return false
    goto = hist[histpos][1]
    goto.label == label.label
end
function is_goto(goto, hist, histpos)
    isgoto(goto) || return false
    label = hist[histpos][1]
    goto.label == label.label
end
function is_unless_goto(goto, hist, histpos)
    isgoto(goto) || return false
    unless = hist[histpos][1]
    goto.label == unless.args[2]
end

save_resolve_func(x) = nothing
save_resolve_func(f::AllFuncs) = f
save_resolve_func{T}(::Type{T}) = T
save_resolve_func(f::Union{GlobalRef, Symbol}) = eval(f)

function iscolon(expr)
    isa(expr, Expr) && expr.head == :(=) || return false
    lhs, rhs = expr.args
    isa(rhs, Expr) && rhs.head == :call || return false
    save_resolve_func(rhs.args[1]) == colon || return false
    return true
end

function isstart(expr, hist)
    isa(expr, Expr) && expr.head == :(=) || return false
    lhs, rhs = expr.args
    isa(rhs, Expr) && rhs.head == :call || return false
    save_resolve_func(rhs.args[1]) == start || return false
    colon = hist[1][1]
    slot = colon.args[1]
    slotstart = rhs.args[2]
    slot == slotstart || return false
    return true
end
function isnext(expr, hist)
    isa(expr, Expr) && expr.head == :(=) || return false
    lhs, rhs = expr.args
    isa(rhs, Expr) && rhs.head == :call || return false
    save_resolve_func(rhs.args[1]) == next || return false
    colon = hist[1][1]; slot = colon.args[1]
    slotnext = rhs.args[2]
    slot == slotnext || return false
    return true
end
function is_done_unless(expr, hist)
    Sugar.isunless(expr) || return false
    condition = expr.args[1]
    isa(condition, Expr) && condition.head == :call || return false
    save_resolve_func(condition.args[1]) == (!) || return false
    done_expr = condition.args[2]
    isa(done_expr, Expr) && done_expr.head == :call || return false
    save_resolve_func(done_expr.args[1]) == done || return false
    colon = hist[1][1]; slot = colon.args[1]
    slotdone = done_expr.args[2]
    slot == slotdone || return false
    return true
end
function isgetfield(expr, slotnext)
    isa(expr, Expr) && expr.head == :(=) || return false
    lhs, rhs = expr.args
    isa(rhs, Expr) && rhs.head == :call || return false
    save_resolve_func(rhs.args[1]) == getfield || return false
    slot = rhs.args[2]
    slot == slotnext || return false
    return true
end

##########################
# code lowered AST patterns
const for_pattern = (
    iscolon,
    isstart,
    islabelnode, # loop goto label
    is_done_unless, # while condition branch
    anything, # new varnodes etc
    isnext,
    anything, # body
    Greed(islabelnode, 0:1), # optional continue label
    (l,h)-> is_goto(l, h, 3), # goto label, matching first label
    (l,h)-> is_unless_label(l, h, 4) # goto and break
)

const ifelse_pattern = (
    isunless, # if branch
    anything, # if body
    isgoto, # goto to jump over else
    (l, h)-> is_unless_label(l, h, 1),
    anything,   # else body
    (l, h)-> is_goto_label(l, h, 3) # label matching above goto
)
const while_pattern = (
    islabelnode, # loop goto label
    isunless, # while condition branch
    anything, # body
    Greed(islabelnode, 0:1), # optional continue label
    (l, h)-> is_goto(l, h, 1), # goto label, matching first label
    (l, h)-> is_unless_label(l, h, 2) # goto and break
)

const if_pattern = (
    isunless, # condition
    anything, # body
    Greed(is_unless_label, 0:1)
)

const goto_neighbours = (
    isgoto, # goto label directly next to it
    (l, h)-> is_goto_label(l, h, 1)
)

##########################
# LLVM string patterns

# Function call
# %15 = call i8 @julia_mapreduce_sc_impl_60977(%jl_value_t* %0)
llvm_function_call = r"(.*) = call (.*) @(.*)\(.*\)"

##########################
# Assembly string patterns
assembly_source_line = r"Source line: ([1-9][0-9]+)"


text:
mn=rand(3:7)
tsname="InplaceOps"
@testset "$tsname" begin

verbose && println("$tsname - Regular ($mn,$mn)")
@testset "Regular" begin
    DDT=ComplexF64
    RDT=ComplexF64
    A=joMatrix(rand(DDT,mn,mn),DDT=DDT,RDT=RDT)

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=A*x
    my=A*mx
    @! Y=A*x
    @test Y≈y
    @! MY=A*mx
    @test MY≈my

    y=conj(A)*x
    my=conj(A)*mx
    @! Y=conj(A)*x
    @test Y≈y
    @! MY=conj(A)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=transpose(A)*x
    my=transpose(A)*mx
    @! Y=transpose(A)*x
    @test Y≈y
    @! MY=transpose(A)*mx
    @test MY≈my

    y=adjoint(A)*x
    my=adjoint(A)*mx
    @! Y=adjoint(A)*x
    @test Y≈y
    @! MY=adjoint(A)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=A\x
    my=A\mx
    @! Y=A\x
    @test Y≈y
    @! MY=A\mx
    @test MY≈my

    y=conj(A)\x
    my=conj(A)\mx
    @! Y=conj(A)\x
    @test Y≈y
    @! MY=conj(A)\mx
    @test MY≈my

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=transpose(A)\x
    my=transpose(A)\mx
    @! Y=transpose(A)\x
    @test Y≈y
    @! MY=transpose(A)\mx
    @test MY≈my

    y=adjoint(A)\x
    my=adjoint(A)\mx
    @! Y=adjoint(A)\x
    @test Y≈y
    @! MY=adjoint(A)\mx
    @test MY≈my
end

verbose && println("$tsname - joMatrixInplace ($mn,$mn)")
@testset "joMatrixInplace" begin
    DDT=ComplexF64
    RDT=ComplexF64
    a=rand(DDT,mn,mn)
    A=joMatrix(a,DDT=DDT,RDT=RDT)
    AI=joMatrixInplace(a,DDT=DDT,RDT=RDT)

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=A*x
    my=A*mx
    @! Y=AI*x
    @test Y≈y
    @! MY=AI*mx
    @test MY≈my

    y=conj(A)*x
    my=conj(A)*mx
    @! Y=conj(AI)*x
    @test Y≈y
    @! MY=conj(AI)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=transpose(A)*x
    my=transpose(A)*mx
    @! Y=transpose(AI)*x
    @test Y≈y
    @! MY=transpose(AI)*mx
    @test MY≈my

    y=adjoint(A)*x
    my=adjoint(A)*mx
    @! Y=adjoint(AI)*x
    @test Y≈y
    @! MY=adjoint(AI)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=A\x
    my=A\mx
    @! Y=AI\x
    @test Y≈y
    @! MY=AI\mx
    @test MY≈my

    y=conj(A)\x
    my=conj(A)\mx
    @! Y=conj(AI)\x
    @test Y≈y
    @! MY=conj(AI)\mx
    @test MY≈my

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=transpose(A)\x
    my=transpose(A)\mx
    @! Y=transpose(AI)\x
    @test Y≈y
    @! MY=transpose(AI)\mx
    @test MY≈my

    y=adjoint(A)\x
    my=adjoint(A)\mx
    @! Y=adjoint(AI)\x
    @test Y≈y
    @! MY=adjoint(AI)\mx
    @test MY≈my
end

verbose && println("$tsname - joLooseMatrixInplace ($mn,$mn)")
@testset "joLooseMatrixInplace" begin
    DDT=ComplexF64
    RDT=ComplexF64
    a=rand(DDT,mn,mn)
    A=joMatrix(a,DDT=DDT,RDT=RDT)
    AI=joLooseMatrixInplace(a,DDT=DDT,RDT=RDT)

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=A*x
    my=A*mx
    @! Y=AI*x
    @test Y≈y
    @! MY=AI*mx
    @test MY≈my

    y=conj(A)*x
    my=conj(A)*mx
    @! Y=conj(AI)*x
    @test Y≈y
    @! MY=conj(AI)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=transpose(A)*x
    my=transpose(A)*mx
    @! Y=transpose(AI)*x
    @test Y≈y
    @! MY=transpose(AI)*mx
    @test MY≈my

    y=adjoint(A)*x
    my=adjoint(A)*mx
    @! Y=adjoint(AI)*x
    @test Y≈y
    @! MY=adjoint(AI)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=A\x
    my=A\mx
    @! Y=AI\x
    @test Y≈y
    @! MY=AI\mx
    @test MY≈my

    y=conj(A)\x
    my=conj(A)\mx
    @! Y=conj(AI)\x
    @test Y≈y
    @! MY=conj(AI)\mx
    @test MY≈my

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=transpose(A)\x
    my=transpose(A)\mx
    @! Y=transpose(AI)\x
    @test Y≈y
    @! MY=transpose(AI)\mx
    @test MY≈my

    y=adjoint(A)\x
    my=adjoint(A)\mx
    @! Y=adjoint(AI)\x
    @test Y≈y
    @! MY=adjoint(AI)\mx
    @test MY≈my
end

verbose && println("$tsname - joLinearFunctionInplace ($mn,$mn)")
@testset "joLinearFunctionInplace" begin
    DDT=ComplexF64
    RDT=ComplexF64
    a=rand(DDT,mn,mn)
    A=joMatrix(a,DDT=DDT,RDT=RDT)
    fa=factorize(a)
    AI=joLinearFunctionInplaceAll(mn,mn,
        (y1,x1)->mul!(y1,a,x1),
        (y2,x2)->mul!(y2,transpose(a),x2),
        (y3,x3)->mul!(y3,adjoint(a),x3),
        (y4,x4)->mul!(y4,conj(a),x4),
        (y5,x5)->ldiv!(y5,fa,x5),
        (y6,x6)->ldiv!(y6,transpose(fa),x6),
        (y7,x7)->ldiv!(y7,adjoint(fa),x7),
        (y8,x8)->ldiv!(y8,factorize(conj(a)),x8),
        DDT,RDT)

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=A*x
    my=A*mx
    @! Y=AI*x
    @test Y≈y
    @! MY=AI*mx
    @test MY≈my

    y=conj(A)*x
    my=conj(A)*mx
    @! Y=conj(AI)*x
    @test Y≈y
    @! MY=conj(AI)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=transpose(A)*x
    my=transpose(A)*mx
    @! Y=transpose(AI)*x
    @test Y≈y
    @! MY=transpose(AI)*mx
    @test MY≈my

    y=adjoint(A)*x
    my=adjoint(A)*mx
    @! Y=adjoint(AI)*x
    @test Y≈y
    @! MY=adjoint(AI)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=A\x
    my=A\mx
    @! Y=AI\x
    @test Y≈y
    @! MY=AI\mx
    @test MY≈my

    y=conj(A)\x
    my=conj(A)\mx
    @! Y=conj(AI)\x
    @test Y≈y
    @! MY=conj(AI)\mx
    @test MY≈my

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=transpose(A)\x
    my=transpose(A)\mx
    @! Y=transpose(AI)\x
    @test Y≈y
    @! MY=transpose(AI)\mx
    @test MY≈my

    y=adjoint(A)\x
    my=adjoint(A)\mx
    @! Y=adjoint(AI)\x
    @test Y≈y
    @! MY=adjoint(AI)\mx
    @test MY≈my
end

verbose && println("$tsname - joLooseLinearFunctionInplace ($mn,$mn)")
@testset "joLooseLinearFunctionInplace" begin
    DDT=ComplexF64
    RDT=ComplexF64
    a=rand(DDT,mn,mn)
    A=joMatrix(a,DDT=DDT,RDT=RDT)
    fa=factorize(a)
    AI=joLooseLinearFunctionInplaceAll(mn,mn,
        (y1,x1)->mul!(y1,a,x1),
        (y2,x2)->mul!(y2,transpose(a),x2),
        (y3,x3)->mul!(y3,adjoint(a),x3),
        (y4,x4)->mul!(y4,conj(a),x4),
        (y5,x5)->ldiv!(y5,fa,x5),
        (y6,x6)->ldiv!(y6,transpose(fa),x6),
        (y7,x7)->ldiv!(y7,adjoint(fa),x7),
        (y8,x8)->ldiv!(y8,factorize(conj(a)),x8),
        DDT,RDT)

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=A*x
    my=A*mx
    @! Y=AI*x
    @test Y≈y
    @! MY=AI*mx
    @test MY≈my

    y=conj(A)*x
    my=conj(A)*mx
    @! Y=conj(AI)*x
    @test Y≈y
    @! MY=conj(AI)*mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=transpose(A)*x
    my=transpose(A)*mx
    @! Y=transpose(AI)*x
    @test Y≈y
    @! MY=transpose(AI)*mx
    @test MY≈my

    y=adjoint(A)*x
    my=adjoint(A)*mx
    @! Y=adjoint(AI)*x
    @test Y≈y
    @! MY=adjoint(AI)*mx
    @test MY≈my

    y=conj(A)\x
    my=conj(A)\mx
    @! Y=conj(AI)\x
    @test Y≈y
    @! MY=conj(AI)\mx
    @test MY≈my

    x=rand(RDT,mn)
    mx=rand(RDT,mn,mn)
    Y=zeros(DDT,mn)
    MY=zeros(DDT,mn,mn)

    y=A\x
    my=A\mx
    @! Y=AI\x
    @test Y≈y
    @! MY=AI\mx
    @test MY≈my

    x=rand(DDT,mn)
    mx=rand(DDT,mn,mn)
    Y=zeros(RDT,mn)
    MY=zeros(RDT,mn,mn)

    y=transpose(A)\x
    my=transpose(A)\mx
    @! Y=transpose(AI)\x
    @test Y≈y
    @! MY=transpose(AI)\mx
    @test MY≈my

    y=adjoint(A)\x
    my=adjoint(A)\mx
    @! Y=adjoint(AI)\x
    @test Y≈y
    @! MY=adjoint(AI)\mx
    @test MY≈my
end

end


text:
# file for all point cloud drawing related functions

function reconstruct(dc::DepthCamera, depth::Array{Float64})
  s = dc.skip
  depth_sampled = depth[1:s:end,1:s:end]
  # assert(depth_sampled.shape == self.xs.shape)
  r,c = size(dc.xs)

  ret = Array{Float64,3}(r,c,3)
  ret[:,:,1] = dc.xs .* depth_sampled
  ret[:,:,2] = dc.ys .* depth_sampled
  ret[:,:,3] = depth_sampled
  return ret
end

function prepcolordepthcloud!( cvid::Int,
      X::Array;
      rgb::Array=Array{<:Colorant,2}(),
      skip::Int=4,
      maxrange::Float64=4.5 )
  #
  pointcloud = nothing
  pccols = nothing
  havecolor = size(rgb,1) > 0
  if typeof(X) == Array{Float64,3}
    r,c,h = size(X)
    Xd = X[1:skip:r,1:skip:c,:]
    rd,cd,hd = size(Xd)
    mask = Xd[:,:,:] .> maxrange
    Xd[mask] = Inf

    rgbss = havecolor ? rgb[1:skip:r,1:skip:c] : nothing
    # rgbss = rgb[1:4:r,1:4:c,:]./255.0
    pts = Vector{Vector{Float64}}()
    pccols = Vector()
    for i in 1:rd, j in 1:cd
      if !isnan(Xd[i,j,1]) && Xd[i,j,3] != Inf
        push!(pts, vec(Xd[i,j,:]) )
        havecolor ? push!(pccols, rgbss[i,j] ) : nothing
        # push!(pccols, RGB(rgbss[i,j,3], rgbss[i,j,2], rgbss[i,j,1]) )
      end
    end
    pointcloud = PointCloud(pts)
  elseif typeof(X) == Array{Array{Float64,1},1}
    pointcloud = PointCloud(X)
    pccols = rgb # TODO: refactor
  elseif size(X,1)==0
    return nothing
  else
    error("dont know how to deal with data type=$(typeof(X)),size=$(size(X))")
  end
  if havecolor
    pointcloud.channels[:rgb] = pccols
  else
    #submap colors
    smc = submapcolor(cvid, length(X))
    pointcloud.channels[:rgb] = smc
  end
  return pointcloud
end

function drawpointcloud!(vis, # ::DrakeVisualizer.Visualizer
                         poseswithdepth::Dict,
                         vsym::Symbol,
                         pointcloud,
                         va,
                         param::Dict,
                         sesssym::Symbol;
                         # imshape=(480,640),
                         wTb::CoordinateTransformations.AbstractAffineMap=
                               Translation(0,0,0.0) ∘ LinearMap(
                               CoordinateTransformations.Quat(1.0, 0, 0, 0))   )
                         # bTc::CoordinateTransformations.AbstractAffineMap=
                         #       Translation(0,0,0.6) ∘ LinearMap(
                         #       CoordinateTransformations.Quat(0.5, -0.5, 0.5, -0.5))  )
  #

  pcsym = Symbol(string("pc_", va != "none" ? va : "ID"))
  setgeometry!(vis[sesssym][pcsym][vsym][:pose], Triad())
  settransform!(vis[sesssym][pcsym][vsym][:pose], wTb) # also updated as parallel track
  setgeometry!(vis[sesssym][pcsym][vsym][:pose][:cam], Triad())
  settransform!(vis[sesssym][pcsym][vsym][:pose][:cam], param["bTc"] )
  setgeometry!(vis[sesssym][pcsym][vsym][:pose][:cam][:pc], pointcloud )

  # these poses need to be update if the point cloud is to be moved
  if !haskey(poseswithdepth,vsym)
    thetype = typeof(vis[sesssym][pcsym][vsym][:pose])
    poseswithdepth[vsym] = Vector{ thetype }()
  end
  push!(poseswithdepth[vsym], vis[sesssym][pcsym][vsym][:pose])

  nothing
end

"""
    $(SIGNATURES)
Draw point cloud on pose.
xTc -> pose to camera transform
""" #TODO: confirm xTc or cTx
function visPointCloudOnPose!(vis::Visualizer, x::Symbol, pointcloud::PointCloud, sessionId, xTc::SE3 = SE3([0,0,0],I))::Nothing
	# TODO: Cleanup and make parameter.
	material = PointsMaterial(size=0.02) # Used to set size of particles
	setobject!(vis[Symbol(sessionId)][:poses][x][:pc], pointcloud, material)
	trans = Translation(xTc.t[1], xTc.t[2], xTc.t[3])∘LinearMap(Quat(xTc.R.R))
	settransform!(vis[Symbol(sessionId)][:poses][x][:pc], trans)
	return nothing
end
function visPointCloudOnPose!(botvis::BotVis2, x::Symbol, pointcloud::PointCloud, sessionId, xTc::SE3 = SE3([0,0,0],I))::Nothing
   visPointCloudOnPose!(botvis.vis, x, pointcloud, sessionId, xTc)
end

function drawPointCloudonPose!(botvis::BotVis2, x::Symbol, pointcloud::PointCloud, xTc::SE3 = SE3([0,0,0],I))::Nothing
    @warn "drawPointCloudonPose! decprecated, use visPointCloudOnPose! instead."
    visPointCloudOnPose!(botvis, x, pointcloud, xTc)
end

"""
    $(SIGNATURES)

Main plugin callback function
"""
function pointCloudPlugins(vis::MeshCat.Visualizer,
                               params::Dict{Symbol, Any},
                               rose_fgl)
    # Assume static list for now...
	if !haskey(params, :dataEntries)
		@info "Getting data entries to cache..."
		entries = getDataEntriesForSession()
		params[:dataEntries] = entries
	end
	if !haskey(params, :depthsDrawn)
		params[:depthsDrawn] = Dict{String, Any}()
	end
	if !haskey(params, :robotConfig)
		params[:robotConfig] = getRobotConfig()
	end

	@show allVarLabels = collect(keys(params[:cachevars]))

	for varLabel in string.(allVarLabels)
		@info "Rendering potential data from $varLabel..."
		if haskey(params[:dataEntries], varLabel)
			dataSet = params[:dataEntries][varLabel]
			dataKeys = Dict(map(i -> i.id, dataSet) .=> dataSet)
			if haskey(dataKeys, "Depth") && haskey(dataKeys, "Sensor") #kQi
				# TODO: Check formats and dimensions and everything that is going to blow up here...
				# literally... everything.
				depthElem = getRawData(getNode(varLabel), dataKeys["Depth"])
				depth = collect(reinterpret(UInt16, base64decode(depthElem)))
				depth = reshape(depth, (640, 480))

				sensorElem = getRawData(getNode(varLabel), dataKeys["Sensor"])
				sensorElem = JSON.parse(sensorElem)
				kQi = map(a -> Float64(a), sensorElem["kQi"]) #Vector{Float64}
				kTc = (SE3([0,0,0], Quaternion(kQi)))
				trans = Translation([0,0,0])∘LinearMap(Quat(kTc.R.R))

				dcam = Arena.CameraModel(640, 480, 387.205, [322.042, 238.544])
                pointcloud = cloudFromDepthImageClampZ(depth, dcam, trans, maxrange=3f0, clampz = [-1f0,1f0]) # , colmap=repeatedColorMap
                Arena.visPointCloudOnPose!(vis, Symbol(varLabel), pointcloud, params[:sessionId])
			else
				@info " --- Skipping because: Depth exists? $(haskey(dataKeys, "Depth") ? "TRUE" : "FALSE"); Sensor exists? $(haskey(dataKeys, "Sensor") ? "TRUE" : "FALSE")"
			end
		else
			@info " --- No data entries for this label..."
		end
	end

	# ID of depth data ~ Depth
	# Format of depth data - hopefully UInt16
	# Dimensions - width x height
	# Camera transform - intrinsic = camModel
	# Camera transform - extrinsic = kQi
    return nothing
end


text:
function reciprocal_cycle_len(n::Int)
	seen=Dict{Int,Int}()
    x=1
    index=1
	while true

		if x in keys(seen)
			return index-seen[x]+1
		else
			seen[x]=index
			x=x*10 %n
		end
		index+=1
	end
end

function compute()
	ans=0
	for i =1:1000
		ans=max(ans,reciprocal_cycle_len(i))
	end

	return ans

end

println(compute())


text:
# Query adjacency matrix in a database table
using PyPlot

# Create a starting set of vertices
Nv0 = 100
v0 = ceil.(Int,10000*rand(Nv0,1))

# Convert to string list.
v0str = string.(v0)

# Get degrees of vertices.
Adeg = str2num(TadjDeg[v0str,:])

# Select vertices in an out degree range.
degMin = 5;  degMax = max(100,minimum(getval(Adeg[:,"OutDeg,"]))+1)
v1str = getrow((Adeg[:,"OutDeg,"] > degMin) < degMax )

# Get vertex neighbors.
A = logical(Tadj[v1str,:])

figure()
spy(A)
xlabel("end vertex")
ylabel("start vertex");


text:
"""
$(TYPEDEF)
"""
mutable struct Point3D
    # Coordinates
    x::Float64
    y::Float64
    z::Float64

    # Velocity at those coordinates
    ix::Int
    iy::Int
end

"""
$(SIGNATURES)
"""
function dist(p1::Point3D, p2::Point3D)
    dx = p2.x - p1.x
    dy = p2.y - p1.y
    return sqrt(dx * dx + dy * dy)
end

"""
$(SIGNATURES)

A midpoint is the average between two points
"""
midpoint(p1::Point3D, p2::Point3D) =
    Point3D(0.5p1.x + 0.5p2.x, 0.5p1.y + 0.5p2.y, 0.5p1.z + 0.5p2.z, 0, 0)

"""
$(SIGNATURES)

Barycentre of three points.
"""
centroid(p1::Point3D, p2::Point3D, p3::Point3D) = Point3D(
    (p1.x + p2.x + p3.x) / 3.0,
    (p1.y + p2.y + p3.y) / 3.0,
    (p1.z + p2.z + p3.z) / 3.0,
    0,
    0,
)

"""
$(TYPEDEF)
"""
struct Triangle
    pt1::Int64
    pt2::Int64
    pt3::Int64

    Triangle() = new(0, 0, 0)
    Triangle(pt1, pt2, pt3) = new(pt1, pt2, pt3)
end

"""
$(TYPEDEF)
"""
struct Mesh
    nodes::Vector{Point3D}
    nodeArray::Matrix{Point3D}
    triangles::Vector{Triangle}
    width::Int
    height::Int
end

"""
$(SIGNATURES)
"""
function triangle_area(mesh::Mesh, index::Int)
    triangle = mesh.triangles[index]
    pt1 = mesh.nodes[triangle.pt1]
    pt2 = mesh.nodes[triangle.pt2]
    pt3 = mesh.nodes[triangle.pt3]

    return triangle_area(pt1, pt2, pt3)
end

"""
$(SIGNATURES)
"""
function triangle_area(p1::Point3D, p2::Point3D, p3::Point3D)
    a = dist(p1, p2)
    b = dist(p2, p3)
    c = dist(p3, p1)
    s = (a + b + c) / 2.0

    return sqrt(s * (s - a) * (s - b) * (s - c))
end

max_stars_repo_name: Mehul-Bagaria/causticsEngineering

max_issues_repo_name: Mehul-Bagaria/causticsEngineering

max_forks_repo_name: Mehul-Bagaria/causticsEngineering


text:
module PyImage
    using PyCall, PyPlot

    const PIL = pyimport("PIL")
    const torchvision = pyimport("torchvision")
    const plt = PyPlot.matplotlib.pyplot
    const io = pyimport("io")

    export transforms
    const transforms = torchvision.transforms

    # 画像をファイルから読み込み
    open(filepath::String)::PyObject = PIL.Image.open(filepath)

    # bytesから画像生成
    ## byte: UInt8型1次元配列
    open(bytes::Array{UInt8, 1})::PyObject = PIL.Image.open(io.BytesIO(bytes))

    # PIL.Imageオブジェクトに対して変換処理を行う関数を生成
    ## @param compose: PIL.transformsの配列（変換方法の指定）
    ## @return image::PyObject（PIL.Image） -> data::Array{AbstractFloat, 3}（色数×高さ×幅）
    make_transformer(compose::Array{PyObject, 1})::(PyObject -> Array{AbstractFloat, 3}) = begin
        transformer = transforms.Compose(hcat(compose, [transforms.ToTensor()]))
        (image::PyObject) -> transformer(image).numpy()
    end

    # 画像表示
    show(image::Any)::PyObject = plt.imshow(image)
end

max_stars_repo_name: amenoyoya/pytorch-learning

max_issues_repo_name: amenoyoya/pytorch-learning

max_forks_repo_name: amenoyoya/pytorch-learning


text:
@testset "filter.jl" begin
    println("             > filter.jl")
    m_out  = 1024
    m_in   = 1024

    state = BinomialState(128, m_out, m_in)
    model = BinomialGridModel(
        m_out,
        1:5,
        LinRange(0.05,0.95,5),
        LinRange(0.1,2,5),
        LinRange(0.05,2,5),
        LinRange(0.05,2,5)
    )

    fstate = NestedParticleState(state, model)
    filter = NestedParticleFilter(12)
    obs    = BinomialObservation(0.3f0, 0.1f0)

    if true
        println("")
        println("Benchmarking one filter update step: should take less than 10ms")
        display(@benchmark CUDA.@sync update!($fstate, $obs, $filter))
        println("")
        println("")
    end
end


text:
function weak_phs1(N,a,b)
	N1 = N+1   # variables x1, e1 and v1
	N2 = N   	#    variables x2, e2 and v2

	x1,w1,P = lglnodes(N1-1,a,b)    # discretization of x1 variables
	if N2 > 1
		x2,w2 = lgwt(N2,a,b)
	else
		x2,w2 = (a+b)/2, (b-a);
	end
	#M = massmatrix(ei,xi,we);
	#xquad, wquad = lgwt(Ne,a,b)\
	#M = massmatrix(ei,xi, xquad, wquad);
	M1 = massmatrix(x1,x1, w1);
	M2 = massmatrix([x2],[x2], [w2]);

	D = dermatrix(x1,x2,1)'*M2 ;
	p0 = map(i->leg_pol(a,x1,i), 1:length(x1))
	pL = map(i->leg_pol(b,x1,i), 1:length(x1))

    phJ = [zeros(N1,N1) D;
	   -D' zeros(N2, N2)];
    phB = [p0 -pL; zeros(N2,2)];
    phD = zeros(2,2);

	Q = blkdiag(inv(M1),inv(M2));
	#pflow = Collocation_points(x1, w1)
	#peffort = Collocation_points(x2, w2)

	p = Phs(phJ, phB, phD, Q)
	#p.Q = Q
	#p.disc_data = Discrete_domain(pflow, peffort)
	#func(x :: Array) = (0.5*x'*p.Q*x)[1]
	#p.Hamiltonian = func
	#p.GradHam = ForwardDiff.gradient(p.Hamiltonian)
	p
end

function weak_phs2(N,a,b)
	N1 = N   # variables x1, e1 and v1
	N2 = N+1   	#    variables x2, e2 and v2

	x2,w2,P = lglnodes(N2-1,a,b)    # discretization of x1 variables
	if N1 > 1
		x1,w1 = lgwt(N1,a,b)
	else
		x1,w1 = (a+b)/2, (b-a);
	end
	#M = massmatrix(ei,xi,we);
	#xquad, wquad = lgwt(Ne,a,b)\
	#M = massmatrix(ei,xi, xquad, wquad);
	M1 = massmatrix([x1],[x1], [w1]);
	M2 = massmatrix(x2,x2, w2);

	D = dermatrix(x2,x1,1)'*M1 ;
	p0 = map(i->leg_pol(a,x2,i), 1:length(x2))
	pL = map(i->leg_pol(b,x2,i), 1:length(x2))

    phJ = [zeros(N1,N1) -D';
	   D zeros(N2, N2)];
    phB = [zeros(N1,2); p0 -pL];
    phD = zeros(2,2);

		Q = blkdiag(inv(M1),inv(M2));

	#pflow = Collocation_points(x1, w1)
	#peffort = Collocation_points(x2, w2)

	p = Phs(phJ, phB, phD, Q)
	#p.Q = Q
	#p.disc_data = Discrete_domain(pflow, peffort)
	#func(x :: Array) = (0.5*x'*p.Q*x)[1]
	#p.Hamiltonian = func
	#p.GradHam = ForwardDiff.gradient(p.Hamiltonian)
	p
end

#coupled_gyrator(ph1 :: Phs, ports1, ph2 :: Phs, ports2, couple_matrix)
function finelem(Nelem,Nint,a,b)
	dx = (b-a)/Nelem;

	p = weak_phs1(Nint,0,dx)
	for i = 2:Nelem
		if (i%2 == 1)
			pn = weak_phs1(Nint,0,dx)
		else
			pn = weak_phs2(Nint,0,dx)
		end
		p = PortHamiltonian.coupled_gyrator(p,[2],pn,[1],[-1])
	end
	p
end

max_stars_repo_name: flavioluiz/PortHamiltonian

max_issues_repo_name: flavioluiz/PortHamiltonian

max_forks_repo_name: flavioluiz/PortHamiltonian


text:
module ScoreDrivenModels

using Distributions, Optim, SpecialFunctions
using LinearAlgebra, Printf

import Base: length, deepcopy, show

const SCALINGS  = [0.0, 1/2, 1.0]
const BIG_NUM = 1e8
const SMALL_NUM = 1e-8

# Core files
include("abstracts.jl")
include("utils.jl")
include("link_functions.jl")
include("score.jl")
include("model.jl")
include("initial_params.jl")
include("simulate.jl")
include("diagnostics.jl")
include("univariate_score_driven_recursion.jl")
include("MLE.jl")
include("prints.jl")

# Optimization methods
include("opt_methods/common_methods.jl")
include("opt_methods/LBFGS.jl")
include("opt_methods/IPNewton.jl")
include("opt_methods/NelderMead.jl")

# Distributions
include("distributions/non_native_dists.jl")
include("distributions/common_interface.jl")
include("distributions/beta.jl")
include("distributions/betafour.jl")
include("distributions/exponential.jl")
include("distributions/gamma.jl")
include("distributions/logitnormal.jl")
include("distributions/lognormal.jl")
include("distributions/negativebinomial.jl")
include("distributions/normal.jl")
include("distributions/poisson.jl")
include("distributions/tdist.jl")
include("distributions/tdistlocationscale.jl")
include("distributions/weibull.jl")

end


text:
using POMDPModels
using POMDPs
using POMDPModelTools
using StaticArrays
using Random
using DiscreteValueIteration
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationRandomStrategy
using POMDPLinter: show_requirements, get_requirements
using Test

@testset "all" begin

    println("Testing Requirements")
    gifa = LocalGIFunctionApproximator(RectangleGrid([0.0, 1.0], [0.0, 1.0]))
    @test_skip @requirements_info LocalApproximationRandomSolver(gifa) SimpleGridWorld()
    show_requirements(get_requirements(POMDPs.solve, (LocalApproximationRandomSolver(gifa), SimpleGridWorld())))
end


text:
using Reexport

@reexport module JointDistanceDistribution
    import HypothesisTests: OneSampleTTest, pvalue
    export OneSampleTTest, pvalue
    include("joint_distance_distribution.jl")
end