text:
# This file contains the necessary ingredients to create a PackageManager for BinDeps
using BinDeps

type EnvManager{T} <: BinDeps.PackageManager
    packages::Vector{Compat.UTF8String}
end

"Manager for root environment"
const Manager = EnvManager{Symbol(PREFIX)}

function Base.show{T}(io::IO, manager::EnvManager{T})
    write(io, "Conda packages: ", join(manager.packages, ", "))
end

BinDeps.can_use(::Type{EnvManager}) = true

function BinDeps.package_available{T}(manager::EnvManager{T})
    pkgs = manager.packages
    # For each package, see if we can get info about it. If not, fail out
    for pkg in pkgs
        if !exists(pkg, T)
            return false
        end
    end
    return true
end

BinDeps.libdir{T}(m::EnvManager{T}, ::Any) = lib_dir(T)
BinDeps.bindir{T}(m::EnvManager{T}, ::Any) = bin_dir(T)

BinDeps.provider{T, S<:String}(::Type{EnvManager{T}}, packages::Vector{S}; opts...) = EnvManager{T}(packages)
BinDeps.provider{T}(::Type{EnvManager{T}}, packages::String; opts...) = EnvManager{T}([packages])

function BinDeps.generate_steps(dep::BinDeps.LibraryDependency, manager::EnvManager, opts)
    pkgs = manager.packages
    if isa(pkgs, AbstractString)
        pkgs = [pkgs]
    end
    ()->install(pkgs, manager)
end

function install{T}(pkgs, manager::EnvManager{T})
    for pkg in pkgs
        add(pkg, T)
    end
end


text:
import HDF5, MPI

import RIOPA: IO

struct HDF5IOTag <: IO.IOTag end

function perform_step(::HDF5IOTag, ds::DataSet)
    worldrank = MPI.Comm_rank(MPI.COMM_WORLD)
    stepname = ds.cfg.basename * "_" * lpad(ds.curr_step, 5, '0')
    for i = 1:length(ds.streams)
        streamCfg = ds.cfg.streams[i]
        pathname = stepname * "/" * streamCfg.name
        mkpath(pathname)
        filename = "D_" * lpad(worldrank, 5, '0') * ".h5"
        HDF5.h5write(pathname * "/" * filename, "data", ds.streams[i].vec)
    end
end

IO.add("HDF5", HDF5IOTag())

ext: jl


text:
abstract type TradingStructures end
### TODO (later) use a suitable Blotter/FinancialSeries abstract type

"Order type"
struct Order <: TradingStructures
  id::String # unique order id for submission/tracking
  quantity::Int64 # abs. value of position change targeted
  price::Float64 # limit price, set to NaN for market order
  side::Symbol # :buy :sell (no need for :short <=> {position = 0 with sell})
  ordertype::Symbol # :market :limit
  status::Symbol # :pending :complete :cancelled
  trigger::String # custom string for tracking what triggered the order

  function Order(id::String, quantity::Int64, price::Float64,
                 side::Symbol, ordertype::Symbol, status::Symbol,
                 trigger::String)
    (side == :buy || side == :sell) || error("Unknown order side")
    (ordertype == :market || ordertype == :limit) || error("Unsupported order type")
    if !(status == :pending || status == :complete || status == :cancelled)
      error("Unsupported order status")
    end
    new(id, quantity, price, side, ordertype, status, trigger)
  end
end

"Empty order: no quantity"
emptyorder() = Order("empty", 0, NaN, :buy, :market, :complete, "init")

"Check if order status is `:pending`"
ispending(orde::Order) = orde.status == :pending

"Change order status to `:cancelled`"
function setcancelled!(orde::Order)
  orde.status = :cancelled
  return
end

"Signed position change in the Order object"
function getorderposchg(orde::Order)
  if orde.side == :buy
    return orde.quantity
  end
  if orde.side == :sell
    return -orde.quantity
  end
  error("Unknown order side")
end

const Blotter = @compat Dict{DateTime,Tuple{Int64,Float64}}

const Targ = @compat Tuple{Int64,Vector{Float64}}

const OHLC = @compat Tuple{DateTime,Vector{Float64}}

"""
Initialize empty blotter as an associative collection
`DateTime => (Qty::Int64, FillPrice::Float64)`
"""
emptyblotter() = Blotter()

# "Completed (closed) trade type."
# immutable ClosedTrade
#   topen::DateTime
#   tclose::DateTime
#   popen::Float64
#   pclose::Float64
#   tradeside::Symbol
#   tradequantity::Int64
# end

# """
# Array-like type keeping all open trades,
# thus holding current (actual) position information with open prices.
# Multiple entries may appear when an order is filled in parts.
# """
# immutable OpenTrades <: TradingStructures
#   quantities::Vector{Int64} # abs. value of outstanding position
#   openprices::Vector{Float64}
#   sides::Vector{Symbol} # :long :short

#   ### TODO inner constructor checks (similar to TimeArray in TimeSeries.jl)
# end

# "OpenTrades object without outstanding position"
# emptyopentrades() = OpenTrades([0], [NaN], [:long])

# "Get total outstanding position from `OpenTrades`-object."
# getopenposition(otrades::OpenTrades) = sum(otrades.quantities)

### use array of Trade-objects instead of what is below?
### (append closed trades to it)

# """
# Array type keeping trades information, both open and completed.
# Since it contains open trades, it holds
# current (actual) position information with open prices.
# Multiple open entries may appear when an order is filled in parts.
# """
# immutable TradesArray <: TradingStructures
#   entertimes::Vector{DateTime}
#   exittimes::Union(Vector{DateTime}, Symbol) # :open of exit DateTime
#   quantities::Vector{Int64} # abs. value of outstanding position
#   openprices::Vector{Float64}
#   closeprices::Vector{Float64}
#   sides::Vector{Symbol} # :long :short

#   ### TODO inner constructor checks (similar to TimeArray in TimeSeries.jl)
# end

# "Empty TradesArray object"
# emptytrades() = TradesArray(Array(DateTime, 0), Array(DateTime, 0),
#                             Array(Int64, 0),
#                             Array(Float64, 0), Array(Float64, 0),
#                             Array(Symbol, 0))

# "Get total outstanding position from `TradesArray`-object."
# getopenposition(trarr::TradesArray) = sum(otrades.quantities)


text:
export @any, @all

macro any(comp)
    if !isa(comp, Expr) || (comp.head != :comprehension)
        error("@any may only precede a comprehension")
    end
    compgen = comp.args[1]
    loop = Expr(:for, esc(compgen.args[2]),
                          quote
                              if $(esc(compgen.args[1]))
                                  result = true
                                  break
                              end
                          end)
    quote
        result = false
        $loop
        result
    end
end

macro all(comp)
    if !isa(comp, Expr) || (comp.head != :comprehension)
        error("@all may only precede a comprehension")
    end
    compgen = comp.args[1]
    loop = Expr(:for, esc(compgen.args[2]),
                          quote
                              if !$(esc(compgen.args[1]))
                                  result = false
                                  break
                              end
                          end)
    quote
        result = true
        $loop
        result
    end
end

@generated function blend{T1,T2}(b::StaticArray{Bool}, a1::StaticArray{T1}, a2::StaticArray{T2})
    if !(size(b) == size(a1) == size(a2))
        error("Dimensions must match. Got sizes $(size(b)), $(size(a1)), and $(size(a2))")
    end
    newtype = :(similar_type($a1, promote_type(T1, T2)))
    exprs = [:(b[$j] ? a1[$j] : a2[$j]) for j = 1:length(b)]
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

(::Type{SVector})(x::Vector) = SVector(x...)
(::Type{MVector})(x::Vector) = MVector(x...)

# Inspired by SparseVectors.jl
import Base: length, size, nnz, countnz
import Base.SparseArrays: nonzeros, nonzeroinds
export viewcol, nonzeroinds

typealias CVecView{T} SubArray{T,1,Vector{T},Tuple{UnitRange{Int64}},true}

immutable SparseVectorView{Tv,Ti<:Integer} <: AbstractSparseVector{Tv,Ti}
    n::Int                  # the number of elements
    nzind::CVecView{Ti}     # the indices of nonzeros
    nzval::CVecView{Tv}     # the values of nonzeros

    function SparseVectorView(n::Integer, nzind::CVecView{Ti}, nzval::CVecView{Tv})
        n >= 0 || throw(ArgumentError("The number of elements must be non-negative."))
        size(nzind) == size(nzval) || throw(DimensionMismatch("The lengths of nzind and nzval are inconsistent."))
        new(convert(Int, n), nzind, nzval)
    end
end
SparseVectorView{Tv,Ti}(n::Integer, nzind::CVecView{Ti}, nzval::CVecView{Tv}) = SparseVectorView{Tv,Ti}(n, nzind, nzval)

function viewcol(x::SparseMatrixCSC, j::Integer)
    1 <= j <= x.n || throw(BoundsError())
    r1 = convert(Int, x.colptr[j])
    r2 = convert(Int, x.colptr[j+1]) - 1
    rgn = r1:r2
    SparseVectorView(x.m, view(x.rowval, rgn), view(x.nzval, rgn))
end
length(x::SparseVectorView) = x.n
size(x::SparseVectorView) = (x.n,)
nnz(x::SparseVectorView) = length(x.nzval)
countnz(x::SparseVectorView) = countnz(x.nzval)
nonzeros(x::SparseVectorView) = x.nzval
nonzeroinds(x::SparseVectorView) = x.nzind

# angles
mod2piF{T}(x::T) = mod(x, 2*T(pi))
function adiff(x,y)
    x = mod2piF(x)
    y = mod2piF(y)
    abs(x - y) <= pi ? x - y : (x > y ? x - y - 2pi : x - y + 2pi)
end

include("vec2Dutils.jl")


text:
# Try to select a reasonable set of sigma values that should be used
# for the adaptive shrinkage model based on the values of x (the noisy
# observations) and s (the standard error in the observations). The
# return value is a vector of sigma values.
function autoselectmixsd(x::Array{Float64,1},
                         s::Array{Float64,1} = ones(size(x));
                         gridmult::Float64 = 1.4,
                         nv::Int = 0)

  # Get the number of samples.
  n = length(x);

  # Check input "s"---it should be the same length.
  if length(s) != n
    throw(ArgumentError("Arguments \"x\" and \"s\" should have the same" *
                        "length"))
  elseif any(s .<= 0)
    throw(ArgumentError("All elements of \"s\" should be positive"))
  end

  # Check input "gridmult".
  if gridmult < 0
    throw(ArgumentError("Input \"gridmult\" should be non-negative"))
  end

  # Check input "nv".
  if !(nv == 0 || nv > 1)
    throw(ArgumentError("Input \"nv\" should be 0, or greater than 1"))
  end

  # Determine the minimum and maximum sigma settings.
  smin = sqrt(minimum(s.^2))/10;
  if all(x.^2 .<= s.^2)

    # This deals with the occasional boundary case.
    smax = 10 * smin;
  else

    # This is, roughly, the largest value you'd want to use.
    smax = 2 * sqrt(maximum(x.^2 - s.^2));
  end

  # Choose the grid of sigmas.
  if nv > 0
    return vcat(0,logspace(log10(smin),log10(smax),nv - 1));
  elseif gridmult == 0
    return vcat(0,smax)
  else
    m = ceil(Int,log2(smax/smin)/log2(gridmult));
    return vcat(0,logspace(log10(smin),log10(smax),m + 1));
  end
end

# Compute the n x k conditional likelihood matrix, where n is the
# number of samples and k is the number of mixture components, for the
# case when the likelihood is univariate normal and prior is a mixture
# of univariate normals.
#
# Entry (i,j) of the conditional likelihood matrix is equal to
# N(0,s[i]^2 + sd[j]^2), the normal density with zero mean and
# variance s[i]^2 + sd[j]^2.
#
# If normalize = true, each row of the likelihood matrix is divided by
# the largest value in the row. After normalization, the largest value
# in each row is 1.
function normlikmatrix(x::Array{Float64,1},
                       s::Array{Float64,1} = ones(size(x));
                       sd::Array{Float64,1} = autoselectmixsd(x,s),
                       normalizerows::Bool = true)

  # Get the number of samples (n) and the number of prior mixture
  # components (k).
  n = length(x);
  k = length(sd);

  # Check input "s"---it should be the same length at x.
  if length(s) != n
    throw(ArgumentError("Arguments \"x\" and \"s\" should have the same" *
                        "length"))
  elseif any(s .<= 0)
    throw(ArgumentError("All elements of \"s\" should be positive"))
  end

  # Check input "sd".
  if any(sd .< 0)
    throw(ArgumentError("All elements of \"sd\" should be non-negative"))
  end

  # Compute the n x k matrix of standard deviations.
  S = sqrt.((s.^2) .+ (sd.^2)');

  # Compute the log-densities, and normalize the rows, if requested.
  L = -(x./S).^2/2 - log.(S) .- log(2*pi)/2;
  if normalizerows

    # This is the same as
    #
    #   L = L - repmat(maximum(L,2),1,k);
    #
    # but uses memory more efficiently to complete the operation.
    L = broadcast(-,L,maximum(L,dims = 2));
  end
  return exp.(L)
end

# Implements the "logspace" function that was available in an
# earlier version of Julia.
function logspace(a, b, n)
  return exp10.(range(a,stop = b,length = n))
end

max_stars_repo_name: stephenslab/mixsqp-paper

max_issues_repo_name: stephenslab/mixsqp-paper

max_forks_repo_name: stephenslab/mixsqp-paper


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using StaticArrays: SVector, SMatrix

"""
    plane_from_points(points::SMatrix{D, N, P}}) ->  centroid, normal, local_to_world transform

Points to be fitted are assumed to be stored by column in the `points` matrix.
Estimate the best fitting plane for a set of points in 3D.
`D` is the dimension of the plane.
`N` is the number of points to fit.
`P` is the number type used to represent points.
"""
function plane_from_points(points::SMatrix{D, N, P}) where {D, N, P <: Real}
    center = Statistics.mean(points, dims = 2) #compute average of columns of point matrix.

    u, _, _ = svd(points .- center)

    #always want a rotation matrix.
    if det(u) < 0
        u = SMatrix{D, D}(u[:, 1:(end - 1)]..., -u[:, end]...) #change sign of last singular vector to turn into rotation matrix
    end

    normal = u[:, 3]             # The two largest singular vectors lie in the plane that is the best fit to the points, i.e., that accounts for the largest fraction of variance in the set of points. The smallest singular vector is perpendicular to this plane.

    # make sure the normal is pointing consistently to positive Z direction of local coordinate frame
    if dot(normal, unitZ3()) < P(0.0)
        normal = normal * P(-1.0)
    end

    return SVector(center), SVector(normal), u     # convert from SMatrix to SVector, return u matrix to use as local coordinate frame for the plane.
end

function plane_from_points(points::AbstractMatrix{P}) where {P <: Real}
    dims = size(points)
    temp = MMatrix{dims[1], dims[2], P}(points)
    return plane_from_points(SMatrix(temp))
end

"""only returns real roots"""
function quadraticroots(a::T, b::T, c::T) where {T <: Real}
    temp = b^2 - 4 * a * c

    if temp < zero(T)
        return nothing # no real roots so no ray cylinder intersection
    end

    radical = sqrt(temp)
    if b >= zero(T)
        radicalterm = -b - radical
        x1 = radicalterm / (2a)
        x2 = (2c) / radicalterm
    else
        radicalterm = -b + radical
        x1 = (2c) / radicalterm
        x2 = radicalterm / (2a)
    end

    return SVector{2, T}(x1, x2)
end

replprint(a) = show(IOContext(stdout), "text/plain", a)

@inline samepoint(pt1, pt2) = isapprox(pt1, pt2, atol = 1e-10)

function password()
    return String(
        collect(rand(('1':'9'..., 'A':'Z'..., 'a':'z'...)) for i in 1:14),
    )
end

# 1 if val>= 0 and -1 otherwise. Unlike Math.Sign which returns 0 if val==0

@inline sign(val::T) where {T <: Real} = val >= zero(T) ? 1 : -1

@inline function NaNsafeatan(x::T, y::T)::T where {T <: Real}
    if y == zero(T)
        if x == zero(T)
            return zero(T)
        elseif x > zero(T)
            return T(π / 2)
        else
            return -T(π / 2)
        end
    else
        return atan(x, y)
    end
end

@inline function NaNsafeasin(x::T)::T where {T <: Real}
    if x == zero(T)
        return zero(T)
    elseif x == one(T)
        return T(π / 2)
    elseif x == -one(T)
        return -T(π / 2)
    else
        return asin(x)
    end
end

@inline function NaNsafeacos(x::T)::T where {T <: Real}
    if x == zero(T)
        return T(π / 2)
    elseif x == one(T)
        return zero(T)
    elseif x == -one(T)
        return T(π)
    else
        return acos(x)
    end
end

# some place holders for package level function names.
# these names need to exist before any internal module can override them.
# NOTE is this really required?? Just rexport/bring into scope...
# function origin end


text:
#=
CF: according to git grep: this is not called.
Quite likely, it has been replaced by the other reduction functions
################################################################################
#
#  Size reduction
#
################################################################################

function _reduce_size(x::Array{T, 1}, prec::Int = 64) where T
  K = _base_ring(x[1])

  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  conlog = conjugates_arb_log(x[1], prec)

  A = zero_matrix(parent(conlog[1]), length(x), rr)

  B = zero_matrix(FlintZZ, nrows(A), ncols(A))

  for i in 1:rr
    A[1, i] = conlog[i]
  end

  Ar = base_ring(A)

  for i in 1:nrows(A)
    for j in 1:ncols(A)
      b, y = unique_integer(ceil(ldexp(A[i, j], 64)))
      @assert b
      B[i, j] = y
    end
  end

  L, U = lll_with_transform(B)
end

=#

################################################################################
#
#  Reduce units using LLL
#
################################################################################

function scaled_log_matrix(u::Array{T, 1}, pr::Int = 32) where T

  r,s = signature(_base_ring(u[1]))
  A = zero_matrix(FlintZZ, length(u), r + s)
  @vprint :UnitGroup 1 "starting prec in scaled_log_matrix: $pr\n"

  for i in 1:length(u)
    c = conjugates_arb_log(u[i], pr)
    for k in 1:length(c)
      #@show T
      @assert radiuslttwopower(c[k], -pr)
    end

    if any(x->radius(x) > 1e-9, c)  # too small...
      @vprint :UnitGroup 1 "increasing prec in scaled_log_matrix, now: $pr\n"
      pr *= 2
      if pr > 2^30
        error("cannot do lll on units")
        break
      end
      return scaled_log_matrix(u, pr)
    end
    for j in 1:length(c)
      tt = fmpz()
      t = ccall((:arb_mid_ptr, :libarb), Ptr{arf_struct}, (Ref{arb}, ), c[j])
      l = ccall((:arf_get_fmpz_fixed_si, :libarb), Int, (Ref{fmpz}, Ptr{arf_struct}, Int), tt, t, -pr)
      A[i, j] = tt
    end
  end
  return A, pr
end

function row_norm(A::fmpz_mat, i::Int)
  return sum([A[i,j]^2 for j=1:ncols(A)])
end

function row_norms(A::fmpz_mat)
  return [row_norm(A, i) for i=1:nrows(A)]
end

function reduce(u::Array{T, 1}, prec::Int = 32) where T
  @vprint :UnitGroup 1 "prec in reduce, now: $prec\n"
  r = length(u)
  if r == 0
    return u
  end

  while true
    A, prec = scaled_log_matrix(u, prec)

    L, U = lll_with_transform(A)
    @vprint :UnitGroup 2 "reducing units by $U\n"
    pA = prod(row_norms(A))
    pL = prod(row_norms(L))
    @vprint :UnitGroup 1 "reducing norms of logs from 2^$(nbits(pA)) -> 2^$(nbits(pL)), rat is $(Float64(1.0*pA//pL))\n"
    u = transform(u, transpose(U))
    if nbits(pL) >= nbits(pA)
    #  u = [compact_presentation(x, decom = factor(1*maximal_order(base_ring(x)))) for x = u]
      return u
    end
    @vprint :UnitGroup 1 "trying to reduce further...\n"
  end
end


text:
@testset "Poisson loss" begin
    N = 6
    img = abs.(randn((N, N)))

    @test 0.22741127776021886 ≈ poisson_aux([1.,2.], [3.,4.])

    @test all(0 .≈ Zygote.gradient(poisson_aux, img, img)[1])

    @test all( -1 .≈ Zygote.gradient(poisson_aux, [1.], [2.])[1])
end

@testset "Gaussian Loss" begin
    N = 6
    img = abs.(randn((N, N)))

    gauss = Gauss()
    @test 0 ≈ gauss(img, img)

    @test 0 ≈ gauss_aux(img, img)

    @test all(0 .≈ Zygote.gradient(gauss_aux, img, img)[1])

end

@testset "Scaled Gauss" begin
    N = 6
    img = abs.(randn((N, N)))

    scaled_gauss = ScaledGauss(0)

    @test 3.4094379124341003 ≈ scaled_gauss([5.], [2.])
    @test 3.4094379124341003 ≈ scaled_gauss_aux([5.], [2.], read_var=0)

    @test all( -0.75 .≈ Zygote.gradient((a, b) -> scaled_gauss_aux(a, b, read_var=1), [1.], [2.])[1])
end

@testset "Anscombe Loss" begin
    N = 6
    img = abs.(randn((N, N)))

    anscombe = Anscombe(1.0)
    @test 0 ≈ anscombe(img, img)

    @test 0 ≈ anscombe_aux(img, img,b=1)

    @test all(0 .≈ Zygote.gradient(im1 -> anscombe_aux(im1,img, b=1), img)[1])

end


text:
# This loads the PBMC8 dataset that is used for later tests

checkDir()

#create genData and data folder and change dir to dataPath
cwd = pwd()
if occursin("jenkins", homedir()) || "TRAVIS" in keys(ENV)
    genDataPath = mktempdir()
    dataPath = mktempdir()
else
    if !occursin("test", cwd)
        cd("test")
        cwd = pwd()
    end
    dataFolders = ["genData", "data"]
    for dir in dataFolders
        if !isdir(dir)
            mkdir(dir)
        end
    end
    genDataPath = cwd * "/genData"
    dataPath = cwd * "/data"
end

refDataPath = cwd * "/refData"
cd(dataPath)

# fetch the required data for testing and download the zip archive and unzip it
dataFiles = ["PBMC8_metadata.xlsx", "PBMC8_panel.xlsx", "PBMC8_fcs_files.zip"]
for f in dataFiles
    if !isfile(f)
        download("http://imlspenticton.uzh.ch/robinson_lab/cytofWorkflow/" * f, f)
        if occursin(".zip", f)
            run(`unzip PBMC8_fcs_files.zip`)
        end
    else
    end
end

# verify the data consistency using the stored checksums
fileNames = readdir()
csDict = Dict{String,Any}()
for f in fileNames
    if f[end-3:end] == ".fcs" || f[end-4:end] == ".xlsx"
        cs = bytes2hex(sha256(f))
        csDict[f] = cs
    end
end
csTest = JSON.parsefile(cwd * "/checkSums/csTest.json")
if csDict != csTest
    @error "Downloaded dataset does not match expectations, perhaps it is corrupted?" csDict csTest
    error("dataset checksum error")
end

md = DataFrame(XLSX.readtable("PBMC8_metadata.xlsx", "Sheet1", infer_eltypes = true)...)
panel = DataFrame(XLSX.readtable("PBMC8_panel.xlsx", "Sheet1", infer_eltypes = true)...)

antigens = panel[panel[:, :Lineage].==1, :Antigen]
_, fcsParams = loadFCSHeader(md[1, :file_name])
_, fcsAntigens = getMarkerNames(getMetaData(fcsParams))
cleanNames!(antigens)
cleanNames!(fcsAntigens)

di = loadFCSSet(:fcsData, md[:, :file_name], [myid()])

#prepare the data a bit
dselect(di, fcsAntigens, antigens)
cols = Vector(1:length(antigens))
dtransform_asinh(di, cols, 5)
dscale(di, cols)

pbmc8_data = distributed_collect(di)
undistribute(di)

@testset "load-time columns normalization" begin
    di = loadFCSSet(
        :fcsData,
        md[:, :file_name],
        [myid()],
        postLoad = selectFCSColumns(antigens),
    )
    dtransform_asinh(di, cols, 5)
    dscale(di, cols)
    @test pbmc8_data == distributed_collect(di)
    undistribute(di)
end


text:
# Container operations

# TODO: Macroize these definitions

function Base.push!(dv::DataVector, v::NAtype)
    resize!(dv.data, length(dv.data) + 1)
    push!(dv.na, true)
    return v
end

function Base.push!{S, T}(dv::DataVector{S}, v::T)
    push!(dv.data, v)
    push!(dv.na, false)
    return v
end

function Base.pop!(dv::DataVector)
    d, m = pop!(dv.data), pop!(dv.na)
    if m
        return NA
    else
        return d
    end
end

function Base.unshift!{T}(dv::DataVector{T}, v::NAtype)
    ccall(:jl_array_grow_beg, Void, (Any, UInt), dv.data, 1)
    unshift!(dv.na, true)
    return v
end

function Base.unshift!{S, T}(dv::DataVector{S}, v::T)
    unshift!(dv.data, v)
    unshift!(dv.na, false)
    return v
end

function Base.shift!{T}(dv::DataVector{T})
    d, m = shift!(dv.data), shift!(dv.na)
    if m
        return NA
    else
        return d
    end
end

function Base.splice!(dv::DataVector, inds::(@compat Union{Integer, UnitRange{Int}}))
    v = dv[inds]
    deleteat!(dv.data, inds)
    deleteat!(dv.na, inds)
    v
end

function Base.splice!(dv::DataVector, inds::(@compat Union{Integer, UnitRange{Int}}), ins::AbstractVector)
    # We cannot merely use the implementation in Base because this
    # needs to handle NA in the replacement vector
    v = dv[inds]
    m = length(ins)
    a = dv.data
    if m == 0
        deleteat!(a, inds)
        deleteat!(dv.na, inds)
        return v
    end

    n = length(a)
    f = first(inds)
    l = last(inds)
    d = length(inds)

    if VERSION >= v"0.5.0-dev+5022"
        if m < d
            delta = d - m
            i = (f - 1 < n - l) ? f : (l - delta + 1)
            Base._deleteat!(a, i, delta)
        elseif m > d
            delta = m - d
            i = (f - 1 < n - l) ? f : (l + 1)
            Base._growat!(a, i, delta)
        end
    else
        if m < d
            delta = d - m
            if f-1 < n-l
                Base._deleteat_beg!(a, f, delta)
            else
                Base._deleteat_end!(a, l-delta+1, delta)
            end
        elseif m > d
            delta = m - d
            if f-1 < n-l
                Base._growat_beg!(a, f, delta)
            else
                Base._growat_end!(a, l+1, delta)
            end
        end
    end

    for k = 1:m
        if !isna(ins, k)
            if isa(ins, DataVector)
                a[f+k-1] = ins.data[k]
            elseif isa(ins, PooledDataVector)
                a[f+k-1] = ins.pool[ins.refs[k]]
            else
                a[f+k-1] = ins[k]
            end
        end
    end

    splice!(dv.na, inds, isna(ins))
    v
end

function Base.deleteat!(dv::DataVector, inds)
    deleteat!(dv.data, inds)
    deleteat!(dv.na, inds)
    dv
end

# TODO: should this be an AbstractDataVector, so it works with PDV's?
function Base.map(f::Function, dv::DataVector)
    n = length(dv)
    res = DataArray(Any, n)
    for i in 1:n
        res[i] = f(dv[i])
    end
    return res
end

function Base.push!{T,R}(pdv::PooledDataVector{T,R}, v::NAtype)
    push!(pdv.refs, zero(R))
    return v
end

function Base.push!{S,R,T}(pdv::PooledDataVector{S,R}, v::T)
    v = convert(S,v)
    push!(pdv.refs, getpoolidx(pdv, v))
    return v
end

Base.pop!(pdv::PooledDataVector) = pdv.pool[pop!(pdv.refs)]

function Base.unshift!{T,R}(pdv::PooledDataVector{T,R}, v::NAtype)
    unshift!(pdv.refs, zero(R))
    return v
end

function Base.unshift!{S,R,T}(pdv::PooledDataVector{S,R}, v::T)
    v = convert(S,v)
    unshift!(pdv.refs, getpoolidx(pdv, v))
    return v
end

Base.shift!(pdv::PooledDataVector) = pdv.pool[shift!(pdv.refs)]

Base.reverse(x::AbstractDataVector) = x[end:-1:1]

function Base.splice!(pdv::PooledDataVector, inds::(@compat Union{Integer, UnitRange{Int}}))
    v = pdv[inds]
    deleteat!(pdv.refs, inds)
    v
end

function Base.splice!(pdv::PooledDataVector, inds::(@compat Union{Integer, UnitRange{Int}}), ins::AbstractVector)
    v = pdv[inds]
    splice!(pdv.refs, inds, [getpoolidx(pdv, v) for v in ins])
    v
end

Base.deleteat!(pdv::PooledDataVector, inds) = (deleteat!(pdv.refs, inds); pdv)

function Base.append!(da::AbstractDataVector, items::AbstractVector)
    oldn = length(da)
    itn = length(items)
    resize!(da, oldn+itn)
    da[oldn+1:end] = items[1:itn]
    da
end

function Base.sizehint!(da::DataVector, newsz::Integer)
    sizehint!(da.data, newsz)
    sizehint!(da.na, newsz)
end
Base.sizehint!(pda::PooledDataVector, newsz::Integer) =
    sizehint!(pda.refs, newsz)

# Pad a vector with NA's

function padNA(dv::AbstractDataVector,
               front::Integer,
               back::Integer)
    n = length(dv)
    res = similar(dv, front + n + back)
    for i in 1:n
        res[i + front] = dv[i]
    end
    return res
end


text:
using ConvexBodyProximityQueries
using CurveProximityQueries
using Plots
import Random: seed!

pgfplots()

seed!(22)
poly = randpoly([-0.5, 0.4]; n=20, scale=0.5)

curve = Bernstein([[0.944223, 0.177264], [0.46193, 0.105887], [0.791169, 0.172386], [0.256546, 0.463735], [0.70966, 0.3173], [0.664713, 0.830956], [0.952688, 0.736242], [0.00763912, 0.97189], [0.270595, 0.827442], [0.16636, 0.996615], [0.381073, 0.543077], [0.225381, 0.754049], [0.570043, 0.638543], [0.254696, 0.988563]])

pts = closest_points(poly, curve)

plot(poly,
     fillalpha=1.0,
     fillcolor=Gray(0.7),
    )
plot!(curve,
      linewidth=2.0,
      linecolor=:black,
     )
plot!(pts,
      linestyle=:dash,
      markersize=8.0,
      linewidth=2.0,
     )
savefig("out/fig1a.pdf")


text:
function abstract_word_constructors_test(::Type{Wo}) where Wo

    @testset "constructors: $Wo" begin

        @test Wo([1,2]) isa Wo
        @test Wo([1,2]) isa KnuthBendix.AbstractWord
        @test Wo{Int}([1,2]) isa Wo{Int}

        @test one(Wo{Int}) isa KnuthBendix.AbstractWord{Int}
        @test one(Wo{Int}) isa Wo{Int}

        @test similar(Wo([1,2])) isa Wo

        @test Word(Wo([1,2])) isa Word
        @test Wo(Word([1,2])) isa Wo
    end
end

function abstract_word_basic_functions_test(::Type{Wo}) where Wo

    @testset "basic functionality: $Wo" begin
        @test one(Wo{Int}) isa Wo
        @test one(Wo{UInt16}) isa Wo

        w = Wo([1,2])
        W = Word{Int}([1,2])

        @test W  == w
        @test W !== w
        @test one(w) == one(W)
        @test !isone(w)
        @test !isone(W)
        @test isone(one(w))
        @test isone(one(W))

        @test hash(w) isa UInt
        @test hash(w) == hash(W)
        @test hash(w, UInt(1)) != hash(w, UInt(0))
        @test hash(w) != hash(one(w))

        @test length(Set(Any[w,W])) == 1

        @test deepcopy(w) == w
        @test deepcopy(w) !== w
    end
end

function abstract_word_push_pop_append_test(::Type{Wo}) where Wo

    @testset "push!,pop!,append!,prepend!: $Wo" begin
        w = one(Wo{UInt16}); W = deepcopy(w);

        @test push!(w, 1) == Wo([1])
        @test w == Wo([1]) && isone(W)
        W = deepcopy(w)

        @test append!(w, [1,2]) == Wo([1,1,2])
        @test w == Wo([1,1,2]) && W == Wo([1])
        W = deepcopy(w)

        @test pop!(w) == 2
        @test w == Wo([1,1]) && W == Wo([1,1,2])
        W = deepcopy(w)

        @test pushfirst!(w, 3) == Wo([3,1,1])
        @test w == Wo([3,1,1]) && W == Wo([1,1])
        W == deepcopy(w)

        @test popfirst!(w) == 3
        @test w == Wo([1,1]) && W == Wo([1,1])
        W == deepcopy(w)

        @test prepend!(w,[1,2]) == Wo([1,2,1,1])
        @test w == Wo([1,2,1,1]) && W == Wo([1,1])
        W = deepcopy(w)

        @test resize!(w, 3) == Wo([1,2,1])
        @test w == Wo([1,2,1]) && W == Wo([1,2,1,1])
        W = deepcopy(w)

        @test resize!(w, 30) isa Wo{eltype(w)}
        @test length(w) == 30 && length(W) == 3
        @test resize!(w, 3) == W
        @test w == W

        @test_throws ArgumentError resize!(w, -4)
    end
end

function abstract_word_indexing_test(::Type{Wo}) where Wo

    @testset "indexing: $Wo" begin
        w = Wo([1,2]); W = deepcopy(w);

        @test collect(W) isa Vector{eltype(w)}
        @test collect(push!(W, 3)) == push!(collect(w), 3)
        @test W[end] == 3

        WW = deepcopy(W)

        @test collect(pushfirst!(W, 4)) == pushfirst!(collect(WW), 4)
        @test W[1] == 4

        @test W[2:3] isa Wo
        @test W[1:1] == Wo([4])
        @test W[end:end] isa Wo
        @test W[end:end] == Wo([3])
        @test W[2:2][1] == W[2]

        W_arr = collect(W)
        lw = length(W_arr)

        @test append!(W, Wo([6])) == Wo([W_arr; 6])
        @test prepend!(W, Wo([5])) == Wo([5; W_arr; 6])
        @test collect(W) isa Vector{eltype(w)}
        @test collect(W) == [5; W_arr; 6]
        @test W[1] == 5
        @test length(W) == lw + 2
        @test_throws BoundsError W[-1]
        @test_throws BoundsError W[lw+3]

        @test @view(W[2:3]) isa KnuthBendix.AbstractWord
    end
end

function abstract_word_arithmetic_test(::Type{Wo}) where Wo

    @testset "arithmetic: $Wo" begin

        @test Wo([1, 2]) * Wo([2, 3]) == Wo([1, 2, 2, 3])
        @test Wo([1, 2]) * Wo{Int}([2, 3]) == Wo{Int}([1, 2, 2, 3])
        @test Wo{Int}([1, 2]) * Wo([1, 2]) == Wo{Int}([1, 2, 1, 2])

        w = Wo([1,2])
        @test w^2 == Wo([collect(w); collect(w)])
        w_arr = collect(w)

        u1 = Word([1,2,3,4])
        u2 = Word([1,2])
        u3 = Word([4,1,2,3,4])
        u4 = Word([1])
        u5 = Word([2,3,4])

        @test KnuthBendix.longestcommonprefix(u1, u2) == 2
        @test KnuthBendix.longestcommonprefix(u1, u1) == 4
        @test KnuthBendix.lcp(u3, u2) == 0

        @test  KnuthBendix.isprefix(u2, u1)
        @test !KnuthBendix.isprefix(u2, u3)
        @test  KnuthBendix.issuffix(u5, u3)
        @test !KnuthBendix.issuffix(u5, u2)

        @test occursin(u2, u1) == true
        @test occursin(u2, u3) == true
        @test occursin(u1, u2) == false
        @test occursin(u4, u2) == true

        @test findnext(isequal(2), u3, 1)  == 3
        @test findnext(isequal(2), u3, 4)  === nothing
    end
end

function abstract_word_conformance_test(::Type{Wo}) where Wo
    @testset "KnuthBendix.AbstractWord conformance test: $Wo" begin
        abstract_word_constructors_test(Wo)
        abstract_word_push_pop_append_test(Wo)
        abstract_word_basic_functions_test(Wo)
        abstract_word_indexing_test(Wo)
        abstract_word_arithmetic_test(Wo)
    end
end


text:
export RK4

"""
    RK4()

The classical [Runge-Kutta](https://en.wikipedia.org/wiki/Runge–Kutta_methods)
predictor of order 4.
"""
struct RK4 <: AbstractPredictor end
struct RK4Cache{T} <: AbstractPredictorCache
    dt::Vector{T}
    k1::Vector{T}
    k2::Vector{T}
    k3::Vector{T}
    k4::Vector{T}
end

function cache(::RK4, H, x, ẋ, t)
    RK4Cache(dt(H, x, t), copy(ẋ), copy(ẋ), copy(ẋ), copy(ẋ))
end
#
function predict!(
    xnext,
    cache::RK4Cache,
    H::HomotopyWithCache,
    x,
    t,
    Δt,
    ẋ,
    Jac::JacobianMonitor,
)
    dt, mk₁, mk₂, mk₃, mk₄ = cache.dt, cache.k1, cache.k2, cache.k3, cache.k4
    n = length(xnext)
    @inbounds for i = 1:n
        mk₁[i] = -ẋ[i]
        xnext[i] = x[i] - 0.5Δt * mk₁[i]
    end

    minus_ẋ!(mk₂, H, xnext, t + 0.5Δt, Jac, dt)
    @inbounds for i = 1:n
        xnext[i] = x[i] - 0.5Δt * mk₂[i]
    end

    minus_ẋ!(mk₃, H, xnext, t + 0.5Δt, Jac, dt)
    @inbounds for i = 1:n
        xnext[i] = x[i] - Δt * mk₃[i]
    end

    minus_ẋ!(mk₄, H, xnext, t + Δt, Jac, dt)
    @inbounds for i = 1:n
        xnext[i] = x[i] -
                   0.16666666666666666Δt * (mk₁[i] + 2 * mk₂[i] + 2 * mk₃[i] + mk₄[i])
    end

    nothing
end

order(::RK4Cache) = 5


text:
include("../src/OnePhase.jl")

#using JLD, OnePhase, NLPModels, CUTEst
using JLD, NLPModels, CUTEst
USE_IPOPT = false
if USE_IPOPT
  using Ipopt
end
include("ipopt.jl")
include("plots.jl")
include("summary.jl")
include("process_results.jl")
include("run_and_store.jl")

function if_mkdir(dir::String)
  if !isdir(dir)
     mkdir(dir)
  end
end;

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
import Base: +, -, *, /, ^, ==, ≈
import Base.abs

struct Interval
    left::Real
    right::Real

    function Interval(num1::Real, num2::Real)
        if num1 > num2
            error("Invalid interval [$num1, $num2], left must be less than right")
        end
        new(Float64(num1), Float64(num2))
    end

    function Interval(a::Real)
        new(Float64(a), Float64(a))
    end

    function Interval(a::Vector{<:Real})
        new(Float64(a[1]), Float64(a[2]))
    end

    function Interval()
        new(NaN, NaN)
    end
end

..(a::Real, b::Real) = Interval(a, b)

function Base.getindex(a::Interval, index::Int64)
    if index == 1 return a.left
    elseif index == 2 return a.right
    else return nothing
    end
end

Base.iterate(a::Interval) = a[1], 2
Base.iterate(a::Interval, state) = state > 2 ? nothing : (a[state], state+1)
Base.vec(a::Interval) = [a.left, a.right]

Base.:(==)(a::Interval, b::Interval) = a.left == b.left && a.right == b.right
Base.:≈(a::Interval, b::Interval) = a.left ≈ b.left && a.right ≈ b.right

Base.isempty(a::Interval) = isnan(a.left) || isnan(a.right)

Base.copy(a::Interval) = Interval(a.left, a.right)

width(a::Interval) = a.right - a.left
mid(a::Interval) = (a.left + a.right) / 2
rad(a::Interval) = (a.right - a.left) / 2
mag(a::Interval) = max(abs(a.left), abs(a.right))
mig(a::Interval) = min(abs(a.left), abs(a.right))
abs(a::Interval) = Interval(mig(a), mag(a))
sqr(a::Interval) = a ^ 2

# INTERVAL ARITHMETIC ========================================================

Base.:-(a::Interval) = Interval(-a.right, -a.left)

Base.:+(a::Real, b::Interval) = Interval(a + b.left, a + b.right)
Base.:+(a::Interval, b::Real) = b + a
Base.:+(a::Interval, b::Interval) = Interval(a.left + b.left, a.right + b.right)
Base.:-(a::Interval, b::Interval) = Interval(a.left - b.right, a.right - b.left)

function Base.:*(a::Interval, b::Interval)
    m1 = a.left * b.left
    m2 = a.left * b.right
    m3 = a.right * b.left
    m4 = a.right * b.right
    left = min(m1, m2, m3, m4)
    right = max(m1, m2, m3, m4)
	Interval(left, right)
end

function Base.:*(a::Real, b::Interval)
    left = a * b.left
    right = a * b.right
    if a < 0
        left, right = right, left
    end
	Interval(left, right)
end
Base.:*(a::Interval, b::Real) = b * a

function Base.inv(a::Interval)
    if a.left == 0
        left = a.left
        right = 1 / a.right
    elseif a.left < 0 && 0 < a.right
        left = 1 / a.left
        right = 1 / a.right
    else
        left = a.right == 0 ? 0 : 1 / a.right
        right = a.left == 0 ? 0 : 1 / a.left
    end
    # left = 1 / a.right
    # right = 1 / a.left
    # println(a.left, ", ", a.right)
    # println(left, ", ", right)
    Interval(left, right)
end
Base.:/(a::Real, b::Interval) = a * inv(b)
Base.:/(a::Interval, b::Interval) = a * (1 / b)
Base.:/(a::Interval, b::Real) = a / Interval(b)

function Base.:^(a::Interval, b::Real)
    # println(a.left, ",", a.right)
    if b < 0
        a = 1 / a
        b = -b
    end

    if b == 0.5
        left = sqrt(a.left)
        right = sqrt(a.right)
    elseif b == 2
        if a.left <= 0 && 0 <= a.right
            left = 0
            right = mag(a)^2
        elseif a.left^2 <= a.right^2
            left = a.left ^ b
            right = a.right ^ b
        elseif a.right^2 <= a.left^2
            left = a.right ^ b
            right = a.left ^ b
        end
    else
        is_integer = floor(Int, b) == b
        b = is_integer ? floor(Int, b) : b
        # if (a.left > 0) || isodd(b)
        #     left = a.left ^ b
        #     right = a.right ^ b
        if (a.right < 0) && iseven(b)
            left = a.right ^ b
            right = a.left ^ b
        elseif (a.left <= 0 && 0 <= a.right) && iseven(b)
            left = 0
            right = mag(a) ^ b
        else
            left = a.left ^ b
            right = a.right ^ b
        end
    end
    left = left == Inf ? 0 : left
    right = right == Inf ? 0 : right
	Interval(left, right)
end

function Base.round(a::Interval; digits=2)
    Interval(round(a.left, digits=digits), round(a.right, digits=digits))
end


text:
struct MyInterp{T,N,A<:AbstractArray{T,N}}
    data::A
end

@inline (itp::MyInterp{T,N})(i::Vararg{<:Number,N}) where {T,N} = expand(itp, i)

@inline function expand(itp::MyInterp, ipre::Tuple{Vararg{Number,L}}, ipost::Vararg{Integer,M}) where {L,M}  # force specialization
    ifront, ilast = Base.front(ipre), ipre[end]
    im, ip = floor(ilast), ceil(ilast)
    return (ip - ilast)*expand(itp, ifront, unsafe_trunc(Int, im), ipost...) +
           (ilast - im)*expand(itp, ifront, unsafe_trunc(Int, ip), ipost...)
end

@inline expand(itp::MyInterp, ::Tuple{}, ipost::Vararg{Integer,N}) where N =
    @inbounds itp.data[CartesianIndex(ipost)]


text:
using Dates, LibGit2
import Base: download, show

# function Base.show(io::IO, x::Union{Float64,Float32})
#     Base.Grisu._show(io, round(x, sigdigits=4), Base.Grisu.SHORTEST, 0,
#                      get(io, :typeinfo, Any) !== typeof(x), false)
# end

Counter()  = zeros(Int,2)
Accuracy() = (lemma=Counter(), tag=Counter(), complete=Counter())
F1Scores() = (score=Dict{String,Float64}(), total=Dict{String,Float64}())
F1Metric() = (precision=F1Scores(), recall=F1Scores())

"""
    printLog(f::IO, str...)

    Logger for Morse. Prints current time in eachline and flush the stream
"""
printLog(f::IO, str...) = (println(f,Dates.Time(now()),": ",str...); flush(f);)

"""
    printConfig(f::IO,o::Dict{Symbol,Any})

    Prints the configuration dictionary
"""
function printConfig(f::IO,o::Dict{Symbol,Any})
    printLog(f,"Configuration: ")
    for (k,v) in o; v!== nothing && println(f, k, " => " , v); end
    flush(f)
end

"""
    download(dataset::Type{<:DataSet})

    downloads `DataSet` files.
"""
function download(dataset::Type{UDDataSet})
    if !isdir(dir("data","ud-treebanks-v2.1"))
        download("https://lindat.mff.cuni.cz/repository/xmlui/bitstream/handle/11234/1-2515/ud-treebanks-v2.1.tgz",dir("data","ud-treebanks-v2.1.tgz"))
        run(`tar -xf $(dir("data","ud-treebanks-v2.1.tgz")) -C $(dir("data"))`)
        rm(dir("data","ud-treebanks-v2.1.tgz"))
    end
end

function download(dataset::Type{TRDataSet}; path=dir("data","TrMor2018"))
    if !isdir(path)
        repo = LibGit2.clone("https://github.com/ai-ku/TrMor2018", path)
        LibGit2.checkout!(repo,"2b50593e756b36ba7fc3f8b16a9f9f92589710fc")
        run(`ln -s  $(joinpath(path,"TrMor2006","trmor2006.train")) $(joinpath(path,"TrMor2016","trmor2016.train"))`)
    end
end

const server_url ="people.csail.mit.edu/deniz/models/morse/"

function download(model::Type{T}, format; vers="2018", lemma=true, lang="tr") where T
    flang = format===TRDataSet ? string("TR-tr",vers) : string("UD-",lang)
    mname = string("bestModel.",T,"_lemma_",lemma,"_lang_",flang,"_size_full","_params.jld2")
    lpath = dir("checkpoints",mname)
    if !isfile(lpath)
        mpath = string(server_url, mname)
        download(mpath,lpath)
    end
    lpath
end

function trained(x...;o...)
    lpath = download(x...;o...)
    model,_,vocabulary,parser = loadModel(lpath)
    return model,vocabulary,parser
end

"""
    loadModel(fname::AbstractString)

    Loads from saved model file. Return tuple of (model,config,vocab,parser)
"""
function loadModel(fname::AbstractString)
    f = KnetLayers.load(fname)
    prms,opts,vocab,parser = f["model"], f["opts"], f["vocab"], f["parser"]
    ModelType = eval(Meta.parse(opts[:modelType]))
    model = ModelType(opts,vocab)
    for (wm,wl) in zip(params(model),prms)
        copyto!(wm.value,wl)
    end
    return model,opts,vocab,parser
end

"""
    saveModel(fname::AbstractString, model, config, vocab, parser)

    Saves model, config, vocab, parser to a file.
"""
saveModel(fname::AbstractString, model::Model, config::Dict, v::Vocabulary, p::Parser) =
    KnetLayers.save(fname, "model", model, "opts", config, "vocab", v, "parser", p)

"""
   setoptim!(M::Model, optimizer)

   Set optimizer for model parameters
"""
setoptim!(M::Model, optimizer) =
    for p in params(M); p.opt = deepcopy(optimizer); end

"""
   lrdecay!(M::Model, decay::Real)

   Decays learning rate
"""
lrdecay!(M::Model, decay::Real) =
    for p in params(M); p.opt.lr = p.opt.lr*decay; end

"""
   splitdata(data::Vector, r::Vector{<:Real})

   Splits an array according to r.
   r can be Integer array in which every element specifies the exact size of the splits
   r can be Float array in which every element specifies the ratio of the splits
"""
function splitdata(data::Vector, r::Vector{<:Real})
    ls = isa(r[1],AbstractFloat) ? floor.(Int, r .* length(data)) : r #Ratio or Number
    map(rng->data[rng[1]:rng[2]], zip(cumsum(ls) .- ls .+ 1, cumsum(ls)))
end

"""
    createsplits(data::Vector, seed::Integer, tsize::Integer)

    If there is 1 or 2 splits makes it 3 (trn,dev,test)
    according to `seed` and `tsize` inputs.
"""
function createsplits(data::Vector,seed::Integer,tsize::Integer)
    l  = length(data)
    tl = length(data[1])
    shuffle!(MersenneTwister(seed),data[1])
    if l==1
        if tsize < tl
            testsize = min((tl-tsize) ÷ 2, floor(Int,tl*0.1))
            data = splitdata(data[1],[tsize,testsize,testsize])
        else
            data = splitdata(data[1],[0.8,0.1,0.1])
        end
    elseif l==2
        if tsize < tl
            testsize = min(tl-tsize, floor(Int,tl*0.1), 2.5length(data[2]))
            data = push!(splitdata(data[1],[tsize,testsize]),data[2])
        else
            devratio = min(floor(Int,length(data[1])*0.1), 2.5length(data[2])) / length(data[1])
            data = push!(splitdata(data[1],[0.9,devratio]),data[2])
        end
    else
        if tsize < tl
            data[1] = data[1][1:tsize]
        end
    end
    return data
end

"""
   prepareData(files::Vector{<:AbstractString}, dtype::Type{<:DataSet};
                           seed=31, tsize=typemax(Int), withLemma=true)

   preparesData for training and evaluation for datafiles.
   returns encodedData, vocabulary, parser
"""
function prepareData(files::Vector{<:AbstractString}, dtype::Type{<:DataSet};
                     vers=2018, seed=31, tsize=typemax(Int), withLemma=true, parseAll=false)
    parser = Parser{dtype}(Val(vers))
    data = createsplits(parseFile.(files; p=parser, withLemma=withLemma, parseAll=parseAll), seed, tsize)
    vocab = Vocabulary(data)
    return encode.(data, v=vocab), vocab, parser
end

"""
    StringAnalysis

    Keeps analysis of a word with String fields.
"""
struct StringAnalysis
    word::String
    lemma::String
    tags::Vector{String}
end

# Sequential
function StringAnalysis(word::Vector{Int}, seq::Vector{Int}; v::Vocabulary, isgold=false)
    lemma, tags = String[], String[]
    for index in seq
        if index == v.specialIndices.eow
            push!(tags,v.specialTokens.eow); break
        end
        !isgold && index ∈ v.specialIndices && continue
        str = v.tags[index]
        length(str)==1 ? push!(lemma,str) : push!(tags,str)
    end
    StringAnalysis(join(v.chars[word]), join(lemma), tags)
end

# Discriminative
function StringAnalysis(word::Vector{Int}, output::Int; v::Vocabulary, isgold=false)
    StringAnalysis(join(v.chars[word]), "",
                   push!(split(v.comptags[output],'+'),v.specialTokens.eow))
end

# Vocabulary is needed for digit masking
function getLabels(pred::StringAnalysis, gold::StringAnalysis; v::Vocabulary)
    if length(gold.lemma) != length(pred.lemma)
        return (false, pred.tags == gold.tags)
    end
    for (gc,pc) in zip(gold.lemma, pred.lemma)
        isdigit(gc) && continue
        gc != pc && return (false,pred.tags == gold.tags)
    end
    return (true, pred.tags == gold.tags)
end

function makeFormat(a::StringAnalysis, p::Parser{UDDataSet})
    lemma = length(a.lemma) == 0 ? "X" : join(a.lemma)
    tags = a.tags[1:end-1]
    posTag = length(tags) == 0 ? "X" : tags[1]
    morphFeats = length(tags) < 2 ? p.unkToken : join(tags[2:end],p.tagsSeperator)
    return (a.word, lemma, posTag, morphFeats)
end

printFormat(predS,parser; f=stdout) =
    for (i,pred) in enumerate(predS) printFormat(f,i,pred,parser) end

function printFormat(f::IO, i::Integer, predStr::NTuple, p::Parser{UDDataSet})
    word, lemma, pos_tag, morph_feats = predStr
    s,u = p.partsSeperator, p.unkToken
    if startswith(pos_tag,"Upos=")
        pos_tag = SubString(pos_tag,6:length(pos_tag))
    end
    write(f,string(i),s,word,s,lemma,s,pos_tag,s,u,s,morph_feats,s,u,s,u,s,u,s,u,'\n')
end

function makeFormat(a::StringAnalysis, p::Parser{TRDataSet})
    lemma  = length(a.lemma) == 0 ? "" : join(a.lemma)
    tags   = a.tags[1:end-1]
    morphFeats = join(tags,p.tagsSeperator)
    return (a.word, lemma, "", morphFeats)
end

function printFormat(f::IO, i::Integer, predStr::NTuple, p::Parser{TRDataSet})
    word, lemma, pos_tag, morph_feats = predStr
    write(f,word,p.partsSeperator,lemma,p.tagsSeperator,morph_feats,'\n')
end

function evaluate(M::Model, data::Vector{SentenceBatch}, v::Vocabulary, p::Parser{T}; file=nothing) where T
    # Initialize metrics, amb: accuracy metric in ambigious words
    loss, acc, amb, F1 = zeros(2), Accuracy(), Accuracy(), F1Metric()

    for d in progress(data)
        # Forward Run
        preds, J = predict(M, d; v=v)
        # Cross Entropy Loss Metric
        if isa(M,Sequential)
            loss .+= [J, sum(d.masks)]
        else
            loss .+= [J, 1] .* length(d.encodedIOs)
        end

        if file !== nothing
            write(file, p.sentenceStart, " = ",
                  join(map(a->join(v.chars[a.chars]), d.encodedIOs)," "),
                  '\n')
        end

        @inbounds for (i,encodedIO) in enumerate(d.encodedIOs)

            analysis = encodedIO.analyses[1]  #correct analysis

            if isa(M,Sequential)
                pred  = StringAnalysis(encodedIO.chars, preds[i,:], v=v)
                gold  = StringAnalysis(encodedIO.chars, d.seqOutputs[i,:], v=v, isgold=true)
            else
                pred  = StringAnalysis(encodedIO.chars, preds[i], v=v)
                gold  = StringAnalysis(encodedIO.chars, d.compositeOutputs[i], v=v, isgold=true)
            end

            tf = getLabels(gold, pred, v=v) #(lemma::Bool, tag::Bool, complete::Bool)

            if analysis.isValid

                F1update!(pred.tags, gold.tags, F1.precision, F1.recall, p=p)

                for (accuracy, mul) in zip((acc,amb), (true,encodedIO.isAmbigous))
                    accuracy.lemma     .+= [Int(tf[1]), 1] .* mul
                    accuracy.tag       .+= [Int(tf[2]), 1] .* mul
                    accuracy.complete  .+= [Int(tf[1]&tf[2]), 1] .* mul
                end
            end
            #Generation
            file !== nothing && printFormat(file, i, makeFormat(pred, p), p)
        end
        T === UDDataSet && file !== nothing && print(file,'\n')
    end
   T === UDDataSet && file !== nothing && print(file,"# text = _end_")
   f1macro, f1micro = F1average(F1.precision, F1.recall)
   return (loss=loss[1]/loss[2], accuracies=percentage(acc), amb=percentage(amb),
            f1macro=f1macro, f1micro=f1micro)
end

percentage(a::NamedTuple) = (lemma=100a.lemma[1]/a.lemma[2], tag=100a.tag[1]/a.tag[2],
                             complete=100a.complete[1]/a.complete[2])

function F1average(precision::NamedTuple, recall::NamedTuple)
    micro_precision = sum(values(precision.score))/sum(values(precision.total))
    micro_recall = sum(values(recall.score))/sum(values(recall.total))
    micro_score = 2(micro_precision * micro_recall)/(micro_precision + micro_recall)

    for k in keys(precision.score)
        precision.score[k] = precision.score[k]/precision.total[k]
    end

    f1_scores = Dict{String,Float64}()
    f1_average = 0.0
    for k in keys(recall.score)
        recall.score[k] = recall.score[k]/recall.total[k]
        if recall.score[k]==0 || !haskey(precision.score,k)
            f1_scores[k] = 0
        else
            f1_scores[k] = 2(precision.score[k] * recall.score[k])/
                             (precision.score[k] + recall.score[k])
        end
        f1_average += recall.total[k] * f1_scores[k]
    end
    f1_average /= sum(values(recall.total))

    return 100(f1_average), 100(micro_score)
end

function F1update!(pred::Vector{String}, gold::Vector{String},
                   precision::NamedTuple, recall::NamedTuple; p::Parser{T}) where T
    if T === UDDataSet
        p_dict = Dict(filter(x->length(x)>1, split.(pred,p.tagValueSeperator)))
        g_dict = Dict(filter(x->length(x)>1, split.(gold,p.tagValueSeperator)))
    else # TRDataSet
        p_dict = Dict(map(x->(x,true),pred))
        g_dict = Dict(map(x->(x,true),gold))
    end
    for (x,y,score,total) in ( (p_dict, g_dict, precision.score, precision.total),
                               (g_dict, p_dict, recall.score, recall.total) )
        for (k,v) in x
            if !haskey(score,k)
                score[k] = total[k] = 0
            end
            if haskey(y,k) && y[k] == v
                score[k] += 1
            end
            total[k] += 1
        end
    end
end

#
# function generate(M::Model, data; v::Vocabulary, p::Parser, prefix="")
#     open(prefix*".generation","w") do f
#         for (i,d) in enumerate(data)
#             preds  = predict(M,d; v=v)
#             write(f, p.sentenceStart, " = ", join(words," "), '\n')
#             for (k,encodedIO) in d.encodedIOs
#                 a = encodedIO.analyses[1] # correct analysis
#                 formatted = makeFormat(StringAnalysis(a.chars, preds[k,:]; v=v), p)
#                 printFormat(f, k, formatted, p)
#             end
#             i != length(data) && write(f,'\n')
#         end
#     end
# end


text:
# Show some associative array math.

# Read CSV file into associative array, get vertices and convert to numbers.
E = ReadCSV(joinpath(Base.source_dir(),"Edge.csv"))
Ev = logical( E[:, StartsWith("V,")] )

# Get orange and green edges
EvO = Ev[StartsWith("O,"),:]
EvG = Ev[StartsWith("G,"),:]

# Compute (empty) vertex adjacency graph.
AvOG = transpose(EvO) * EvG
printFull(AvOG)

# Compute edge adjacency graph.
AeOG = EvO * transpose(EvG)
printFull(AeOG)

# Compute edge adjacency graph preserving keys.
AeOG = CatKeyMul(EvO,transpose(EvG))
printFull(AeOG)


text:
# produce a wireframe using Plots
import Plots: plot

mutable struct PlotData
  title :: AbstractString
  x :: Array{Float64,1}
  y :: Array{Float64,1}
  z :: Array{Float64,1}
  group :: Array{Int,1}
  marker :: Array{Symbol,1}
  markercolor :: Array{Symbol,1}
  markeralpha :: Array{Float64,1}
  markersize :: Array{Float64,1}
  markerstrokewidth :: Array{Float64,1}
  linecolor :: Array{Symbol,1}
  groupcounter :: Int
  PlotData() = new("",[],[],[],[],[],[],[],[],[],[],0)
end

struct PlotElementParameters
  marker :: Symbol
  markercolor :: Symbol
  markeralpha :: Float64
  markersize :: Float64
  markerstrokewidth :: Float64
  linecolor :: Symbol
end
struct PlotScheme
  node :: PlotElementParameters
  segment :: PlotElementParameters
  plane :: PlotElementParameters
  planenode :: PlotElementParameters
  invisiblenode :: PlotElementParameters
end
const plotelementparameters_node=
  PlotElementParameters(:circle,:red,0.3,3.0,0.1, :blue)
const plotelementparameters_segment=
  PlotElementParameters(:none,  :red,0.3,3.0,0.1, :blue)
const plotelementparameters_plane=
  PlotElementParameters(:none,  :red,0.3,3.0,0.1, :green)
const plotelementparameters_planenode=
  PlotElementParameters(:circle,:green,0.3,3.0,0.1, :green)
const plotelementparameters_invisiblenode=
  PlotElementParameters(:cross, :red,0.0,0.0,0.0, :blue)
const defaultplotscheme = PlotScheme(plotelementparameters_node,
                              plotelementparameters_segment,
                              plotelementparameters_plane,
                              plotelementparameters_planenode,
                              plotelementparameters_invisiblenode)

function Base.push!(pd::PlotData, xyz::Array{Float64,1})
  push!(pd.x, xyz[1])
  push!(pd.y, xyz[2])
  push!(pd.z, xyz[3])
  return nothing
end

function Base.push!(pd::PlotData, pep::PlotElementParameters)
  push!(pd.marker,            pep.marker)
  push!(pd.markercolor,       pep.markercolor)
  push!(pd.markeralpha,       pep.markeralpha)
  push!(pd.markersize,        pep.markersize)
  push!(pd.markerstrokewidth, pep.markerstrokewidth)
  push!(pd.linecolor,         pep.linecolor)
  return nothing
end

function pushgroupcounter!(pd::PlotData; incerment::Int=0)
  pd.groupcounter +=  incerment
  push!(pd.group,pd.groupcounter)
  return nothing
end

# push! a group with a single marker into PlotData
function pushmarkergroup!(pd::PlotData,
                          xyz::Array{Float64,1},
                          pep::PlotElementParameters)
  pushgroupcounter!(pd, incerment=1)
  push!(pd,xyz)
  push!(pd,pep)
  return nothing
end

# force same scale on all axis with invisible markers in corners
function pointsatlimits!(pd::PlotData, ps::PlotScheme = defaultplotscheme)
  xlims = [minimum(pd.x) maximum(pd.x)]
  ylims = [minimum(pd.y) maximum(pd.y)]
  zlims = [minimum(pd.z) maximum(pd.z)]
  xcenter = (xlims[1]+xlims[2])/2.0
  ycenter = (ylims[1]+ylims[2])/2.0
  zcenter = (zlims[1]+zlims[2])/2.0
  range = maximum([xlims[2]-xlims[1] ylims[2]-ylims[1] zlims[2]-zlims[1]])
  pushmarkergroup!(pd,
                  [xcenter+range/2.0,ycenter+range/2.0,zcenter+range/2.0],
                  ps.invisiblenode)
  pushmarkergroup!(pd,
                  [xcenter-range/2.0,ycenter-range/2.0,zcenter-range/2.0],
                  ps.invisiblenode)
end

Base.transpose(x::Array{Symbol,1}) = reshape(x,1,length(x))
function plot(e::Element, ps::PlotScheme = defaultplotscheme)
  pd = plotdata(e,ps)
  pointsatlimits!(pd,ps)
  plot(pd.x, pd.y, pd.z, group=pd.group,
    title = pd.title,
    legend = false,
    linecolor = transpose(pd.linecolor),
    marker=transpose(pd.marker),
    markercolor = transpose(pd.markercolor),
    markeralpha = transpose(pd.markeralpha),
    markersize = transpose(pd.markersize),
    markerstrokewidth = transpose(pd.markerstrokewidth))
end

function plotdata(element::Element,ps::PlotScheme = defaultplotscheme)
  pd = PlotData()
  pd.groupcounter = 0
  context = Context(element)
  pd.title = context.title
  for e in traverseTree(element)
    appendplotdata!(pd,e,context,ps)
  end
  return pd
end

appendplotdata!(::PlotData, ::Element, ::Context, ::PlotScheme) =
  nothing
function appendplotdata!(pd::PlotData,
                         node::Node,
                         context::Context,
                         ps::PlotScheme)
  pushmarkergroup!(pd, xyz1(node,context), ps.node)
  return nothing
end
function appendplotdata!(pd::PlotData,
                         segment::Segment,
                         context::Context,
                         ps::PlotScheme)
  appendplotdata!(pd, corners(segment,context), ps.segment)
  return nothing
end
function appendplotdata!(pd::PlotData,
                         plane::UniformPlane,
                         context::Context,
                         ps::PlotScheme)
  appendplotdata!(pd, corners(plane,context), ps.plane)
  pep = ps.planenode
  for xyz1 in nodes_xyz1(plane,context)
    pushmarkergroup!(pd, xyz1, pep)
  end
  return nothing
end
function appendplotdata!(pd::PlotData,
                         corners::Array{Float64,2},
                         pep::PlotElementParameters)
  t = [(1,2,3,4,1,5,6,7,8,5),(2,6),(3,7),(4,8)]
  for g in t
    pd.groupcounter += 1
    push!(pd,pep)
    for p in g
      push!(pd.group,pd.groupcounter)
      push!(pd,corners[:,p])
    end
  end
  return nothing
end

function corners(segment::Segment, context::Context)
  (w,h) = width_height(segment,context)
  (n1_xyz1,n2_xyz1) = nodes_xyz1(segment,context)
  n1_xyz = n1_xyz1[1:3]
  n2_xyz = n2_xyz1[1:3]
  widthvector = wxyz(segment,context)
  normalize!(widthvector)
  halfwidthvector = 0.5*w * widthvector
  mp1 = n1_xyz+halfwidthvector
  mp2 = n1_xyz-halfwidthvector
  mp3 = n2_xyz+halfwidthvector
  mp4 = n2_xyz-halfwidthvector
  lengthvector = n2_xyz-n1_xyz
  normalize!(lengthvector)
  heightvector = cross(widthvector,lengthvector)
  normalize!(heightvector)
  halfheightvector = 0.5*h * heightvector
  c = Array{Float64}(undef,(3,8))
  c[:,1] = mp1+halfheightvector
  c[:,2] = mp1-halfheightvector
  c[:,3] = mp2-halfheightvector
  c[:,4] = mp2+halfheightvector
  c[:,5] = mp3+halfheightvector
  c[:,6] = mp3-halfheightvector
  c[:,7] = mp4-halfheightvector
  c[:,8] = mp4+halfheightvector
  return c
end
function corners(plane::UniformPlane, context::Context)
  (c1,c2,c3,thick) = corners_xyz1_thick(plane,context)
  widthvector = c3[1:3] - c2[1:3]
  lengthvector = c1[1:3] - c2[1:3]
  perpvector = cross(widthvector,lengthvector)
  normalize!(perpvector)
  halfthickperpvector = 0.5*thick *perpvector
  c = Array{Float64}(undef,(3,8))
  c4 = widthvector + c1[1:3]
  c[:,1] = c1[1:3] - halfthickperpvector
  c[:,2] = c2[1:3] - halfthickperpvector
  c[:,3] = c3[1:3] - halfthickperpvector
  c[:,4] = c4 - halfthickperpvector
  c[:,5] = c1[1:3] + halfthickperpvector
  c[:,6] = c2[1:3] + halfthickperpvector
  c[:,7] = c3[1:3] + halfthickperpvector
  c[:,8] = c4 + halfthickperpvector
  return c
end


text:
function load_constrained(op, u₁loop, u₂loop, forprefetch = false)
    loopdeps = loopdependencies(op)
    dependsonu₁ = u₁loop ∈ loopdeps
    if u₂loop === Symbol("##undefined##")
        if forprefetch
            dependsonu₁ || return false
        end
        # unrolleddeps = [ u₁loop ]
    else
        dependsonu₂ = u₂loop ∈ loopdeps
        if forprefetch
            (dependsonu₁ & dependsonu₂) || return false
        end
        # unrolleddeps = [ u₁loop, u₂loop ]
    end
    unrolleddeps = Symbol[]
    dependsonu₁ && push!(unrolleddeps, u₁loop)
    dependsonu₂ && push!(unrolleddeps, u₂loop)
    any(opp -> isload(opp) && all(in(loopdependencies(opp)), unrolleddeps), parents(op))
end
function check_if_remfirst(ls, ua)
    usorig = ls.unrollspecification[]
    @unpack u₁, u₁loopsym, u₂loopsym, u₂max = ua
    u₁loop = getloop(ls, u₁loopsym)
    u₂loop = getloop(ls, u₂loopsym)
    if isstaticloop(u₁loop) && (usorig.u₁ != u₁)
        return true
    end
    if isstaticloop(u₂loop) && (usorig.u₂ != u₂max)
        return true
    end
    false
end
function sub_fmas(ls::LoopSet, op::Operation, ua::UnrollArgs)
    @unpack u₁, u₁loopsym, u₂loopsym, u₂max = ua
    !(load_constrained(op, u₁loopsym, u₂loopsym) || check_if_remfirst(ls, ua))
end

struct FalseCollection end
Base.getindex(::FalseCollection, i...) = false
function parent_unroll_status(op::Operation, u₁loop::Symbol, u₂loop::Symbol, ::Nothing)
    map(opp -> isunrolled_sym(opp, u₁loop), parents(op)), FalseCollection()
end
# function parent_unroll_status(op::Operation, u₁loop::Symbol, u₂loop::Symbol, ::Nothing)
#     vparents = parents(op);
#     parent_names = Vector{Symbol}(undef, length(vparents))
#     parents_u₁syms = Vector{Bool}(undef, length(vparents))
#     # parents_u₂syms = Vector{Bool}(undef, length(vparents))
#     for i ∈ eachindex(vparents)
#         parent_names[i], parents_u₁syms[i], _ = variable_name_and_unrolled(vparents[i], u₁loop, u₂loop, nothing)
#     end
#     parent_names, parents_u₁syms, FalseCollection()
# end
function parent_unroll_status(op::Operation, u₁loop::Symbol, u₂loop::Symbol, u₂iter::Int)
    vparents = parents(op);
    # parent_names = Vector{Symbol}(undef, length(vparents))
    parents_u₁syms = Vector{Bool}(undef, length(vparents))
    parents_u₂syms = Vector{Bool}(undef, length(vparents))
    for i ∈ eachindex(vparents)
        # parent_names[i], parents_u₁syms[i], parents_u₂syms[i] = variable_name_and_unrolled(vparents[i], u₁loop, u₂loop, u₂iter)
        parents_u₁syms[i], parents_u₂syms[i] = isunrolled_sym(vparents[i], u₁loop, u₂loop, u₂iter)
    end
    # parent_names, parents_u₁syms, parents_u₂syms
    parents_u₁syms, parents_u₂syms
end

# """
#     Requires a parents_op argument, because it may `===` parents(op), due to previous divergence, e.g. to handle unrolling.
# """
# function isreducingidentity!(q::Expr, op::Operation, parents_op::Vector{Operation}, U::Int, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol, suffix)
#     vparents = copy(parents_op) # don't mutate the original!
#     for (i,opp) ∈ enumerate(parents_op)
#         @show opp vectorized ∈ loopdependencies(opp), vectorized ∈ reducedchildren(opp) # must reduce
#         @show vectorized, loopdependencies(opp), reducedchildren(opp) # must reduce
#         if vectorized ∈ loopdependencies(opp) || vectorized ∈ reducedchildren(opp) # must reduce
#             loopdeps = [l for l ∈ loopdependencies(opp) if l !== vectorized]
#             @show opp
#             reductinstruct = reduction_to_scalar(instruction(opp))

#             reducedparent = Operation(
#                 opp.identifier, gensym(opp.variable), opp.elementbytes, Instruction(:LoopVectorization, reductinstruct), opp.node_type,
#                 loopdeps, opp.reduced_deps, opp.parents, opp.ref, opp.reduced_children
#             )
#             pname,   pu₁,  pu₂ = variable_name_and_unrolled(opp, u₁loop, u₂loop, suffix)
#             rpname, rpu₁, rpu₂ = variable_name_and_unrolled(reducedparent, u₁loop, u₂loop, suffix)
#             @assert pu₁ == rpu₁ && pu₂ == rpu₂
#             if rpu₁
#                 for u ∈ 0:U-1
#                     push!(q.args, Expr(:(=), Symbol(rpname,u), Expr(:call, lv(reductinstruct), Symbol(pname,u))))
#                 end
#             else
#                 push!(q.args, Expr(:(=), rpname, Expr(:call, lv(reductinstruct), pname)))
#             end
#             vparents[i] = reducedparent
#         end
#     end
#     vparents
# end

function add_loopvalue!(instrcall::Expr, loopval::Symbol, vectorized::Symbol, u::Int)
    if loopval === vectorized
        if isone(u)
            push!(instrcall.args, Expr(:call, lv(:valadd), VECTORWIDTHSYMBOL, _MMind(Expr(:call, lv(:staticp1), loopval))))
        else
            push!(instrcall.args, Expr(:call, lv(:valmuladd), VECTORWIDTHSYMBOL, u, _MMind(Expr(:call, lv(:staticp1), loopval))))
        end
    else
        push!(instrcall.args, Expr(:call, lv(:vadd), loopval, u + 1))
    end
end
function add_loopvalue!(instrcall::Expr, loopval, ua::UnrollArgs, u::Int)
    @unpack u₁loopsym, u₂loopsym, vectorized, suffix = ua
    if u > 0 && loopval === u₁loopsym #parentsunrolled[n]
        add_loopvalue!(instrcall, loopval, vectorized, u)
    elseif !isnothing(suffix) && suffix > 0 && loopval === u₂loopsym
        add_loopvalue!(instrcall, loopval, vectorized, suffix)
    elseif loopval === vectorized
        push!(instrcall.args, _MMind(Expr(:call, lv(:staticp1), loopval)))
    else
        push!(instrcall.args, Expr(:call, lv(:staticp1), loopval))
    end
end

function lower_compute!(
    q::Expr, op::Operation, ls::LoopSet, ua::UnrollArgs, mask::Union{Nothing,Symbol,Unsigned} = nothing,
)
    @unpack u₁, u₁loopsym, u₂loopsym, vectorized, suffix = ua
    var = name(op)
    instr = instruction(op)
    parents_op = parents(op)
    nparents = length(parents_op)
    mvar, u₁unrolledsym, u₂unrolledsym = variable_name_and_unrolled(op, u₁loopsym, u₂loopsym, suffix)
    opunrolled = u₁unrolledsym || u₁loopsym ∈ loopdependencies(op)
    # parent_names, parents_u₁syms, parents_u₂syms = parent_unroll_status(op, u₁loop, u₂loop, suffix)
    parents_u₁syms, parents_u₂syms = parent_unroll_status(op, u₁loopsym, u₂loopsym, suffix)
    tiledouterreduction = if isnothing(suffix)
        -1
    else
        suffix_ = Symbol(suffix, :_)
        isouterreduction(op)
    end
    # parenttiled = isunrolled_sym.(parents_op, tiled, unrolled)
    # if instr.instr === :identity && name(first(parents_op)) === var && isone(length(parents_op))
    #     if (u₁unrolledsym == first(parents_u₁syms)) && ((!isnothing(suffix)) == parents_u₂syms[1])
    #         return
    #     end
    # end
    # unrollsym = isunrolled_sym(op, unrolled)
    if !opunrolled && any(parents_u₁syms) # TODO: Clean up this mess, refactor the naming code, putting it in one place and have everywhere else use it for easy equivalence.
        parents_op = copy(parents_op) # don't mutate the original!
        for i ∈ eachindex(parents_u₁syms)
            parents_u₁syms[i] || continue
            parents_u₁syms[i] = false
            parentop = parents_op[i]
            i == tiledouterreduction && isconstant(parentop) && continue
            newparentop = Operation(
                parentop.identifier, gensym(parentop.variable), parentop.elementbytes, parentop.instruction, parentop.node_type,
                parentop.dependencies, parentop.reduced_deps, parentop.parents, parentop.ref, parentop.reduced_children
            )
            parentname = mangledvar(parentop)
            newparentname = mangledvar(newparentop)
            parents_op[i] = newparentop
            if parents_u₂syms[i]
                parentname = Symbol(parentname, suffix_)
                newparentname = Symbol(newparentname, suffix_)
            end
            if isconstant(newparentop)
                push!(q.args, Expr(:(=), newparentname, Symbol(parentname, 0)))
            else
                for u ∈ 0:u₁-1
                    push!(q.args, Expr(:(=), Symbol(newparentname, u), Symbol(parentname, u)))
                end
                reduce_expr!(q, newparentname, Instruction(reduction_to_single_vector(instruction(newparentop))), u₁)
                push!(q.args, Expr(:(=), newparentname, Symbol(newparentname, 0)))
            end
        end
    end
    # cache unroll and tiling check of parents
    # not broadcasted, because we use frequent checks of individual bools
    # making BitArrays inefficient.
    # parentsyms = [opp.variable for opp ∈ parents(op)]
    Uiter = opunrolled ? u₁ - 1 : 0
    isreduct = isreduction(op)
    if !isnothing(suffix) && isreduct# && (iszero(suffix) || (ls.unrollspecification[].u₂ - 1 == suffix))
        # instrfid = findfirst(isequal(instr.instr), (:vfmadd, :vfnmadd, :vfmsub, :vfnmsub))
        instrfid = findfirst(isequal(instr.instr), (:vfmadd_fast, :vfnmadd_fast, :vfmsub_fast, :vfnmsub_fast))
        # want to instcombine when parent load's deps are superset
        # also make sure opp is unrolled
        if !isnothing(instrfid) && (opunrolled && u₁ > 1) && sub_fmas(ls, op, ua)
            specific_fmas = Base.libllvm_version >= v"11.0.0" ? (:vfmadd, :vfnmadd, :vfmsub, :vfnmsub) : (:vfmadd231, :vfnmadd231, :vfmsub231, :vfnmsub231)
            # specific_fmas = (:vfmadd231, :vfnmadd231, :vfmsub231, :vfnmsub231)
            instr = Instruction(specific_fmas[instrfid])
        end
    end
    reduceddeps = reduceddependencies(op)
    vecinreduceddeps = isreduct && vectorized ∈ reduceddeps
    maskreduct = !isnothing(mask) && vecinreduceddeps #any(opp -> opp.variable === var, parents_op)
    # if vecinreduceddeps && vectorized ∉ loopdependencies(op) # screen parent opps for those needing a reduction to scalar
    #     # parents_op = reduce_vectorized_parents!(q, op, parents_op, U, u₁loopsym, u₂loopsym, vectorized, suffix)
    #     isreducingidentity!(q, op, parents_op, U, u₁loopsym, u₂loopsym, vectorized, suffix) && return
    # end
    # if a parent is not unrolled, the compiler should handle broadcasting CSE.
    # because unrolled/tiled parents result in an unrolled/tiled dependendency,
    # we handle both the tiled and untiled case here.
    # bajillion branches that go the same way on each iteration
    # but smaller function is probably worthwhile. Compiler could theoreically split anyway
    # but I suspect that the branches are so cheap compared to the cost of everything else going on
    # that smaller size is more advantageous.
    isvectorized = vectorized ∈ loopdependencies(op)
    modsuffix = 0
    for u ∈ 0:Uiter
        instrcall = callexpr(instr)
        varsym = if tiledouterreduction > 0 # then suffix !== nothing
            modsuffix = ((u + suffix*(Uiter + 1)) & 3)
            # modsuffix = u
            # modsuffix = suffix # (suffix & 3)
            Symbol(mangledvar(op), modsuffix)
        elseif u₁unrolledsym
            Symbol(mvar, u)
        else
            mvar
        end
        for n ∈ 1:nparents
            if isloopvalue(parents_op[n])
                loopval = first(loopdependencies(parents_op[n]))
                add_loopvalue!(instrcall, loopval, ua, u)
            else
                parent = mangledvar(parents_op[n])
                if n == tiledouterreduction
                    parent = Symbol(parent, modsuffix)
                else
                    if parents_u₂syms[n]
                        parent = Symbol(parent, suffix_)
                    end
                    if parents_u₁syms[n]
                        parent = Symbol(parent, u)
                    end
                end
                if isvectorized && isload(parents_op[n]) && vectorized ∉ loopdependencies(parents_op[n]) && n != tiledouterreduction && !(parents_u₁syms[n] & parents_u₂syms[n])
                    broadcast_parent = Symbol(parent, "##broadcasted##")
                    if (parents_u₂syms[n] && iszero(u)) || (parents_u₁syms[n] && (isnothing(suffix) || iszero(suffix))) || (!(parents_u₁syms[n] | parents_u₂syms[n]) && iszero(u) && (isnothing(suffix) || iszero(suffix)))
                        push!(q.args, Expr(:(=), broadcast_parent, Expr(:call, lv(:vbroadcast), VECTORWIDTHSYMBOL, parent)))
                    end
                    parent = broadcast_parent
                end
                push!(instrcall.args, parent)
            end
        end
        if maskreduct && (u == Uiter || u₁loopsym !== vectorized) # only mask last
            if last(instrcall.args) == varsym
                pushfirst!(instrcall.args, lv(:vifelse))
                insert!(instrcall.args, 3, mask)
            else
                push!(q.args, Expr(:(=), varsym, Expr(:call, lv(:vifelse), mask, instrcall, varsym)))
                continue
            end
        end
        if instr.instr === :identity && isone(length(parents_op))
            push!(q.args, Expr(:(=), varsym, instrcall.args[2]))
        else
            push!(q.args, Expr(:(=), varsym, instrcall))
        end
    end
end


text:
# Autogenerated wrapper script for ZITSOL_1_jll for i686-linux-gnu-libgfortran3
export libZITSOL_1

using LAPACK_jll
using OpenBLAS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ZITSOL_1")
JLLWrappers.@declare_library_product(libZITSOL_1, "libZITSOL_1.so")
function __init__()
    JLLWrappers.@generate_init_header(LAPACK_jll, OpenBLAS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libZITSOL_1,
        "lib/libZITSOL_1.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module CoordinateDescent

using ProximalBase
using DataStructures: nlargest
using SparseArrays
using Statistics
using LinearAlgebra

export
  lasso, sqrtLasso, feasibleLasso!, scaledLasso, scaledLasso!,
  LassoPath, refitLassoPath,
  IterLassoOptions,

  # CD
  CoordinateDifferentiableFunction,
  CDLeastSquaresLoss, CDWeightedLSLoss, CDQuadraticLoss, CDSqrtLassoLoss,
  CDOptions,
  coordinateDescent, coordinateDescent!,

  # var coef
  GaussianKernel, SmoothingKernel, EpanechnikovKernel, evaluate, createKernel,
  locpoly, locpolyl1

include("utils.jl")
include("atom_iterator.jl")
include("cd_differentiable_function.jl")
include("coordinate_descent.jl")
include("lasso.jl")
include("varying_coefficient_lasso.jl")

end


text:
# using Revise
using Test, BifurcationKit, LinearAlgebra, SparseArrays, Setfield, Parameters, ForwardDiff
const BK = BifurcationKit

f1(u, v) = u^2 * v
norminf = x -> norm(x, Inf)

function Fbru!(f, x, p)
	@unpack α, β, D1, D2, l = p
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h
	c1 = D1 / l^2 / h2
	c2 = D2 / l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	# Dirichlet boundary conditions
	f[1]   = c1 * (α	  - 2u[1] + u[2] ) + α - (β + 1) * u[1] + f1(u[1], v[1])
	f[end] = c2 * (v[n-1] - 2v[n] + β / α)			 + β * u[n] - f1(u[n], v[n])

	f[n]   = c1 * (u[n-1] - 2u[n] +  α  )  + α - (β + 1) * u[n] + f1(u[n], v[n])
	f[n+1] = c2 * (β / α  - 2v[1] + v[2])			 + β * u[1] - f1(u[1], v[1])

	for i=2:n-1
		  f[i] = c1 * (u[i-1] - 2u[i] + u[i+1]) + α - (β + 1) * u[i] + f1(u[i], v[i])
		f[n+i] = c2 * (v[i-1] - 2v[i] + v[i+1])			  + β * u[i] - f1(u[i], v[i])
	end
	return f
end

Fbru(x, p) = Fbru!(similar(x), x, p)

function Jbru_sp(x, p)
	@unpack α, β, D1, D2, l = p
	# compute the Jacobian using a sparse representation
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h

	c1 = D1 / p.l^2 / h2
	c2 = D2 / p.l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	diag   = zeros(eltype(x), 2n)
	diagp1 = zeros(eltype(x), 2n-1)
	diagm1 = zeros(eltype(x), 2n-1)

	diagpn = zeros(eltype(x), n)
	diagmn = zeros(eltype(x), n)

	@. diagmn = β - 2 * u * v
	@. diagm1[1:n-1] = c1
	@. diagm1[n+1:end] = c2

	@. diag[1:n]    = -2c1 - (β + 1) + 2 * u * v
	@. diag[n+1:2n] = -2c2 - u * u

	@. diagp1[1:n-1] = c1
	@. diagp1[n+1:end] = c2

	@. diagpn = u * u
	J = spdiagm(0 => diag, 1 => diagp1, -1 => diagm1, n => diagpn, -n => diagmn)
	return J
end

Jbru_ana(x, p) = ForwardDiff.jacobian(z->Fbru(z,p),x)
jet = BK.get3Jet(Fbru, Jbru_ana)

n = 100
par_bru = (α = 2., β = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
	sol0 = vcat(par_bru.α * ones(n), par_bru.β/par_bru.α * ones(n))

# test that the jacobian is well computed
@test Jbru_sp(sol0, par_bru) - Jbru_ana(sol0, par_bru) |> sparse |> nnz == 0

opt_newton = NewtonPar(tol = 1e-11, verbose = false)
	out, hist, flag = newton(jet[1], jet[2], sol0 .* (1 .+ 0.01rand(2n)), par_bru, opt_newton)

opts_br0 = ContinuationPar(dsmin = 0.001, dsmax = 0.1, ds= 0.01, pMax = 1.8, newtonOptions = opt_newton, detectBifurcation = 3, nev = 16, nInversion = 4)
	br, = continuation(jet[1], jet[2], out, (@set par_bru.l = 0.3), (@lens _.l), opts_br0, printSolution = (x, p) -> norm(x, Inf64))
###################################################################################################
# Hopf continuation with automatic procedure
outhopf, = newton(jet[1], jet[2], br, 1; startWithEigen = true, d2F = jet[3])
outhopf, = newtonHopf(jet[1], jet[2], br, 1; startWithEigen = true)
optconthopf = ContinuationPar(dsmin = 0.001, dsmax = 0.15, ds= 0.01, pMax = 6.8, pMin = 0., newtonOptions = opt_newton, maxSteps = 50, detectBifurcation = 2)
outhopfco, = continuationHopf(jet[1], jet[2], br, 1, (@lens _.β), optconthopf; startWithEigen = true, updateMinAugEveryStep = 1, plot = false)

# Continuation of the Hopf Point using Dense method
ind_hopf = 1
# av = randn(Complex{Float64},2n); av = av./norm(av)
# bv = randn(Complex{Float64},2n); bv = bv./norm(bv)
hopfpt = BK.HopfPoint(br, ind_hopf)
bifpt = br.specialpoint[ind_hopf]
hopfvariable = HopfProblemMinimallyAugmented(
					jet[1], jet[2], nothing, nothing,
					(@lens _.l),
					conj.(br.eig[bifpt.idx].eigenvec[:, bifpt.ind_ev]),
					(br.eig[bifpt.idx].eigenvec[:, bifpt.ind_ev]),
					# av,
					# bv,
					opts_br0.newtonOptions.linsolver)

hopfvariable(hopfpt, par_bru) |> norm

Bd2Vec(x) = vcat(x.u, x.p)
Vec2Bd(x) = BorderedArray(x[1:end-2], x[end-1:end])
hopfpbVec(x, p) = Bd2Vec(hopfvariable(Vec2Bd(x),p))

# finite differences Jacobian
Jac_hopf_fdMA(u0, p) = ForwardDiff.jacobian( u-> hopfpbVec(u, p), u0)
# ``analytical'' jacobian
Jac_hopf_MA(u0, p, pb::HopfProblemMinimallyAugmented) = (return (x=u0,params=p ,hopfpb=pb))

rhs = rand(length(hopfpt))
jac_hopf_fd = Jac_hopf_fdMA(Bd2Vec(hopfpt), par_bru)
sol_fd = jac_hopf_fd \ rhs

# create a linear solver
hopfls = BK.HopfLinearSolverMinAug()
tmpVecforσ = zeros(ComplexF64, 2+2n)
sol_ma,  = hopfls(Jac_hopf_MA(hopfpt, par_bru, hopfvariable), BorderedArray(rhs[1:end-2],rhs[end-1:end]), debugArray = tmpVecforσ)

# we test the expression for σp
σp_fd = Complex(jac_hopf_fd[end-1,end-1], jac_hopf_fd[end,end-1])
σp_fd_ana = tmpVecforσ[1]
@test σp_fd ≈ σp_fd_ana rtol = 1e-5

# we test the expression for σω
σω_fd = Complex(jac_hopf_fd[end-1,end], jac_hopf_fd[end,end])
σω_fd_ana = tmpVecforσ[2]
@test σω_fd ≈ σω_fd_ana rtol = 1e-4

# we test the expression for σx
σx_fd = jac_hopf_fd[end-1, 1:end-2] + Complex(0,1) * jac_hopf_fd[end, 1:end-2]
σx_fd_ana = tmpVecforσ[3:end]
@test σx_fd ≈ σx_fd_ana rtol = 1e-3

outhopf, hist, flag = newton(
		jet[1], jet[2], br, 1, Jᵗ = (x, p) -> transpose(Jbru_sp(x, p)))
		# flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", outhopf.p[1], ", ω = ", outhopf.p[end], ", from l = ",hopfpt.p[1],"\n")

outhopf, _, flag, _ = newton((u, p) -> hopfvariable(u, p),
							(x, p) -> Jac_hopf_MA(x, p, hopfvariable),
							hopfpt, par_bru, NewtonPar(verbose = false, linsolver = BK.HopfLinearSolverMinAug()))
	# flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", outhopf.p[1], ", ω = ", outhopf.p[2], " from ", hopfpt.p, "\n")

# version with analytical Hessian = 2 P(du2) P(du1) QU + 2 PU P(du1) Q(du2) + 2PU P(du2) Q(du1)
function d2F(x, p1, du1, du2)
	n = div(length(x),2)
	out = similar(du1)
	@views out[1:n] .= 2 .* x[n+1:end] .* du1[1:n] .* du2[1:n] .+
				2 .* x[1:n] .* du1[1:n] .* du2[n+1:end] .+
				2 .* x[1:n] .* du2[1:n] .* du1[1:n]
	@inbounds for ii=1:n
		out[ii+n] = -out[ii]
	end
	return out
end

outhopf, hist, flag = newton(jet[1], jet[2], br, 1,
		Jᵗ = (x, p) -> transpose(Jbru_sp(x, p)),
		d2F = (x, p1, v1, v2) -> d2F(x, 0., v1, v2))
		# flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", outhopf.p[1], ", ω = ", outhopf.p[end], ", from l = ",hopfpt.p[1],"\n")

br_hopf, u1_hopf = continuation(
			jet[1], jet[2], br, ind_hopf, (@lens _.β),
			ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 6.5, pMin = 0.0, a = 2., theta = 0.4, maxSteps = 3, newtonOptions = NewtonPar(verbose = false)), verbosity = 0, plot = false)

br_hopf, u1_hopf = continuation(jet[1], jet[2], br, ind_hopf, (@lens _.β), ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 6.5, pMin = 0.0, a = 2., theta = 0.4, maxSteps = 3, newtonOptions = NewtonPar(verbose = false)), Jᵗ = (x, p) ->  transpose(Jbru_sp(x, p)), d2F = (x, p1, v1, v2) -> d2F(x, 0., v1, v2), verbosity = 0, plot = false)
#################################################################################################### Continuation of Periodic Orbit
ind_hopf = 1
hopfpt = BK.HopfPoint(br, ind_hopf)

l_hopf = hopfpt.p[1]
ωH	   = hopfpt.p[2] |> abs
M = 20

orbitguess = zeros(2n, M)
phase = []; scalphase = []
vec_hopf = geteigenvector(opt_newton.eigsolver, br.eig[br.specialpoint[ind_hopf].idx][2], br.specialpoint[ind_hopf].ind_ev-1)
for ii=1:M
	t = (ii-1)/(M-1)
	orbitguess[:, ii] .= real.(hopfpt.u + 26*0.1 * vec_hopf * exp(-2pi * complex(0, 1) * (t - .252)))
end

orbitguess_f = vcat(vec(orbitguess), 2pi/ωH) |> vec

# test guess using function
l_hopf, Th, orbitguess2, hopfpt, vec_hopf = BK.guessFromHopf(br, ind_hopf, opt_newton.eigsolver, M, 2.6; phase = 0.252)

poTrap = PeriodicOrbitTrapProblem(jet[1], Jbru_sp, real.(vec_hopf), hopfpt.u, M, 2n	)

jac_PO_fd = BK.finiteDifferences(x -> poTrap(x, (@set par_bru.l = l_hopf + 0.01)), orbitguess_f)
jac_PO_sp = poTrap(Val(:JacFullSparse), orbitguess_f, (@set par_bru.l = l_hopf + 0.01))

# test of the Jacobian for PeriodicOrbit via Finite differences VS the FD associated jacobian
# test jacobian expression for Periodic Orbit solve problem
@test norm(jac_PO_fd - jac_PO_sp, Inf64) < 1e-4

# test various jacobians and methods
jac_PO_sp =  poTrap(Val(:BlockDiagSparse), orbitguess_f, (@set par_bru.l = l_hopf + 0.01))
BK.getTimeDiff(poTrap, orbitguess_f)
# BK.Jc(poTrap, reshape(orbitguess_f[1:end-1], 2n, M), par_bru, reshape(orbitguess_f[1:end-1], 2n, M))
# BK.Jc(poTrap, orbitguess_f, par_bru, orbitguess_f)

# newton to find Periodic orbit
opt_po = NewtonPar(tol = 1e-8, verbose = false, maxIter = 150)
	outpo_f, _, flag = @time newton(
		(x, p) ->  poTrap(x, p),
		(x, p) ->  poTrap(Val(:JacFullSparse),x,p),
		copy(orbitguess_f), (@set par_bru.l = l_hopf + 0.01), opt_po)
	# println("--> T = ", outpo_f[end])
# flag && printstyled(color=:red, "--> T = ", outpo_f[end], ", amplitude = ", BK.amplitude(outpo_f, n, M; ratio = 2),"\n")

newton(poTrap, orbitguess_f, (@set par_bru.l = l_hopf + 0.01), opt_po; linearPO = :FullLU)

# jacobian of the functional
Jpo2 = poTrap(Val(:JacCyclicSparse), orbitguess_f, (@set par_bru.l = l_hopf + 0.01))

# calcul des exposants de Floquet, extract full vector
BK.MonodromyQaD(Val(:ExtractEigenVector), poTrap, orbitguess_f, par_bru, orbitguess_f[1:2n])

# calcul des exposants de Floquet
floquetES = FloquetQaD(DefaultEig())

# continuation of periodic orbits using :BorderedLU linear algorithm
opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.05, ds= 0.001, pMax = 2.3, maxSteps = 3, theta = 0.1, newtonOptions = NewtonPar(verbose = false), detectBifurcation = 1)
	br_pok2, = continuation(
		poTrap, orbitguess_f, (@set par_bru.l = l_hopf + 0.01), (@lens _.l), opts_po_cont; linearPO = :BorderedLU,
		plot = false, verbosity = 0)

# test of simple calls to newton / continuation
deflationOp = DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [zero(orbitguess_f)])
# opt_po = NewtonPar(tol = 1e-8, verbose = false, maxIter = 15)
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.03, ds= 0.01, pMax = 3.0, maxSteps = 3, newtonOptions = (@set opt_po.verbose = false), nev = 2, precisionStability = 1e-8, detectBifurcation = 1)
for linalgo in [:FullLU, :BorderedLU, :FullSparseInplace]
	@show linalgo
	# with deflation
	@time newton(poTrap,
			copy(orbitguess_f), (@set par_bru.l = l_hopf + 0.01), opt_po, deflationOp; linalgo = linalgo, normN = norminf)
	# classic Newton-Krylov
	outpo_f, hist, flag = @time newton(poTrap,
			copy(orbitguess_f), (@set par_bru.l = l_hopf + 0.01), opt_po; linalgo = linalgo, normN = norminf)
	# continuation
	br_pok2, = @time continuation(poTrap,
			copy(orbitguess_f), (@set par_bru.l = l_hopf + 0.01), (@lens _.l),
			opts_po_cont; linearPO = linalgo, verbosity = 0,
			plot = false, normC = norminf)
end

# test of Matrix free computation of Floquet exponents
eil = EigKrylovKit(x₀ = rand(2n))
ls = GMRESKrylovKit()
ls = DefaultLS()
opt_po = NewtonPar(tol = 1e-8, verbose = false, maxIter = 10, linsolver = ls, eigsolver = eil)
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.03, ds= 0.01, pMax = 3.0, maxSteps = 3, newtonOptions = (@set opt_po.verbose = false), nev = 2, precisionStability = 1e-8, detectBifurcation = 2)
br_pok2, upo, = continuation(poTrap, outpo_f, (@set par_bru.l = l_hopf + 0.01), (@lens _.l), opts_po_cont; linearPO = :FullLU, normC = norminf, verbosity = 0)


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

# The order in which files are included is important,
# since certain things need to be defined before others use them

## Interval type

abstract AbstractInterval <: Real

immutable Interval{T<:Real} <: AbstractInterval
    lo :: T
    hi :: T

    function Interval(a::Real, b::Real)

        if a > b
            (isinf(a) && isinf(b)) && return new(a, b)  # empty interval = [∞,-∞]

            throw(ArgumentError("Must have a ≤ b to construct Interval(a, b)."))
        end

        new(a, b)
    end
end

## Outer constructors

Interval{T<:Real}(a::T, b::T) = Interval{T}(a, b)
Interval{T<:Real}(a::T) = Interval(a, a)
Interval(a::Tuple) = Interval(a...)
Interval{T<:Real, S<:Real}(a::T, b::S) = Interval(promote(a,b)...)

## Concrete constructors for Interval, to effectively deal only with Float64,
# BigFloat or Rational{Integer} intervals.
Interval{T<:Integer}(a::T, b::T) = Interval(float(a), float(b))
Interval{T<:Irrational}(a::T, b::T) = Interval(float(a), float(b))

eltype{T<:Real}(x::Interval{T}) = T

## Include files
include("special.jl")
include("macros.jl")
include("conversion.jl")
include("precision.jl")
include("set_operations.jl")
include("arithmetic.jl")
include("functions.jl")
include("trigonometric.jl")
include("hyperbolic.jl")

# Syntax for intervals

a..b = @interval(a, b)

macro I_str(ex)  # I"[3,4]"
    @interval(ex)
end

a ± b = (a-b)..(a+b)