text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

function isapprox(a1::Assembly, a2::Assembly)
    T = isapprox(a1.K, a2.K)
    T &= isapprox(a1.C1, a2.C1)
    T &= isapprox(a1.C2, a2.C2)
    T &= isapprox(a1.D, a2.D)
    T &= isapprox(a1.f, a2.f)
    T &= isapprox(a1.g, a2.g)
    return T
end

function assemble_prehook!(::Problem, ::T) where T<:Number end

function assemble_posthook!(::Problem, ::T) where T<:Number end

# will be deprecated
function assemble!(::Assembly, ::Problem{P}, ::Element, ::Any) where P
    @warn("One must define assemble! function for problem of type $P. " *
          "Not doing anything.")
    return nothing
end

# will be deprecated
function assemble!(assembly::Assembly, problem::Problem{P},
                   elements::Vector{Element}, time) where P
    @warn("This is default assemble! function. Decreased performance can be " *
          "expected without preallocation of memory. One should implement " *
          "`assemble_elements!(problem, assembly, elements, time)` function.")
    for element in elements
        assemble!(assembly, problem, element, time)
    end
    return nothing
end

"""
    assemble_elements!(problem, assembly, elements, time)

Assemble elements for problem.

This should be overridden with own `assemble_elements!`-implementation.
"""
function assemble_elements!(problem::Problem, assembly::Assembly,
                      elements::Vector{Element{E}}, time) where E
    elements2 = convert(Vector{Element}, elements)
    assemble!(assembly, problem, elements2, time)
end

function assemble!(problem::Problem, time)

    assemble_prehook!(problem, time)
    elements = get_elements(problem)
    assembly = get_assembly(problem)

    if !isempty(assembly)
        @warn("Problem assembly is not empty before assembling. This is probably " *
              "causing unexpected results. To remove old assembly, use " *
              "`empty!(problem.assembly)`", typeof(problem), problem.name)
        assemble_posthook!(problem, time)
        return nothing
    end

    if isempty(elements)
        @warn("There is no elements defined in problem. Before assembling a " *
              "problem, elements must be added using " *
              "`add_elements!(problem, elements)`.", typeof(problem), problem.name)
        assemble_posthook!(problem, time)
        return nothing
    end

    first_element = first(elements)
    unknown_field_name = get_unknown_field_name(problem)
    if !haskey(first_element, unknown_field_name)
        #=
        warn("Assembling elements for problem $(problem.name): seems that ",
             "problem is uninitialized. To initialize problem, use ",
             "`initialize!(problem, time)`.")
        info("Initializing problem $(problem.name) at time $time automatically.")
        =#
        initialize!(problem, time)
    end

    for (element_type, elements) in group_by_element_type(elements)
        assemble_elements!(problem, assembly, elements, time)
    end
    assemble_posthook!(problem, time)
    return nothing
end

function assemble!(problem::Problem)
    @warn("assemble!(problem) will be deprecated. Use assemble!(problem, time)")
    assemble!(problem, 0.0)
end

function assemble_mass_matrix!(problem::Problem, time::Float64)
    if !isempty(problem.assembly.M)
        @info("Mass matrix for is already assembled, not assembling.",
              problem.name)
        return nothing
    end
    elements = get_elements(problem)
    for (element_type, elements) in group_by_element_type(get_elements(problem))
        assemble_mass_matrix!(problem::Problem, elements, time)
    end
    return nothing
end

function assemble_mass_matrix!(problem::Problem, elements::Vector{Element{Basis}}, time) where Basis
    nnodes = length(first(elements))
    dim = get_unknown_field_dimension(problem)
    M = zeros(nnodes, nnodes)
    N = zeros(1, nnodes)
    NtN = zeros(nnodes, nnodes)
    ldofs = zeros(Int, nnodes)
    for element in elements
        fill!(M, 0.0)
        for ip in get_integration_points(element, 2)
            detJ = element(ip, time, Val{:detJ})
            rho = element("density", ip, time)
            w = ip.weight*rho*detJ
            eval_basis!(Basis, N, ip)
            N = element(ip, time)
            mul!(NtN, transpose(N), N)
            rmul!(NtN, w)
            for i=1:nnodes^2
                M[i] += NtN[i]
            end
        end
        for (i, j) in enumerate(get_connectivity(element))
            @inbounds ldofs[i] = (j-1)*dim
        end
        for i=1:dim
            add!(problem.assembly.M, ldofs.+i, ldofs.+i, M)
        end
    end
    return
end

"""
    assemble_mass_matrix!(problem, elements::Vector{Element{Tet10}}, time)

Assemble Tet10 mass matrices using special method. If Tet10 has constant metric
if can be integrated analytically to gain performance.
"""
function assemble_mass_matrix!(problem::Problem, elements::Vector{Element{FEMBasis.Tet10}}, time)
    nnodes = length(Tet10)
    dim = get_unknown_field_dimension(problem)
    M = zeros(nnodes, nnodes)
    N = zeros(1, nnodes)
    NtN = zeros(nnodes, nnodes)
    ldofs = zeros(Int, nnodes)

    M_CM = 1.0/2520.0 * [
         6  1  1  1 -4 -6 -4 -4 -6 -6
         1  6  1  1 -4 -4 -6 -6 -4 -6
         1  1  6  1 -6 -4 -4 -6 -6 -4
         1  1  1  6 -6 -6 -6 -4 -4 -4
        -4 -4 -6 -6 32 16 16 16 16  8
        -6 -4 -4 -6 16 32 16  8 16 16
        -4 -6 -4 -6 16 16 32 16  8 16
        -4 -6 -6 -4 16  8 16 32 16 16
        -6 -4 -6 -4 16 16  8 16 32 16
        -6 -6 -4 -4  8 16 16 16 16 32]

    function is_CM(::Element{FEMBasis.Tet10}, X; rtol=1.0e-6)
        isapprox(X[5],  1/2*(X[1]+X[2]); rtol=rtol) || return false
        isapprox(X[6],  1/2*(X[2]+X[3]); rtol=rtol) || return false
        isapprox(X[7],  1/2*(X[3]+X[1]); rtol=rtol) || return false
        isapprox(X[8],  1/2*(X[1]+X[4]); rtol=rtol) || return false
        isapprox(X[9],  1/2*(X[2]+X[4]); rtol=rtol) || return false
        isapprox(X[10], 1/2*(X[3]+X[4]); rtol=rtol) || return false
        return true
    end

    n_CM = 0
    for element in elements
        for (i, j) in enumerate(get_connectivity(element))
            @inbounds ldofs[i] = (j-1)*dim
        end

        X = element("geometry", time)
        rho = element("density", time)
        if is_CM(element, X) && length(rho) == 1
            ip = (1.0/3.0, 1.0/3.0, 1.0/3.0)
            detJ = element(ip, time, Val{:detJ})
            rho = element("density", ip, time)
            CM_s = detJ*rho
            n_CM += 1
            for i=1:dim
                add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, CM_s * M_CM)
            end
        else
            fill!(M, 0.0)
            for ip in get_integration_points(element, 2)
                detJ = element(ip, time, Val{:detJ})
                rho = element("density", ip, time)
                w = ip.weight*rho*detJ
                eval_basis!(FEMBasis.Tet10, N, ip)
                N = element(ip, time)
                mul!(NtN, transpose(N), N)
                rmul!(NtN, w)
                for i=1:nnodes^2
                    M[i] += NtN[i]
                end
            end
            for i=1:dim
                add!(problem.assembly.M, ldofs .+ i, ldofs .+ i, M)
            end
        end
    end
    @info("$n_CM of $(length(elements)) was constant metric.")
    return
end


text:
export fit

# Default solvers

# TODO: in the future, have cases where if the things are too big, take another
# default. Also should check if p > n in which case should do dual stuff (or
# other appropriate alternative)

# Linear, Ridge
_solver(::GLR{L2Loss,<:L2R}, np::NTuple{2,Int}) = Analytical()

# Logistic, Multinomial
_solver(::GLR{LogisticLoss,<:L2R}, 	  np::NTuple{2,Int}) = LBFGS()
_solver(::GLR{MultinomialLoss,<:L2R}, np::NTuple{2,Int}) = LBFGS()

# Lasso, ElasticNet, Logistic, Multinomial
function _solver(glr::GLR{<:SmoothLoss,<:ENR}, np::NTuple{2,Int})
	(is_l1(glr.penalty) || is_elnet(glr.penalty)) && return FISTA()
	@error "Not yet implemented."
end

# Robust, Quantile
_solver(::GLR{<:RobustLoss,<:L2R}, np::NTuple{2,Int}) = LBFGS()

# Fallback NOTE: should revisit bc with non-smooth, wouldn't work probably
# PGD/PSGD depending on how much data there is
_solver(::GLR, np::NTuple{2,Int}) = @error "Not yet implemented."

"""
$SIGNATURES

Fit a generalised linear regression model using an appropriate solver based on
the loss and penalty of the model. A method can, in some cases, be specified.
"""
function fit(glr::GLR, X::AbstractMatrix{<:Real}, y::AVR;
			 solver::Solver=_solver(glr, size(X)))
    check_nrows(X, y)
	n, p = size(X)
	p += Int(glr.fit_intercept)
	# allocate cache for temporary computations of size n/p
	# which are frequent but otherwise un-important so that
	# we can reduce the overall number of allocations
	# these are const Refs defined when the module is loaded
	c = glr.loss isa MultinomialLoss ? maximum(y) : 0
	allocate(n, p, c)
	# effective call to fit routine
    θ = _fit(glr, solver, X, y)
	# de-allocate cache
	deallocate()
	return θ
end

ext: jl


text:
# additional functions for use with LinkedLists.jl

function Base.getindex(lst::LinkedLists.SLinkedList, idx::Int)
    !(0 < idx <= length(lst)) && throw(BoundsError("index $idx"))
    return lst[positiontoindex(idx, lst)]
end

function Base.getindex(lst::LinkedLists.LinkedList, idx::Int)
    !(0 < idx <= length(lst)) && throw(BoundsError("index $idx"))
    return lst[positiontoindex(idx, lst)]
end


text:
# ----------------------------------------------------

"anything with coords in a scene"
abstract SceneItem

# "layers are containers in a scene... nothing to draw directly"
# abstract Layer <: SceneItem

"these are things that need to be drawn"
abstract Shape <: SceneItem

# ----------------------------------------------------

type Scene
  canvas
  items::Vector{SceneItem}
end

Scene(canvas::Cxx.CppPtr) = Scene(canvas, SceneItem[])

const _scenes = Dict{Int, Scene}()

# NOTE: this is the callback that gets called from within MyCanvas::paintEvent
#       ALL drawing should be done before returning
function draw(canvas::Cxx.CppPtr)
  # dump(canvas)

  # get the index of the canvas
  idx = @cxx canvas->getidx()
  # @show idx

  # get the starting coordinates
  # x = @cxx canvas->x()
  # y = @cxx canvas->y()
  width = @cxx canvas->width()
  height = @cxx canvas->height()
  # @show x y width height
  # box = SceneBox(Nullable{SceneBox}(), P2(x,y), P2(width,height))
  box = BoundingBox(0, 0, width, height)
  # @show box

  # grab the scene object
  scene = _scenes[idx]
  # dump(scene)

  # set up the painter
  painter = Painter(@cxxnew QPainter(canvas))
  @cxx (painter.o)->setRenderHint(@cxx(QPainter::Antialiasing), true)

  # draw something
  for item in scene.items
    draw(item, painter, box)
  end

  nothing
end

function Scene()
  # get the next available index
  idx = isempty(_scenes) ? 1 : maximum(keys(_scenes)) + 1

  # create a new canvas, and set the window title
  canvas = @cxxnew MyCanvas(idx)
  @cxx canvas->setWindowTitle(pointer("idx: $idx"))

  # create, store, and return the Scene object
  scene = Scene(canvas)
  _scenes[idx] = scene
end

function Base.display(scene::Scene)
  @cxx scene.canvas->update()
  @cxx scene.canvas->show()
end

Base.push!(scene::Scene, item::SceneItem) = push!(scene.items, item)

# ----------------------------------------------------

# ----------------------------------------------------

# @enum MetricType PercentMetric PixelMetric

# immutable Metric{M}
#   val::Float64
# end

# typealias Px Metric{PixelMetric}
# typealias Pct Metric{PercentMetric}

# Base.convert(::Type{Px}, x::Real) = Px(float(x))
# Base.convert(::Type{Pct}, x::Real) = Pct(float(x))
# Base.promote_rule{T<:Real}(::Type{Px}, ::Type{T}) = Px
# Base.promote_rule{T<:Real}(::Type{Pct}, ::Type{T}) = Pct

# Base.show(io::IO, p::Px) = print(io, "Pixel{", p.val, "}")
# Base.show(io::IO, p::Pct) = print(io, "Percent{", p.val, "}")

# for op in (:+, :-, :*, :/)
#   @eval $op(p1::Px, p2::Px) = Px($op(p1.val, p2.val))
#   @eval $op(p1::Pct, p2::Pct) = Px($op(p1.val, p2.val))
# end

# "This allows combination of percent and pixel metrics"
# immutable Distance
#   pct::Float64
#   px::Float64
# end

# Base.convert(::Type{Distance}, pct::Pct) = Distance(pct.val, 0.0)
# Base.convert(::Type{Distance}, px::Px) = Distance(0.0, px.val)
# Base.convert(::Type{Distance}, x::Real) = Distance(float(x), 0.0)
# Base.promote_rule{T<:Union{Pct,Px,Real}}(::Type{Distance}, ::Type{T}) = Distance

# for op in (:+, :-, :*, :/)
#   @eval $op(d1::Distance, d2::Distance) = Distance($op(d1.pct, d2.pct), $op(d1.px, d2.px))
# end

# typealias DistOrDists Union{Distance, Vector{Distance}}

# ---------------------------------------------------

# "defines a view inside the widget area... all values in pixels"
# immutable ViewBox
#   x::Float64
#   y::Float64
#   w::Float64
#   h::Float64
# end

# px_x(dist::Distance, box::ViewBox) = box.x + (dist.pct * box.w + dist.px)
# px_y(dist::Distance, box::ViewBox) = box.y + (dist.pct * box.h + dist.px)
# px_w(dist::Distance, box::ViewBox) = dist.pct * box.w + dist.px
# px_h(dist::Distance, box::ViewBox) = dist.pct * box.h + dist.px

# ---------------------------------------------------

"holds a list of child scene items and the coordinates"
type View <: SceneItem
  box::BoundingBox
  items::Vector{SceneItem}
end

function draw(view::View, painter, box::BoundingBox)
  for item in view.items
    draw(item, painter, view.box)
  end
end

Base.push!(v::View, item::Shape) = push!(v.items, item)

# ---------------------------------------------------

immutable Pen{T<:Cxx.CppValue}
  o::T
end

function Pen(cstr::AbstractString)
  Pen(@cxx QPen(pointer(cstr)))
end

immutable Brush{T<:Cxx.CppValue}
  o::T
end

function Brush(cstr::AbstractString)
  Brush(@cxx QBrush(pointer(cstr)))
end

# ---------------------------------------------------

type Painter{T<:Cxx.CppPtr}
  o::T
  pen::Nullable{Pen}
  brush::Nullable{Brush}
end

Painter(o) = Painter(o, Nullable{Pen}(), Nullable{Brush}())

function pen!(painter::Painter, pen::Pen)
  !isnull(painter.pen) && get(painter.pen) === pen && return
  painter.pen = Nullable(pen)
  @cxx (painter.o)->setPen(pen.o)
end
function brush!(painter::Painter, brush::Brush)
  !isnull(painter.brush) && get(painter.brush) === brush && return
  painter.brush = Nullable(brush)
  @cxx (painter.o)->setBrush(brush.o)
end

# ---------------------------------------------------

type Ellipse <: Shape
  x::Distance
  y::Distance
  w::Distance
  h::Distance
  pen::Pen
  brush::Brush
end

function draw(item::Ellipse, painter, box::BoundingBox)
  # xy, wh = convertPctToSceneCoords(box, center, radius)
  x = px_x(item.x, box)
  y = px_y(item.y, box)
  w = px_w(item.w, box)
  h = px_h(item.h, box)
  pen!(painter, item.pen)
  brush!(painter, item.brush)
  @cxx (painter.o)->drawEllipse(x-0.5*w, y-0.5*h, w, h)
end

type Ellipses <: Shape
  n::Int
  x::Vector{Distance}
  y::Vector{Distance}
  w::Vector{Distance}
  h::Vector{Distance}
  pen::Vector{Pen}
  brush::Vector{Brush}
end

function draw(item::Ellipses, painter, box::BoundingBox)
  nx = length(item.x)
  ny = length(item.y)
  nw = length(item.w)
  nh = length(item.h)
  np = length(item.pen)
  nb = length(item.brush)
  x = 0.0
  y = 0.0
  w = 0.0
  h = 0.0
  x = nx == 1 ? px_x(item.x[1], box) : 0.0
  y = ny == 1 ? px_y(item.y[1], box) : 0.0
  w = nw == 1 ? px_w(item.w[1], box) : 0.0
  h = nh == 1 ? px_h(item.h[1], box) : 0.0
  np == 1 && @cxx painter->setPen(item.pen[1].o)
  nb == 1 && @cxx painter->setBrush(item.brush[1].o)

  for i in 1:n
    if nx > 1
      x = px_x(item.x[mod1(i,nx)], box)
    end
    if ny > 1
      y = px_y(item.y[mod1(i,ny)], box)
    end
    if nw > 1
      w = px_w(item.w[mod1(i,nw)], box)
    end
    if nh > 1
      h = px_h(item.h[mod1(i,nh)], box)
    end
    if np > 1
      @cxx painter->setPen(item.pen[mod1(i,np)].o)
    end
    if nb > 1
      @cxx painter->setBrush(item.brush[mod1(i,nb)].o)
    end
    @cxx painter->drawEllipse(x-0.5*w, y-0.5*h, w, h)
  end
end

# ---------------------------------------------------

# note: null parent means it's top level

# "This is the box which contains shapes... all child shapes are relative to this box"
# immutable SceneBox
#   parent::Nullable{SceneBox}
#   xy::P2
#   wh::P2
# end

# TODO: a better composition... maybe pass scene w/h through directly?
# TODO: can we structure so that we don't need a parent reference, and the top level converts pct to scene?
#     note: we can probably pull this off if we stay as percentages... if item is 50% of box, and box is 30% of it's
#           parent, then we can just return 15% of the box pct

# "given values where x/y/w/h are on the scale 0->1, convert to box coords"
# function convertPctToBoxScale(box::SceneBox, xy::P2, wh::P2)
#   xy .* box.wh, wh .* box.wh
# end

# "given values where x/y/w/h are on the scale 0->1, convert to scene coords"
# function convertPctToSceneCoords(box::SceneBox, xy::P2, wh::P2)
#   xy, wh = convertPctToBoxScale(box, xy, wh)
#   if isnull(box.parent)
#     return xy, wh
#   else
#     return convertPctToBoxCoords(get(box.parent), xy, wh)
#   end
# end

# ----------------------------------------------------

# "holds a list of child scene items and the coordinates"
# type ViewLayer <: Layer
#   box::SceneBox
#   items::Vector{SceneItem}
# end

# function draw(layer::ViewLayer, painter, box::SceneBox)
#   for item in items
#     draw(item, painter, layer.box)
#   end
# end

# ----------------------------------------------------

# type Ellipse <: Shape
#   center::P2
#   radius::P2
#   pen
#   brush
# end

# function draw(item::Ellipse, painter, box::SceneBox)
#   xy, wh = convertPctToSceneCoords(box, center, radius)
#   @cxx painter->setPen(item.pen)
#   @cxx painter->setBrush(item.brush)
#   @cxx painter->drawEllipse(xy.x, xy.y, wh.w, wh.h)
# end

# push!(scene::Union{Scene,Layer}, item::Ellipse) = push!(scene.items, item)

# ----------------------------------------------------

# type Canvas
#   widget
# end

# function Canvas(idx::Int, title = "HELLO")
#   w = @cxxnew MyCanvas(idx)
#   @cxx w->setWindowTitle(pointer(title))
#   Canvas(w)
# end

# function Pen(cstr::ASCIIString)
#   @cxx QPen(pointer(cstr))
# end

# function Brush(cstr::ASCIIString)
#   @cxx QBrush(pointer(cstr))
# end

# Base.display(c::Canvas) = @cxx (c.widget)->show()

# ----------------------------------------------------

# function render(painter)
#   dump(painter)
#   nothing
# end

# ----------------------------------------------------

type Window
  o
end

function display(win::Window)
  @cxx (win.o)->show()
end

function window(idx::Int = 0; title = "HI", w = 1000, h = 500)
  Window(@cxxnew MyWindow(idx, pointer(title), w, h))
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds Base
# The underlying abstraction layer for DynamicBounds.
# See https://github.com/PSORLab/DynamicBoundsBase.jl
#############################################################################
# src/integrator.jl
# Defines attributes and functions used by integrators.
#############################################################################

"""
$(TYPEDEF)

Abstract supertype indicating the type of value returned at a specific index.
"""
abstract type AbstractBoundLoc end

"""
$(TYPEDEF)

Indicates the lower bound, relaxation, or (sub)gradient should be returned.
"""
struct Lower <: AbstractBoundLoc end

"""
$(TYPEDEF)

Indicates the upper bound, relaxation, or (sub)gradient should be returned.
"""
struct Upper <: AbstractBoundLoc end

"""
$(TYPEDEF)

Indicates the nominal value or (sub)gradient should be returned.
"""
struct Nominal <: AbstractBoundLoc end

"""
$(TYPEDEF)

The variety of the attribute is unspecified.
"""
struct Undefined <: AbstractBoundLoc end

"""
$(TYPEDEF)

Abstract supertype for any integrator that constructs relaxations of
a differential equation problem.
"""
abstract type AbstractDERelaxIntegrator end

"""
$(TYPEDEF)

Abstract supertype for any integrator that constructs relaxations of
an ordinary differential equation problem.
"""
abstract type AbstractODERelaxIntegrator <: AbstractDERelaxIntegrator end

"""
$(TYPEDEF)

Abstract supertype for attributes that can be used to `set` or
`get` attributes (properties) of variables in the model.
"""
abstract type AbstractIntegratorAttribute end

#This allows to use attributes in broadcast calls without the need to
# embed it in a `Ref`
Base.broadcastable(attribute::AbstractIntegratorAttribute) = Ref(attribute)

"""
$(TYPEDEF)

An abstract supertype for index structures needed to reference parameters or time.
"""
abstract type AbstractDynamicIndex end

"""
$(TYPEDEF)

A time index `t` used to access the value at `t.i`-th time point.
"""
struct TimeIndex
    t::Int
end

"""
$(TYPEDEF)

An integrator attribute for the string identifying the integration scheme.
"""
struct IntegratorName <: AbstractIntegratorAttribute end

"""
$(TYPEDEF)

An integrator attribute retrieving the number of parameter variables (`p`) aka `np`.
"""
struct ParameterNumber end

"""
$(TYPEDEF)

An integrator attribute retrieving the number of state variables (`x`) aka `nx`.
"""
struct StateNumber end

"""
$(TYPEDEF)

Set to `true` if integrator should compute derivatives wrt p for the local
ode integration.
"""
struct LocalSensitivityOn <: AbstractIntegratorAttribute end

"""
$(TYPEDEF)

An integrator attribute retrieving the number of points in the support set of (`x`) aka `nt`.
"""
struct SupportNumber end

"""
$(TYPEDEF)

An integrator attribute for the Gradient.
"""
struct Gradient{T<:AbstractBoundLoc} <: AbstractIntegratorAttribute
    index::Int
    time::Float64
end
Gradient{Lower}() = Gradient{Lower}(-1, -Inf)
Gradient{Upper}() = Gradient{Upper}(-1, -Inf)
Gradient{Nominal}() = Gradient{Nominal}(-1, -Inf)
Gradient() = Gradient{Undefined}(-1, -Inf)

"""
Gradient{T}(i::TimeIndex) where {T <: AbstractBoundLoc}

This constructor builds a structure that references the gradient of
the relaxation at time index `i`.
"""
function Gradient{T}(i::TimeIndex) where {T <: AbstractBoundLoc}
    Gradient{T}(i.t, -Inf)
end

"""
Gradient{T}(x::Float64) where {T <: AbstractBoundLoc}

This constructor builds a structure that references the subgradient of
the relaxation at time = `x`.
"""
function Gradient{T}(x::Float64) where {T <: AbstractBoundLoc}
    Gradient{T}(-1, x)
end

"""
$(TYPEDEF)

An integrator attribute for the Subgradient.
"""
struct Subgradient{T<:AbstractBoundLoc} <: AbstractIntegratorAttribute
    index::Int
    time::Float64
end
Subgradient{Lower}() = Subgradient{Lower}(-1, -Inf)
Subgradient{Upper}() = Subgradient{Upper}(-1, -Inf)
Subgradient() = Subgradient{Undefined}(-1, -Inf)

"""
Subgradient{T}(i::TimeIndex)

This constructor builds a structure that references the subgradient of
the relaxation at time index `i`.
"""
Subgradient{T}(i::TimeIndex) where {T <: AbstractBoundLoc} = Subgradient{T}(i.t, -Inf)

"""
Subgradient{T}(x::Float64)

This constructor builds a structure that references the subgradient of
the relaxation at time = `x`.
"""
Subgradient{T}(x::Float64) where {T <: AbstractBoundLoc} = Subgradient{T}(-1, x)

"""
$(TYPEDEF)

An integrator attribute for value of local solution bounds.
"""
struct Value <: AbstractIntegratorAttribute
    index::Int
    time::Float64
end
Value() = Value(-1, -Inf)
Value(i::TimeIndex) = Value(i, -Inf)
Value(x::Float64) = Value(-1, x)

"""
$(TYPEDEF)

An integrator attribute for state bounds.
"""
struct Bound{T<:AbstractBoundLoc} <: AbstractIntegratorAttribute
    index::Int
    time::Float64
end
Bound{Lower}() = Bound{Lower}(-1, -Inf)
Bound{Upper}() = Bound{Upper}(-1, -Inf)
Bound() = Bound{Undefined}(-1, -Inf)

"""
Bound{T}(i::TimeIndex)

This constructor builds a structure that references the state
bounds at time index `i`.
"""
Bound{T}(i::TimeIndex) where {T <: AbstractBoundLoc} = Bound{T}(i.t, -Inf)

"""
Bound{T}(x::Float64)

This constructor builds a structure that references the state
bounds at time = `x`.
"""
Bound{T}(x::Float64) where {T <: AbstractBoundLoc} = Bound{T}(-1, x)

"""
$(TYPEDEF)

An integrator attribute for relaxations.
"""
struct Relaxation{T<:AbstractBoundLoc} <: AbstractIntegratorAttribute
    index::Int
    time::Float64
end
Relaxation{Lower}() = Relaxation{Lower}(-1, -Inf)
Relaxation{Upper}() = Relaxation{Upper}(-1, -Inf)
Relaxation() = Relaxation{Undefined}(-1, -Inf)

"""
Relaxation{T}(i::TimeIndex)

This constructor builds a structure that references the relaxation
at time index `i`.
"""
Relaxation{T}(i::TimeIndex) where {T <: AbstractBoundLoc} = Relaxation{T}(i.t, -Inf)

"""
Relaxation{T}(x::Float64)

This constructor builds a structure that references relaxation at time = `x`.
"""
Relaxation{T}(x::Float64) where {T <: AbstractBoundLoc} = Relaxation{T}(-1, x)

"""
$(TYPEDEF)

An integrator attribute for indicating bounds/relaxations are of numeric
solution.
"""
struct IsNumeric <: AbstractIntegratorAttribute end

"""
$(TYPEDEF)

An integrator attribute for indicating bounds & relaxations are of exact
solution.
"""
struct IsSolutionSet <: AbstractIntegratorAttribute end

@enum(TerminationStatusCode,
      COMPLETED,
      EMPTY,
      NAN,
      RELAXATION_NOT_CALLED,
      NUMERICAL_ERROR,
      LIMIT_EXCEEDED,
      INVALID_OPTION,
      OTHER_ERROR
)

"""
$(TYPEDEF)

A integrator attribute used to query the `TerminationStatusCode` of the integrator
on completion. Current termination status codes are:
- `COMPLETED`: The algorithm terminate successfully with bounds and relaxations.
- `EMPTY`: The algorithm terminated successfully but the solution set was empty for
         some points in `tspan`.
- `NAN`: The algorithm terminated but some values are not a number (usually indicating
         a domain violation was encoutered when computing relaxations).
- `RELAXATION_NOT_CALLED`: The relaxation has not yet been computed.
- `NUMERICAL_ERROR`: A numerical error was encountered.
- `LIMIT_EXCEEDED`: A preset limit was exceeded (number of steps and so on).
- `INVALID_OPTION`: An invalid option was set.
- `OTHER_ERROR`: Another error was encountered.
"""
struct TerminationStatus <: AbstractIntegratorAttribute end

"""
$(TYPEDEF)

A integrator attribute used to access the current parameter value.
"""
struct ParameterValue <: AbstractIntegratorAttribute
    i::Int
end
ParameterValue() = ParameterValue(-1)

"""
$(TYPEDEF)

A integrator attribute used to access the constant parameter values.
"""
struct ConstantParameterValue <: AbstractIntegratorAttribute
    i::Int
end
ConstantParameterValue() = ConstantParameterValue(-1)

"""
$(TYPEDEF)

A integrator attribute used to access the current parameter value.
"""
struct ParameterBound{T <: AbstractBoundLoc} <: AbstractIntegratorAttribute
    i::Int
end
ParameterBound{Lower}() = ParameterBound{Lower}(-1)
ParameterBound{Upper}() = ParameterBound{Upper}(-1)
ParameterBound() = ParameterBound{Undefined}(-1)

"""
$(TYPEDEF)
A integrator attribute used to access independent variable support set.
"""
struct SupportSet{T <: AbstractFloat} <: AbstractIntegratorAttribute
    s::Vector{T}
end
SupportSet{Float64}() = SupportSet{Float64}(Float64[])
SupportSet() = SupportSet{Float64}()

"""
$(TYPEDEF)
Retreives the LocalIntegrator that corresponds to the integration method used
by the integrator.
"""
struct LocalIntegrator <: AbstractIntegratorAttribute end

"""
$(TYPEDEF)
Retreives the problem used by the integrator.
"""
struct AttachedProblem <: AbstractIntegratorAttribute end

"""
$(FUNCTIONNAME)

Computes the relaxation at the current parameter value with the current
parameter and state bounds.
"""
function relax! end

"""
$(FUNCTIONNAME)

Provides a real-value integration at the current value set for each parameter.
"""
function integrate! end

integrate!(d::AbstractODERelaxIntegrator) = integrate!(d, get(d, AttachedProblem()))


text:
module Calibrate
using Optim
import Optim: optimize

export optimize
export eki_step!
export two_player_eki_step!
export eki_loop!
export calibrate!

include("ensemble_kalman_inversion.jl")

end # module

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
"plugin tools"
module Plugin

using ..Options
using Comonicon
using Pkg

"""
add an ion plugin.

# Arguments

- `name`: plugin name, can be url or a package name.
"""
@cast function add(name::String)
    Pkg.add(name)
    ion = Options.read()::Ion
    push!(ion.plugins, name)
    Options.dump(ion)
    return
end

@cast function rm(name::String)
    Pkg.rm(name)
    ion = Options.read()::Ion
    idx = findfirst(isequal(name), ion.plugins)
    deleteat!(ion.plugins, idx)
    Options.dump(ion)
    return
end

end

@cast Plugin


text:
using Redpitaya
using PyPlot
using FFTW

rp = RedPitaya("10.167.6.139")

freq = 25000.0
dec = optimalDecimation(rp,freq)
#freq=25000
numPeriods = 30
freqR = roundFreq(rp,dec,freq)
numSampPerPeriod = numSamplesPerPeriod(rp,dec,freqR)
numSamp = numSampPerPeriod*numPeriods
if numSamp>16384
#	error("To many sampling points for RedPitaya buffer. Reduce Decimation or reduce number of periods")
end

println("Frequency = $freqR Hz")
println("Number Sampling Points per Period: $numSampPerPeriod")

# start sending
send(rp,"GEN:RST")
#send(rp,"SOUR1:BURS:STAT ON")
#change this to SQUARE to get a SQUARE excitation
sendAnalogSignal(rp,1,"SINE",freqR,0.4,offset=0.2,ncycBurst=60)
#send(rp,"SOUR1:VOLT:OFFS 0")
# receive data
trigger = "AWG_NE" # "CH1_PE"
u1 = receiveAnalogSignalWithTrigger(rp, 1, 0, numSamp, dec=dec, delay=0.2, typ="OLD",
                      trigger=trigger, triggerLevel=-0.0, binary=true, raw=false,
                      triggerDelay=numSampPerPeriod)

figure(3)
subplot(2,1,1)
plot(u1)
subplot(2,1,2)
semilogy(abs.(rfft(u1)),"o-b",lw=2)

#Lets check if the data is periodic
figure(4)
plot(1:length(u1),u1,"x-r",lw=2)
plot(length(u1)+1:2*length(u1),u1,"x-b",lw=2)


text:
using Econometrics, StatsBase, DelimitedFiles
using BSON: @load
using BSON: @save

include("lib/transform.jl")
include("lib/create_transformation.jl")

function Transform()
    @load "raw_data.bson" data
    params = data[:,1:nParams]
    statistics = data[:,nParams+1:end]
    info = create_transformation(statistics)
    statistics = transform(statistics, info)
    @save "cooked_data.bson" params statistics
    return info
end


text:
#NOTE: distribute does not work yet on this version

using MPIClusterManagers # to import MPIManager
using Distributed #provides addprocs()

graph = ModelGraph()

n1 = @optinode(graph)
n2 = @optinode(graph)

@variable(n1,0 <= x <= 2)
@variable(n1,0 <= y <= 3)
@variable(n1, z >= 0)
@constraint(n1,x+y+z >= 4)
@objective(n1,Min,y)

@variable(n2,x)
@NLconstraint(n2,ref,exp(x) >= 2)
@variable(n2,z >= 0)
@constraint(n2,z + x >= 4)
@objective(n2,Min,x)

@linkconstraint(graph,n1[:x] == n2[:x])

manager=MPIManager(np=2)
addprocs(manager)

#Setup worker environments
@everywhere using Pkg
@everywhere Pkg.activate((@__DIR__)*"/..")
@everywhere using Plasmo
@everywhere using PipsSolver

#Distribute the graph to workers
julia_workers = collect(values(manager.mpi2j))
remote_references = PipsSolver.distribute(graph,julia_workers,remote_name = :pipsgraph)

# The remote optigraphs can be queried if they are fetched from the other workers
# r1 = fetch(remote_references[1])
# r2 = fetch(remote_references[2])

#Solve with pips-nlp
@mpi_do manager begin
    using MPI
    PipsSolver.pipsnlp_solve(pipsgraph)
end

rank_zero = manager.mpi2j[0] #julia process representing rank 0
solution = fetch(@spawnat(rank_zero, pipsgraph))


text:
# using Pkg
# pkg"activate .."
using TwiliteTimeline
using Documenter

makedocs(;
    modules=[TwiliteTimeline],
    authors="Data Psycho",
    repo="https://github.com/DataPsycho/TwiliteTimeline.jl/blob/{commit}{path}#L{line}",
    sitename="TwiliteTimeline.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DataPsycho.github.io/TwiliteTimeline.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/DataPsycho/TwiliteTimeline.jl",
)


text:
module FastEigen

export FortranArnoldiMethod,
  NativeLanczos,
  NativeEigen,
  NativeLOBPCG,
  NativeArnoldiMethod,
  NativeEigen,
  KrylovMethods,
  ARPACK

export eigs

import KrylovKit,IterativeSolvers,ArnoldiMethod,Arpack,LinearAlgebra

mutable struct FortranArnoldiMethod
  maxiter::Union{Float64,Nothing}
  tol::Union{Float64,Nothing}
  function FortranArnoldiMethod(;maxiter=nothing,tol=nothing)
    instance = new();
    instance.maxiter = maxiter;
    instance.tol = tol;
    return instance;
  end
end

mutable struct NativeLanczos
  maxiter::Union{Float64,Nothing}
  tol::Union{Float64,Nothing}
  function NativeLanczos(;maxiter=nothing,tol=nothing)
    instance = new();
    instance.maxiter = maxiter;
    instance.tol = tol;
    return instance;
  end
end

mutable struct NativeEigen end

mutable struct NativeArnoldiMethod
  maxiter::Union{Float64,Nothing}
  tol::Union{Float64,Nothing}
  function NativeArnoldiMethod(;maxiter=nothing,tol=nothing)
    instance = new();
    instance.maxiter = maxiter;
    instance.tol = tol;
    return instance;
  end
end

mutable struct NativeLOBPCG
  maxiter::Union{Float64,Nothing}
  tol::Union{Float64,Nothing}
  function NativeLOBPCG(;maxiter=200,tol=nothing)
    instance = new();
    instance.maxiter = maxiter;
    instance.tol = tol;
    return instance;
  end
end

struct KrylovMethods end

struct ARPACK end

function eigs(algorithm::KrylovMethods, x::Array{Float64, 2};
              nev::Int64      = 10,
              order::Symbol   = :LR,
              symmetric::Bool = true,
              maxiter::Int64  = 100)
  local λ, X
  if nev <= 20
    λ, X = KrylovKit.eigsolve(x, nev, order;issymmetric = symmetric)
  else
    λ, X = KrylovKit.eigsolve(x, nev, order;issymmetric = symmetric, krylovdim = nev + 5, maxiter = maxiter)
  end
  X = hcat(X...)
  return λ[1:nev], X[:, 1:nev]
end

function eigs(algorithm::ARPACK, x::Array{Float64, 2};
              nev::Int64      = 10,
              order::Symbol   = :LM,
              symmetric::Bool = true,
              maxiter::Int64  = 100,
              tol = 0.0)
  local λ, X
  λ, X = Arpack.eigs(x, nev = nev, which = order, tol = tol, maxiter = maxiter)
  # @show(λ)
  return λ, X
end

function eigs(algorithm::NativeEigen,x::Array{Float64,2};
              nev::Integer=6,eigen_vectors::Bool=true, order::Symbol=:LR)
  local eigen_val_id;
  eigen_decomp = LinearAlgebra.eigen(x);
  if order == :LR
    eigen_val_id = Base.partialsortperm(eigen_decomp.values,1:nev,rev=true);
  elseif order == :SR
    eigen_val_id = Base.partialsortperm(eigen_decomp.values,1:nev,rev=false)
  end
  if eigen_vectors == true
    return eigen_decomp.values[eigen_val_id],eigen_decomp.vectors[:,eigen_val_id]
  end
  return eigen_decomp.values[eigen_val_id];
end

function eigs(algorithm::NativeLOBPCG,x::Array{Float64,2};
              nev::Integer=6,eigen_vectors::Bool=true, order::Symbol=:LR)
  local eigen_decomp;
  if order == :LR
    eigen_decomp = IterativeSolvers.lobpcg(x,true, nev);
  end
  if order == :SR
    eigen_decomp = IterativeSolvers.lobpcg(x,false, nev);
  end
  if eigen_vectors == true
    return eigen_decomp.λ, eigen_decomp.X;
  end
  return eigen_decomp.λ
end

end


text:
"""
    NBodyEquation{N,O}(orbital, operator::NBodyOperator[, factor::NBodyTerm])

Equation for an `orbital`, acted upon by an operator, which may be a
single-particle operator, or an N-body operator, contracted over all
coordinates but one, and optionally multiplied by an
[`NBodyTerm`](@ref), corresponding to overlaps/matrix elements between
other orbitals.
"""
struct NBodyEquation{QO<:OneBodyOperator,O}
    orbital::O
    operator::QO
    factor::NBodyTerm
    NBodyEquation(orbital::O, operator::QO,
                  factor::NBodyTerm=one(NBodyTerm)) where {QO<:OneBodyOperator,O} =
        new{QO,O}(orbital, operator, factor)
end

Base.iszero(::NBodyEquation) = false

Base.:(*)(nbe::NBodyEquation, factor::NBodyTerm) =
    NBodyEquation(nbe.orbital, nbe.operator, nbe.factor*factor)

function Base.show(io::IO, eq::NBodyEquation)
    show(io, eq.factor, show_sign=true)
    eq.orbital isa Conjugate && write(io, "⟨$(conj(eq.orbital))|")
    show(io, eq.operator)
    !(eq.orbital isa Conjugate) && write(io, "|$(eq.orbital)⟩")
end

"""
    LinearCombinationEquation(equations)

A type representing a linear combination of
[`NBodyEquation`](@ref)s. Typically arises when varying a multi-term
energy expression.
"""
struct LinearCombinationEquation
    equations::Vector{<:NBodyEquation}
end

Base.zero(::Type{LinearCombinationEquation}) = LinearCombinationEquation(NBodyEquation[])
Base.zero(::LinearCombinationEquation) = zero(LinearCombinationEquation)
Base.iszero(eq::LinearCombinationEquation) =
    isempty(eq.equations)

function Base.show(io::IO, eq::LinearCombinationEquation)
    if iszero(eq)
        write(io, "0")
        return
    end
    noprintterms = 2

    terms = if get(io, :limit, false) && length(eq.equations) > 2noprintterms
        vcat(eq.equations[1:noprintterms], "…", eq.equations[end-(noprintterms-1):end])
    else
        eq.equations
    end
    write(io, join(string.(terms), " "))
end

export NBodyEquation, LinearCombinationEquation


text:
using Plots, Random, Measures

function path(rng, alpha, n=5000)
    x, y = 0.0, 0.0
    xDat, yDat = [], []
    for _ in 1:n
        flip = rand(rng,1:4)
        if flip == 1
            x += 1
        elseif flip == 2
            y += 1
        elseif flip == 3
            x -= (2+alpha)*rand(rng)
        elseif flip == 4
            y -= (2+alpha)*rand(rng)
        end
        push!(xDat,x)
        push!(yDat,y)
    end
    return xDat, yDat
end

alphaRange = [0.2, 0.21, 0.22]

default(xlabel = "x", ylabel = "y", xlims=(-150,50), ylims=(-250,50))
p1 = plot(path(MersenneTwister(27), alphaRange[1]), c=:blue)
p1 = plot!(path(MersenneTwister(27), alphaRange[2]), c=:red)
p1 = plot!(path(MersenneTwister(27), alphaRange[3]), c=:green)

rng = MersenneTwister(27)
p2 = plot(path(rng, alphaRange[1]), c=:blue)
p2 = plot!(path(rng, alphaRange[2]), c=:red)
p2 = plot!(path(rng, alphaRange[3]), c=:green)

plot(p1, p2, size=(800, 400), legend=:none, margin=5mm)


text:
function [X,J,dXdx,theta]=JacobiSphereDG(DG,F,Grid)
Rad=Grid.Rad;
ksi=DG.xw;
eta=DG.xw;
n=DG.OrdPoly+1;
lam=zeros(n,n);
theta=zeros(n,n);
X=zeros(n,n,3);
for j=1:n
  ksi2=eta(j);
  for i=1:n
    ksi1=ksi(i);
    X1=0.25*(F.P(1,1).*(1-ksi1)*(1-ksi2)+...
      F.P(1,2).*(1+ksi1)*(1-ksi2)+...
      F.P(1,3).*(1+ksi1)*(1+ksi2)+...
      F.P(1,4).*(1-ksi1)*(1+ksi2));
    X2=0.25*(F.P(2,1).*(1-ksi1)*(1-ksi2)+...
      F.P(2,2).*(1+ksi1)*(1-ksi2)+...
      F.P(2,3).*(1+ksi1)*(1+ksi2)+...
      F.P(2,4).*(1-ksi1)*(1+ksi2));
    X3=0.25*(F.P(3,1).*(1-ksi1)*(1-ksi2)+...
      F.P(3,2).*(1+ksi1)*(1-ksi2)+...
      F.P(3,3).*(1+ksi1)*(1+ksi2)+...
      F.P(3,4).*(1-ksi1)*(1+ksi2));
    X(i,j,1)=X1;
    X(i,j,2)=X2;
    X(i,j,3)=X3;
    r=sqrt(X1^2+X2^2+X3^2);
    f=Rad/r;
    X1=X1/r;
    X2=X2/r;
    X3=X3/r;
    [lam(i,j),theta(i,j)]=cart2sphere(X1,X2,X3);

    % DD=[-sin(lam) cos(lam) 0
    %   0        0     1];
    %
    % sinlam=sin(lam);
    % coslam=cos(lam);
    % sinth=sin(theta);
    % costh=cos(theta);
    % a11=sinlam*sinlam*costh*costh+sinth*sinth;
    % a12=-sinlam*coslam*costh*costh;
    % a13=-coslam*sinth*costh;
    % a21=a12;
    % a22=coslam*coslam*costh*costh+sinth*sinth;
    % a23=-sinlam*sinth*costh;
    % a31=-coslam*sinth;
    % a32=-sinlam*sinth;
    % a33=costh;
    % A=[a11 a12 a13
    %    a21 a22 a23
    %    a31 a32 a33];
    %
    % B=[Grid.Nodes(F.N(1)).P(1) Grid.Nodes(F.N(2)).P(1) Grid.Nodes(F.N(3)).P(1) Grid.Nodes(F.N(4)).P(1)
    %    Grid.Nodes(F.N(1)).P(2) Grid.Nodes(F.N(2)).P(2) Grid.Nodes(F.N(3)).P(2) Grid.Nodes(F.N(4)).P(2)
    %    Grid.Nodes(F.N(1)).P(3) Grid.Nodes(F.N(2)).P(3) Grid.Nodes(F.N(3)).P(3) Grid.Nodes(F.N(4)).P(3)];
    % C=0.25*[-1+ksi2  -1+ksi1
    %          1-ksi2  -1-ksi1
    %          1+ksi2   1+ksi1
    %         -1-ksi2   1-ksi1];
    % J=f*DD*A*B*C;
    % D=abs(det(J));
    % X=[X1 X2 X3]*Rad;
  end
end
minlam=min(min(lam));
maxlam=max(max(lam));
if maxlam-minlam>pi
  for i=1:n
    for j=1:n
      if lam(i,j)>pi
        lam(i,j)=lam(i,j)-pi;
      else
        lam(i,j)=lam(i,j)+pi;
      end
    end
  end
end
lam=Rad*lam; %.*cos(theta);
costheta=cos(theta);
theta=theta*Rad;
dXdx=zeros(n,n,2,2);
dXdx(:,:,1,1)=(DG.DS*lam).*costheta;
dXdx(:,:,2,1)=DG.DS*theta;
dXdx(:,:,1,2)=(lam*DG.DS').*costheta;
dXdx(:,:,2,2)=theta*DG.DS';
J=dXdx(:,:,1,1).*dXdx(:,:,2,2)-dXdx(:,:,1,2).*dXdx(:,:,2,1);

end


text:
s = """
f(a,b,c,d) = a & b ? c :
    (d | c ? b : a)
"""
msgs = lintstr(s)
@test length(msgs)==2
@test msgs[1].code == :I475
@test msgs[1].variable == "&"
@test contains(msgs[1].message, "bit-wise in a boolean context. (&,|) do not have " *
    "short-circuit behavior")
@test msgs[2].code == :I475
@test msgs[2].variable == "|"
@test contains(msgs[2].message, "bit-wise in a boolean context. (&,|) do not have " *
    "short-circuit behavior")


text:
function _plot(iddata::IdData, x::Matrix{Float64},
               name::Vector{String}, unit::Vector{String},
               line::Symbol; kwargs...)
    t = get_time_vector(iddata)
    nx = size(x, 1)

    yguide = [merge_name_unit(name[i], unit[i]) for i=1:nx]
    xguide = [merge_name_unit("t", iddata.time_unit) for i=1:nx]
    grid_width = Int(ceil(sqrt(nx)))

    p = [(i <= nx - grid_width) ?
         Plots.plot(t, x[i, :],
                    color="black",
                    yguide=yguide[i],
                    label=string("measured ", name[i]),
                    line=line) :
         Plots.plot(t, x[i, :],
                    color="black",
                    xguide=xguide[i],
                    yguide=yguide[i],
                    label=string("measured ", name[i]),
                    line=line)
         for i=1:nx]

    Plots.plot(p..., layout=nx; kwargs...)
end

function _plot_prediction!{R<:Real, S<:Real}(iddata::IdData,
                                             pred::Predictor,
                                             Θ::AbstractVector{R},
                                             y0::AbstractMatrix{S}=pred.y0;
                                             label::String="predicted",
                                             kwargs...)

    y = predict(pred, Θ, y0)
    t = iddata.t_start + iddata.ts*(pred.start-1:pred.start+pred.length-2)

    ny = size(y, 1)

    p = [Plots.plot!(t, y[i, :],
        label=string(label, " ", iddata.output_name[i]))
        for i=1:ny]

    Plots.plot!(p..., layout=ny; kwargs...)
end

merge_name_unit(name::String, unit::String) = isempty(unit) ? name : string(name," (", unit, ")")

"""
    plot_input(iddata::IdData)

Plot each one of the system inputs.
"""
function plot_input(iddata::IdData; kwargs...)
    u = get_input(iddata)
    _plot(iddata, u, iddata.input_name, iddata.input_unit, :steppost; kwargs...)
end

function plot_input(e::PredictorError, args...; kwargs...)
    plot_input(e.iddata, args...; kwargs...)
end

"""
    plot_output(iddata::IdData)

Plot each one of the system outputs.
"""
function plot_output(iddata::IdData; kwargs...)
    y = get_output(iddata)
    _plot(iddata, y, iddata.output_name, iddata.output_unit, :line; kwargs...)
end

function plot_output{R<:Real, S<:Real}(iddata::IdData, pred::Predictor,
                                       Θ_list::Vector{Vector{R}},
                                       y0_list::Vector{Matrix{S}}=fill(pred.y0,
                                                                       length(Θ_list));
                                       label::Vector{String}=fill("predicted",
                                                                  length(Θ_list)),
                                       kwargs...)
    n_model_types = length(Θ_list)

    plot_output(iddata; kwargs...)
    p = _plot_prediction!(iddata, pred, Θ_list[1], y0_list[1], label=label[1]; kwargs...)
    for i=2:n_model_types
        p = _plot_prediction!(iddata, pred, Θ_list[i], y0_list[i], label=label[i]; kwargs...)
    end
    p
end

function plot_output{R<:Real, S<:Real}(iddata::IdData, pred::Predictor,
                                       Θ::Vector{R},
                                       y0::Matrix{S}=pred.y0;
                                       label::String="predicted",
                                       kwargs...)
    plot_output(iddata; kwargs...)
    _plot_prediction!(iddata, pred, Θ, y0, label=label; kwargs...)
end

function plot_output(e::PredictorError, args...; kwargs...)
    plot_output(e.iddata, e.pred, args...; kwargs...)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

_default_min_ess(samples::DensitySampleVector) = minimum(unshaped(bat_eff_sample_size(samples).result))

function dist_samples_pvalue(
    rng::AbstractRNG, dist::Distribution, samples::DensitySampleVector;
    nsamples::Integer = floor(Int, _default_min_ess(samples)),
    ess::Integer = floor(Int, _default_min_ess(samples)),
)
    samples_v = bat_sample(rng, samples, OrderedResampling(nsamples = ess)).result.v
    samples_dist_logpdfs = logpdf.(Ref(dist), samples_v)
    ref_samples = bat_sample(rng, dist, IIDSampling(nsamples = nsamples)).result
    ref_dist_logpdfs = ref_samples.logd
    samples_dist_logpdfs, ref_dist_logpdfs

    # KS and AD have trouble with large number of samples on 32-bit systems:
    #HypothesisTests.pvalue(HypothesisTests.ApproximateTwoSampleKSTest(samples_dist_logpdfs, ref_dist_logpdfs))
    #HypothesisTests.pvalue(HypothesisTests.KSampleADTest(Vector(samples_dist_logpdfs), Vector(ref_dist_logpdfs)))
    # So use custom KS-calculation instead:
    ks_pvalue(fast_ks_delta(samples_dist_logpdfs, ref_dist_logpdfs), length(samples_dist_logpdfs), length(ref_dist_logpdfs))
end

function dist_samples_pvalue(dist::Distribution, samples::DensitySampleVector; kwargs...)
    dist_samples_pvalue(bat_rng(), dist, samples; kwargs...)
end

function fast_ks_delta(unsorted_x::AbstractVector{<:Real}, unsorted_y::AbstractVector{<:Real})
    x = sort(unsorted_x)
    y = sort(unsorted_y)
    idxs_x = eachindex(x)
    idxs_y = eachindex(y)
    n_x = length(idxs_x)
    n_y = length(idxs_y)
    i0_x::Int = first(idxs_x)
    i0_y::Int = first(idxs_y)
    i_x::Int = i0_x
    i_y::Int = i0_y
    max_delta::Float64 = 0.0
    while i_x in idxs_x && i_y in idxs_y
        x_ge_y = x[i_x] >= y[i_y]
        delta = abs((i_x - i0_x) / n_x - (i_y - i0_y) / n_y)
        max_delta = max(delta, max_delta)
        i_x = ifelse(!x_ge_y, i_x + 1, i_x)
        i_y = ifelse(x_ge_y, i_y + 1, i_y)
    end
    return max_delta
end

function ks_pvalue(delta::Real, n_x::Integer, n_y::Integer)
    nx, ny = float(n_x), float(n_y)
    n = nx * ny / (nx + ny)
    ccdf(HypothesisTests.Kolmogorov(), sqrt(n) * delta)
end


text:
# Notational note
# JMW's dx_history <=> NW's S
# JMW's dgr_history <=> NW's Y

# Here alpha is a cache that parallels betas
# It is not the step-size
# q is also a cache
function twoloop!(s::Vector,
                  gr::Vector,
                  rho::Vector,
                  dx_history::Matrix,
                  dgr_history::Matrix,
                  m::Integer,
                  pseudo_iteration::Integer,
                  alpha::Vector,
                  q::Vector,
                  precon)
    # Count number of parameters
    n = length(s)

    # Determine lower and upper bounds for loops
    lower = pseudo_iteration - m
    upper = pseudo_iteration - 1

    # Copy gr into q for backward pass
    copy!(q, gr)

    # Backward pass
    for index in upper:-1:lower
        if index < 1
            continue
        end
        i = mod1(index, m)
        @inbounds alpha[i] = rho[i] * vecdot(view(dx_history, :, i), q)
        @simd for j in 1:n
            @inbounds q[j] -= alpha[i] * dgr_history[j, i]
        end
    end

    # Copy q into s for forward pass
    # apply preconditioner if precon != nothing
    # (Note: preconditioner update was done outside of this function)
    A_ldiv_B!(s, precon, q)

    # Forward pass
    for index in lower:1:upper
        if index < 1
            continue
        end
        i = mod1(index, m)
        @inbounds beta = rho[i] * vecdot(view(dgr_history, :, i), s)
        @simd for j in 1:n
            @inbounds s[j] += dx_history[j, i] * (alpha[i] - beta)
        end
    end

    # Negate search direction
    scale!(s, -1)

    return
end

macro lbfgstrace()
    quote
        if tracing
            dt = Dict()
            if o.extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(gr)
                dt["Current step size"] = alpha
            end
            grnorm = vecnorm(gr, Inf)
            update!(tr,
                    iteration,
                    f_x,
                    grnorm,
                    dt,
                    o.store_trace,
                    o.show_trace,
                    o.show_every,
                    o.callback)
        end
    end
end

immutable LBFGS{T} <: Optimizer
    m::Int
    linesearch!::Function
    P::T
    precondprep!::Function
end

LBFGS(; m::Integer = 10, linesearch!::Function = hz_linesearch!,
      P=nothing, precondprep! = (P, x) -> nothing) =
    LBFGS(Int(m), linesearch!, P, precondprep!)

function optimize{T}(d::DifferentiableFunction,
                     initial_x::Vector{T},
                     mo::LBFGS,
                     o::OptimizationOptions)
    # Print header if show_trace is set
    print_header(o)

    # Maintain current state in x and previous state in x_previous
    x, x_previous = copy(initial_x), copy(initial_x)

    # Count the total number of iterations
    iteration = 0
    pseudo_iteration = 0

    # Track calls to function and gradient
    f_calls, g_calls = 0, 0

    # Count number of parameters
    n = length(x)

    # Maintain current gradient in gr and previous gradient in gr_previous
    gr, gr_previous = Array(T, n), Array(T, n)

    # Store a history of changes in position and gradient
    rho = Array(T, mo.m)
    dx_history, dgr_history = Array(T, n, mo.m), Array(T, n, mo.m)

    # The current search direction
    s = Array(T, n)

    # Buffers for use in line search
    x_ls, g_ls = Array(T, n), Array(T, n)

    # Store f(x) in f_x
    f_x_previous, f_x = NaN, d.fg!(x, gr)
    f_calls, g_calls = f_calls + 1, g_calls + 1
    copy!(gr_previous, gr)

    # Keep track of step-sizes
    alpha = alphainit(one(T), x, gr, f_x)

    # TODO: How should this flag be set?
    mayterminate = false

    # Maintain a cache for line search results
    lsr = LineSearchResults(T)

    # Buffers for new entries of dx_history and dgr_history
    dx, dgr = Array(T, n), Array(T, n)

    # Buffers for use by twoloop!
    twoloop_q, twoloop_alpha = Array(T, n), Array(T, mo.m)

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = o.store_trace || o.show_trace || o.extended_trace || o.callback != nothing
    @lbfgstrace

    # Assess multiple types of convergence
    x_converged, f_converged = false, false
    g_converged = vecnorm(gr, Inf) < o.g_tol

    # Iterate until convergence
    converged = g_converged
    while !converged && iteration < o.iterations
        # Increment the number of steps we've had to perform
        iteration += 1
        pseudo_iteration += 1

        # update the preconditioner
        mo.precondprep!(mo.P, x)

        # Determine the L-BFGS search direction
        twoloop!(s, gr, rho, dx_history, dgr_history, mo.m, pseudo_iteration,
                 twoloop_alpha, twoloop_q, mo.P)

        # Refresh the line search cache
        dphi0 = vecdot(gr, s)
        if dphi0 > 0.0
            pseudo_iteration = 1
            @simd for i in 1:n
                @inbounds s[i] = -gr[i]
            end
            dphi0 = vecdot(gr, s)
        end

        clear!(lsr)
        push!(lsr, zero(T), f_x, dphi0)

        # Determine the distance of movement along the search line
        alpha, f_update, g_update =
          mo.linesearch!(d, x, s, x_ls, g_ls, lsr, alpha, mayterminate)
        f_calls, g_calls = f_calls + f_update, g_calls + g_update

        # Maintain a record of previous position
        copy!(x_previous, x)

        # Update current position
        @simd for i in 1:n
            @inbounds dx[i] = alpha * s[i]
            @inbounds x[i] = x[i] + dx[i]
        end

        # Maintain a record of the previous gradient
        copy!(gr_previous, gr)

        # Update the function value and gradient
        f_x_previous = f_x
        f_x = d.fg!(x, gr)
        f_calls, g_calls = f_calls + 1, g_calls + 1

        # Measure the change in the gradient
        @simd for i in 1:n
            @inbounds dgr[i] = gr[i] - gr_previous[i]
        end

        # Update the L-BFGS history of positions and gradients
        rho_iteration = 1 / vecdot(dx, dgr)
        if isinf(rho_iteration)
            # TODO: Introduce a formal error? There was a warning here previously
            break
        end
        dx_history[:, mod1(pseudo_iteration, mo.m)] = dx
        dgr_history[:, mod1(pseudo_iteration, mo.m)] = dgr
        rho[mod1(pseudo_iteration, mo.m)] = rho_iteration

        x_converged,
        f_converged,
        g_converged,
        converged = assess_convergence(x,
                                       x_previous,
                                       f_x,
                                       f_x_previous,
                                       gr,
                                       o.x_tol,
                                       o.f_tol,
                                       o.g_tol)

        @lbfgstrace
    end

    return MultivariateOptimizationResults("L-BFGS",
                                           initial_x,
                                           x,
                                           Float64(f_x),
                                           iteration,
                                           iteration == o.iterations,
                                           x_converged,
                                           o.x_tol,
                                           f_converged,
                                           o.f_tol,
                                           g_converged,
                                           o.g_tol,
                                           tr,
                                           f_calls,
                                           g_calls)
end


text:
using Ion
using Test

@testset "release" begin
    include("commands/release.jl")
end


text:
"""
    complete(g::ADiGraph)

Returns a digraph containing both the edges `(u, v)`
of `g` and their reverse `(v, u)`. See also [`complete!`](@ref).
"""
function complete(g::ADiGraph)
    h = copy(g)
    complete!(h)
    return h
end

"""
    complete!(g::ADiGraph)

For each edge `(u, v)` in `g`, adds to `g` its reverse, i.e. `(v, u)`.
"""
function complete!(g::ADiGraph)
    edgs  = collect(edges(g))
    for e in edgs
        add_edge!(g, dst(e), src(e))
    end
    return g
end

"""
    complement(g)

Produces the [graph complement](https://en.wikipedia.org/wiki/Complement_graph)
of a graph.
"""
function complement(g::T) where T<:AGraph
    gnv = nv(g)
    h = T(gnv)
    for i=1:gnv
        for j=i+1:gnv
            if !has_edge(g, i, j)
                add_edge!(h,i,j)
            end
        end
    end
    return h
end

function complement(g::T) where T<:ADiGraph
    gnv = nv(g)
    h = T(gnv)
    for i=1:gnv
        for j=1:gnv
            if i != j && !has_edge(g,i,j)
                add_edge!(h,i,j)
            end
        end
    end
    return h
end

"""
    union(g, h)

Merges graphs `g` and `h` by taking the set union of all vertices and edges.
"""
function union(g::T, h::T) where T<:AGraphOrDiGraph
    gnv = nv(g)
    hnv = nv(h)

    r = copy(g)
    add_vertices!(r, max(gnv, hnv) - gnv)
    for e in edges(h)
        add_edge!(r, e)
    end
    return r
end

"""
    blockdiag(g, h)

Produces a graph with ``|V(g)| + |V(h)|`` vertices and ``|E(g)| + |E(h)|``
edges.

Put simply, the vertices and edges from graph `h` are appended to graph `g`.
"""
function blockdiag(g::T, h::T) where T<:AGraphOrDiGraph
    gnv = nv(g)
    r = T(gnv + nv(h))
    for e in edges(g)
        add_edge!(r, e)
    end
    for e in edges(h)
        add_edge!(r, gnv+src(e), gnv+dst(e))
    end
    return r
end

"""
    intersect(g, h)

Produces a graph with edges that are only in both graph `g` and graph `h`.

Note that this function may produce a graph with 0-degree vertices.
"""
function intersect(g::T, h::T) where T<:AGraphOrDiGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(min(gnv, hnv))
    for e in intersect(edges(g),edges(h))
        add_edge!(r,e)
    end
    return r
end

"""
    difference(g, h)

Produces a graph with all the edges in graph `g` that are not in graph `h`.

Note that this function may produce a graph with 0-degree vertices.
"""
function difference(g::G, h::T) where {T<:AGraphOrDiGraph,
                                       G<:AGraphOrDiGraph}

    @assert (is_directed(g) && is_directed(h)) ||
            (!is_directed(g) && !is_directed(h))

    gnv = nv(g)
    hnv = nv(h)

    r = G(gnv)
    for e in edges(g)
        !has_edge(h, e) && add_edge!(r,e)
    end
    return r
end

"""
    symmetric_difference(g, h)

Produces a graph with edges from graph `g` that do not exist in graph `h`,
and vice versa.

Note that this function may produce a graph with 0-degree vertices.
"""
function symmetric_difference(g::T, h::T) where T<:AGraphOrDiGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(max(gnv, hnv))
    for e in edges(g)
        !has_edge(h, e) && add_edge!(r, e)
    end
    for e in edges(h)
        !has_edge(g, e) && add_edge!(r, e)
    end
    return r
end

"""
    join(g, h)

Merges graphs `g` and `h` using `blockdiag` and then adds all the edges between
 the vertices in `g` and those in `h`.
"""
function join(g::AGraph, h::AGraph)
    r = blockdiag(g, h)
    for i=1:nv(g)
        for j=nv(g)+1:nv(g)+nv(h)
            add_edge!(r, i, j)
        end
    end
    return r
end

"""
    crosspath(g::AGraph, n::Integer)

Replicate `n` times `g` and connect each vertex with its copies in a path.
"""
crosspath(g::G, n) where {G<:AGraph} = cartesian_product(g, PathGraph(n, G))

"""
    cartesian_product(g, h)

Returns the (cartesian product)[https://en.wikipedia.org/wiki/Cartesian_product_of_graphs] of `g` and `h`
"""
function cartesian_product(g::G, h::G) where G<:AGraphOrDiGraph
    z = G(nv(g) * nv(h))
    id(i, j) = i + (j-1) * nv(g)

    for j=1:nv(h)
        for (i1, i2) in edges(g)
            add_edge!(z, id(i1,j), id(i2,j))
        end
    end

    for (j1, j2) in edges(h)
        for i=1:nv(g)
            add_edge!(z, id(i,j1), id(i,j2))
        end
    end
    return z
end

"""
    tensor_product(g, h)

Returns the (tensor product)[https://en.wikipedia.org/wiki/Tensor_product_of_graphs] of `g` and `h`
"""
function tensor_product(g::G, h::G) where G<:AGraphOrDiGraph
    z = G(nv(g)*nv(h))
    id(i, j) = i + (j-1)*nv(g)
    for (j1, j2) in edges(h)
        for (i1, i2) in edges(g)
            add_edge!(z, id(i1, j1), id(i2, j2))
        end
    end
    return z
end

## subgraphs ###

"""
    subgraph(g, vlist) -> sg, vlist
    subgraph(g, elist) -> sg, vlist

Returns the subgraph of `g` induced by the vertices in  `vlist` or by the edges
in `elist`, along with `vlist` itself (a newly created vector for the second method).

The returned graph has `length(vlist)` vertices, with the new vertex `i`
corresponding to the vertex of the original graph in the `i`-th position
of `vlist`.

For easy subgraph creation also `g[vlist]` or `g[elist]` can be used.

If `g` is a network, vector and edge properties won't be converved
`sg`. You can preserve properties using the [`subnetwork`](@ref) method.

### Usage Examples:
```julia
g = CompleteGraph(10)
sg, vlist = subgraph(g, 5:8)
@assert g[5:8] == sg
@assert nv(sg) == 4
@assert ne(sg) == 6
@assert vm[4] == 8

sg, vlist = subgraph(g, [2,8,3,4])
@asssert sg == g[[2,8,3,4]]

elist = [Edge(1,2), Edge(3,4), Edge(4,8)]
sg, vlist = subgraph(g, elist)
@asssert sg == g[elist]
```
"""
function subgraph(g::G, vlist::AbstractVector{V}) where {G<:AGraphOrDiGraph,V<:Integer}
    allunique(vlist) || error("Vertices in subgraph list must be unique")
    h = G(length(vlist))
    newvid = Dict{V, V}()
    for (i,v) in enumerate(vlist)
        newvid[v] = i
    end

    vset = Set(vlist)
    _build_subraph!(h, g, vset, newvid)

    return h, vlist
end

function _build_subraph!(h::AGraphOrDiGraph, g, vset, newvid)
    for s in vset
        for d in out_neighbors(g, s)
            d in vset && add_edge!(h, newvid[s], newvid[d])
        end
    end
end

function subgraph(g::G, elist) where G<:AGraphOrDiGraph
    h = G()
    V = vertextype(h)
    newvid = Dict{V, V}()
    vlist = Vector{V}()

    for e in elist
        u, v = src(e), dst(e)
        for i in (u,v)
            if !haskey(newvid, i)
                add_vertex!(h)
                newvid[i] = nv(h)
                push!(vlist, i)
            end
        end
        add_edge!(h, newvid[u], newvid[v])
    end
    return h, vlist
end

"""
    subnetwork(g, vlist) -> sg, vlist
    subnetwork(g, elist) -> sg, vlist

Equivalent to [`subgraph`](@ref) but preserves vertex and edge properties
when `g` is a network.
"""
function subnetwork(g::G, vlist::AbstractVector{V}) where {G<:ANetOrDiNet,V<:Integer}
    allunique(vlist) || error("Vertices in subgraph list must be unique")
    h = G(length(vlist))
    newvid = Dict{V, V}()
    for (i,v) in enumerate(vlist)
        newvid[v] = i
    end

    vset = Set(vlist)
    _build_subnetwork!(h, g, vset, newvid)

    return h, vlist
end

function _build_subnetwork!(h::ANetOrDiNet, g, vset, newvid)
    #sound right not to copy graph properties of g
    for (name, prop) in vprop(g)
        vprop!(h, name, valtype(prop))
    end
    for (name, prop) in eprop(g)
        eprop!(h, name, valtype(prop))
    end

    for s in vset
        i = newvid[s]
        for (name, prop) in vprop(g)
            vprop(h, name)[i] = prop[s]
        end
        for e in out_edges(g, s)
            d = dst(e)
            if d in vset
                j = newvid[d]
                ok, enew = add_edge!(h, i, j)
                !ok && continue
                for (name, prop) in eprop(g)
                    eprop(h, name)[enew] = prop[e]
                end
            end
        end
    end
end

function subnetwork(g::G, elist) where G<:ANetOrDiNet
    h = G()
    V = vertextype(h)
    newvid = Dict{V, V}()
    vlist = Vector{V}()

    for (name, prop) in vprop(g)
        vprop!(h, name, valtype(prop))
    end
    for (name, prop) in eprop(g)
        eprop!(h, name,  valtype(prop))
    end

    for e in elist
        u, v = src(e), dst(e)
        for s in (u,v)
            if !haskey(newvid, s)
                add_vertex!(h)
                newvid[s] = nv(h)
                push!(vlist, s)
                for (name, prop) in vprop(g)
                    vprop(h, name)[nv(h)] = prop[s]
                end
            end
        end
        ok, enew = add_edge!(h, newvid[u], newvid[v])
        !ok && continue
        for (name, prop) in eprop(g)
            eprop(h, name)[enew] = prop[e]
        end
    end
    return h, vlist
end

"""
    g[iter]

Returns the subgraph induced by the vertex or edge iterable `iter`.
Equivalent to [`subgraph`](@ref)`(g, iter)[1]` or [`subnetwork`](@ref)`(g, iter)[1]`
for networks.
"""
getindex(g::AGraphOrDiGraph, iter) = subnetwork(g, iter)[1]
# in julia 0.5 always gets dispatched to this (julia bug)
subnetwork(g::AGraphOrDiGraph, list) = subgraph(g, list)

"""
    egonet(g, v::Int, d::Int; dir=:out)

Returns the subgraph of `g` induced by the neighbors of `v` up to distance
`d`. If `g` is a `DiGraph` the `dir` optional argument specifies
the edge direction the edge direction with respect to `v` (i.e. `:in` or `:out`)
to be considered. This is equivalent to [`subgraph`](@ref)`(g, neighborhood(g, v, d, dir=dir))[1].`
"""
egonet(g::AGraphOrDiGraph, v::Integer, d::Integer; dir=:out) =  g[neighborhood(g, v, d, dir=dir)]

# The following operators allow one to use a Erdos.Graph as a matrix in
# eigensolvers for spectral ranking and partitioning.
# """Provides multiplication of a graph `g` by a vector `v` such that spectral
# graph functions in [GraphMatrices.jl](https://github.com/jpfairbanks/GraphMatrices.jl) can utilize Erdos natively.
# """
function *(g::AGraph, v::Vector{T}) where T<:Number
    length(v) == nv(g) || error("Vector size must equal number of vertices")
    y = zeros(T, nv(g))
    for e in edges(g)
        i = src(e)
        j = dst(e)
        y[i] += v[j]
        y[j] += v[i]
    end
    return y
end

function *(g::ADiGraph, v::Vector{T}) where T<:Number
    length(v) == nv(g) || error("Vector size must equal number of vertices")
    y = zeros(T, nv(g))
    for e in edges(g)
        i = src(e)
        j = dst(e)
        y[i] += v[j]
    end
    return y
end

"""sum(g,i) provides 1:in_degree or 2:out_degree vectors"""
function sum(g::AGraphOrDiGraph, dim::Int)
    dim == 1 && return in_degree(g, vertices(g))
    dim == 2 && return out_degree(g, vertices(g))
    error("Graphs are only two dimensional")
end

size(g::AGraphOrDiGraph) = (nv(g), nv(g))

"""`size(g,i)` provides 1:nv or 2:nv else 1 """
size(g::AGraph,dim::Int) = (dim == 1 || dim == 2) ? nv(g) : 1

"""`sum(g)`` provides the number of edges in the graph"""
sum(g::AGraphOrDiGraph) = ne(g)

"""
    sparse(g)

Equivalent to [`adjacency_matrix`](@ref).
"""
sparse(g::AGraphOrDiGraph) = adjacency_matrix(g)

#arrayfunctions = (:eltype, :length, :ndims, :size, :strides, :issymmetric)
eltype(g::AGraphOrDiGraph) = vertextype(g)
length(g::AGraphOrDiGraph) = nv(g)*nv(g)
ndims(g::AGraphOrDiGraph) = 2
issymmetric(g::AGraphOrDiGraph) = !is_directed(g)

"""
    contract!(g, vs)
    contract!(g, v1, v2, ....)

Merge the vertices in `vs` into a unique vertex.
"""
contract!(g::AGraphOrDiGraph, vs::AbstractVector) = contract!(g, Set(vs))
contract!(g::AGraphOrDiGraph, vs::Integer...) = contract!(g, Set(vs))

function contract!(g::AGraph, vs::Set)
    n = add_vertex!(g)
    for v in vs
        for u in neighbors(g, v)
            if u ∉ vs
                add_edge!(g, n, u)
            end
        end
    end
    return rem_vertices!(g, vs)
end

function contract!(g::ADiGraph, vs::Set)
    n = add_vertex!(g)
    for v in vs
        for u in out_neighbors(g, v)
            if u ∉ vs
                add_edge!(g, n, u)
            end
        end
        for u in in_neighbors(g, v)
            if u ∉ vs
                add_edge!(g, u, n)
            end
        end
    end
    return rem_vertices!(g, vs)
end


text:
__precompile__(true)

"""
    RootsAndPoles.jl

`RootsAndPoles` is a Julia implementation of the Global complex Roots and Poles Finding
(GRPF) algorithm.

Matlab code is available under MIT license at https://github.com/PioKow/GRPF.

# References

[^1]: P. Kowalczyk, “Global complex roots and poles finding algorithm based on phase
    analysis for propagation and radiation problems,” IEEE Transactions on Antennas and
    Propagation, vol. 66, no. 12, pp. 7198–7205, Dec. 2018, doi: 10.1109/TAP.2018.2869213.
"""
module RootsAndPoles

#==
NOTE: Some variable conversions from the original GRPF papers to this code:

| Paper | Code |
|-------|------|
|   𝓔   |   E  |
|   𝐶   |   C  |
|   ϕ   |  phi |
==#

import Base
import Base.Threads.@threads
using LinearAlgebra
using VoronoiDelaunay
import VoronoiDelaunay: getx, gety, DelaunayEdge, DelaunayTriangle

# NOTE: `max_coord` and `min_coord` are provided by `VoronoiDelaunay.jl`
# We are even more conservative than going from `max_coord` to `min_coords` because it is
# possible to run into floating point issues at the very limits
const MAXCOORD = nextfloat(max_coord, -10)
const MINCOORD = nextfloat(min_coord, 10)

"""
    GRPFParams

Mutable struct for holding values used by `RootsAndPoles` to stop iterating or split
Delaunay triangles.

Default values are provided by `GRPFParams()`.

# Fields

- `maxiterations::Int = 100`: the maximum number of refinement iterations before `grpf`
    returns.
- `maxnodes::Int = 500000`: the maximum number of Delaunay tessalation nodes before `grpf`
    returns.
- `skinnytriangle::Int = 3`: maximum ratio of the longest to shortest side length of
    Delaunay triangles before they are split during `grpf` refinement iterations.
- `tess_sizehint::Int = 5000`: provide a size hint to the total number of expected nodes in
    the Delaunay tesselation. Setting this number approximately correct can improve
    performance.
- `tolerance::Float64 = 1e-9`: maximum allowed edge length of the tesselation defined in the
    `origcoords` domain before returning.
- `multithreading::Bool = false`: use `Threads.@threads` to run the user-provided function
    `fcn` across the `DelaunayTriangulation`.
"""
mutable struct GRPFParams
    maxiterations::Int
    maxnodes::Int
    skinnytriangle::Int
    tess_sizehint::Int
    tolerance::Float64
    multithreading::Bool

    function GRPFParams(maxiterations, maxnodes, skinnytriangle, tess_sizehint, tolerance,
                        multithreading)
        tess_sizehint > maxnodes && @warn "GRPFParams `tess_sizehint` is greater than `maxnodes`"
        new(maxiterations, maxnodes, skinnytriangle, tess_sizehint, tolerance, multithreading)
    end
end
GRPFParams() = GRPFParams(100, 500000, 3, 5000, 1e-9, false)

"""
    GRPFParams(tess_sizehint, tolerance, multithreading=false)

Convenience function for creating a `GRPFParams` object with the most important parameters,
`tess_sizehint`, `tolerance`, and `multithreading`.
"""
GRPFParams(tess_sizehint, tolerance, multithreading=false) =
    GRPFParams(100, 500000, 3, tess_sizehint, tolerance, multithreading)

function Base.isequal(a::GRPFParams, b::GRPFParams)
    for n in fieldnames(GRPFParams)
        isequal(getfield(a,n), getfield(b,n)) || return false
    end
    return true
end
Base.:(==)(a::GRPFParams, b::GRPFParams) = isequal(a,b)

struct PlotData end

"""
    Geometry2Function{T}

Store conversion coefficients from the `VoronoiDelaunay` domain to the `origcoords`
function domain.

`Geometry2Function` instances can be called as a function, e.g.
`z_functiondomain = g2f(z_voronoidelaunay)`.

`Geometry2Function` structs are created internally to `RootsAndPoles` with the
appropriate scaling parameters but are returned when `grpf` is called with the `PlotData()`
argument.
"""
struct Geometry2Function{T}
    rmin::T
    rmax::T
    imin::T
    imax::T
end
(f::Geometry2Function)(z) = geom2fcn(z, f.rmin, f.rmax, f.imin, f.imax)
(f::Geometry2Function)(x, y) = geom2fcn(x, y, f.rmin, f.rmax, f.imin, f.imax)
Base.eltype(f::Geometry2Function{T}) where T = T

"""
    ScaledFunction{F,G<:Geometry2Function}

Store conversion coefficients from the `VoronoiDelaunay` domain to the `origcoords` domain
so that a `ScaledFunction` can be called in place of the original function when providing an
argument in the `VoronoiDelaunay` domain.
"""
struct ScaledFunction{F,G<:Geometry2Function}
    fcn::F
    g2f::G
end
(f::ScaledFunction)(z) = f.fcn(f.g2f(z))

# These files need the above structs defined
include("VoronoiDelaunayExtensions.jl")
include("utils.jl")
include("coordinate_domains.jl")

export rectangulardomain, diskdomain, grpf, PlotData, getplotdata, GRPFParams

"""
    quadrant(val)::Int8

Convert complex function value `val` to quadrant number.

| Quadrant |       Phase       |
|:--------:|:-----------------:|
|    1     | 0 ≤ arg f < π/2   |
|    2     | π/2 ≤ arg f < π   |
|    3     | π ≤ arg f < 3π/2  |
|    4     | 3π/2 ≤ arg f < 2π |
"""
@inline function quadrant(val)::Int8
    # This function correponds to `vinq.m`

    rv, iv = reim(val)
    if (rv > 0) & (iv >= 0)
        return 1
    elseif (rv <= 0) & (iv > 0)
        return 2
    elseif (rv < 0) & (iv <= 0)
        return 3
    else
        # (rv >= 0) & (iv < 0)
        return 4
    end
end

"""
    assignquadrants!(quadrants, nodes, f, multithreading=false)

Evaluate function `f` for [`quadrant`](@ref) at `nodes` and fill `quadrants` in-place.

Each element of `quadrants` corresponds to the index of `IndexablePoint2D` in `nodes`.
"""
function assignquadrants!(
    quadrants::Vector{<:Integer},
    nodes::Vector{IndexablePoint2D},
    f::ScaledFunction{F,G},
    multithreading=false
    ) where {F,G}

    if multithreading
        @threads for ii in eachindex(nodes)
            p = @inbounds nodes[ii]
            quadrants[getindex(p)] = quadrant(f(p))
        end
    else
        for ii in eachindex(nodes)
            p = @inbounds nodes[ii]
            quadrants[getindex(p)] = quadrant(f(p))
        end
    end
    return nothing
end

"""
    candidateedges!(E, tess, quadrants)

Fill in candidate edges `E` that contain a phase change of 2 quadrants.

Any root or pole is located at the point where the regions described by four different
quadrants meet. Since any triangulation of the four nodes located in the four different
quadrants requires at least one edge of ``|ΔQ| = 2``, then all such edges are potentially in
the vicinity of a root or pole.

`E` is not sorted.
"""
function candidateedges!(
    E::Vector{DelaunayEdge{IndexablePoint2D}},
    tess::DelaunayTessellation2D{IndexablePoint2D},
    quadrants::Vector{<:Integer}
    )

    # Better performance compared to `delaunayedges()`
    # see: https://github.com/JuliaGeometry/VoronoiDelaunay.jl/issues/47
    @inbounds for ix in 2:tess._last_trig_index
        tr = tess._trigs[ix]
        isexternal(tr) && continue

        # precalculate
        atr, btr, ctr = geta(tr), getb(tr), getc(tr)
        atri, btri, ctri = getindex(atr), getindex(btr), getindex(ctr)

        ix_na = tr._neighbour_a
        if (ix_na > ix) || isexternal(tess._trigs[ix_na])
            ΔQ = mod(quadrants[btri] - quadrants[ctri], 4)  # phase difference
            if ΔQ == 2
                edge = DelaunayEdge(btr, ctr)
                push!(E, edge)
            end
        end

        ix_nb = tr._neighbour_b
        if (ix_nb > ix) || isexternal(tess._trigs[ix_nb])
            ΔQ = mod(quadrants[atri] - quadrants[ctri], 4)
            if ΔQ == 2
                edge = DelaunayEdge(atr, ctr)
                push!(E, edge)
            end
        end

        ix_nc = tr._neighbour_c
        if (ix_nc > ix) || isexternal(tess._trigs[ix_nc])
            ΔQ = mod(quadrants[atri] - quadrants[btri], 4)
            if ΔQ == 2
                edge = DelaunayEdge(atr, btr)
                push!(E, edge)
            end
        end
    end

    return nothing
end

"""
    candidateedges!(E, tess, quadrants, ::PlotData)

Return candidate edges `E` and the `phasediffs` as a `Vector{Int8}` across each edge.
"""
function candidateedges!(
    E::Vector{DelaunayEdge{IndexablePoint2D}},
    tess::DelaunayTessellation2D{IndexablePoint2D},
    quadrants::Vector{<:Integer},
    ::PlotData
    )

    phasediffs = Vector{Int8}()

    for edge in delaunayedges_fast(tess)
        nodea, nodeb = geta(edge), getb(edge)
        idxa, idxb = getindex(nodea), getindex(nodeb)

        # NOTE: To match Matlab, force `idxa` < `idxb`
        # (order doesn't matter for `ΔQ == 2`, which is the only case we care about)
        # if idxa > idxb
        #     idxa, idxb = idxb, idxa
        # end

        @inbounds ΔQ = mod(quadrants[idxa] - quadrants[idxb], Int8(4))  # phase difference
        if ΔQ == 2
            push!(E, edge)
        end

        push!(phasediffs, ΔQ)
    end

    return E, phasediffs
end

"""
    zone(triangle, edge_idxs)

Return zone `1` or `2` for `triangle`.

Zone `1` triangles have more than one node in `edge_idxs`, whereas zone `2` triangles have
only a single node.
"""
function zone(
    triangle::DelaunayTriangle{IndexablePoint2D},
    edge_idxs::Vector{<:Integer}
    )

    na = geta(triangle)
    nb = getb(triangle)
    nc = getc(triangle)

    nai = getindex(na)
    nbi = getindex(nb)
    nci = getindex(nc)

    zone2 = false
    for idx in edge_idxs
        # with Julia 1.4.2: | is faster than || here
        if (nai == idx) | (nbi == idx) | (nci == idx)
            if !zone2
                # we need to keep searching b/c it might be zone 1
                zone2 = true
            else
                # zone1 = true
                return 1
            end
        end
    end

    return zone2 ? 2 : 0
end

"""
    uniqueindices!(idxs, edges)

Compute unique indices `idxs` in-place of all nodes in `edges`.
"""
function uniqueindices!(
    idxs::Vector{<:Integer},
    edges::Vector{DelaunayEdge{IndexablePoint2D}}
    )
    empty!(idxs)
    for i in eachindex(edges)
        @inbounds ei = edges[i]
        push!(idxs, getindex(geta(ei)))
        push!(idxs, getindex(getb(ei)))
    end
    sort!(idxs)  # calling `sort!` first makes `unique!` more efficient
    unique!(idxs)
    return nothing
end

"""
    zone1newnodes!(newnodes, triangles, g2f, tolerance)

Add nodes to `newnodes` in-place if they are in zone 1, i.e. triangles that had more than
one node.
"""
function zone1newnodes!(
    newnodes::Vector{IndexablePoint2D},
    triangles::Vector{DelaunayTriangle{IndexablePoint2D}},
    g2f::Geometry2Function,
    tolerance
    )

    triangle1 = triangles[1]
    n1a = geta(triangle1)
    n1b = getb(triangle1)
    push!(newnodes, (n1a+n1b)/2)

    @inbounds for ii = 1:length(triangles)-1
        triangle = triangles[ii]
        na = geta(triangle)
        nb = getb(triangle)
        nc = getc(triangle)

        addnewnode!(newnodes, nb, nc, g2f, tolerance)
        addnewnode!(newnodes, nc, na, g2f, tolerance)
        addnewnode!(newnodes, geta(triangles[ii+1]), getb(triangles[ii+1]), g2f, tolerance)
    end
    te = triangles[end]
    na = geta(te)
    nb = getb(te)
    nc = getc(te)
    addnewnode!(newnodes, nb, nc, g2f, tolerance)
    addnewnode!(newnodes, nc, na, g2f, tolerance)

    # Remove the first of `newnodes` if the edge is too short
    distance(g2f(n1a), g2f(n1b)) < tolerance && popfirst!(newnodes)
    return nothing
end

@inline function addnewnode!(
    newnodes::Vector{IndexablePoint2D},
    node1::IndexablePoint2D,
    node2::IndexablePoint2D,
    g2f::Geometry2Function,
    tolerance
    )

    if distance(g2f(node1), g2f(node2)) > tolerance
        avgnode = (node1+node2)/2
        @inbounds for ii in eachindex(newnodes)
            distance(newnodes[ii], avgnode) < 2*eps() && return nothing
        end
        push!(newnodes, avgnode)  # only executed if we haven't already returned
    end
    return nothing
end

"""
    zone2newnodes!(newnodes, triangle, skinnytriangle)

Add node to `newnodes` for zone 2 ("skinny") triangles.

`skinnytriangle` is the maximum allowed ratio of the longest to shortest side length.
"""
@inline function zone2newnode!(
    newnodes::Vector{IndexablePoint2D},
    triangle::DelaunayTriangle{IndexablePoint2D},
    skinnytriangle
    )

    na = geta(triangle)
    nb = getb(triangle)
    nc = getc(triangle)

    # For skinny triangle check, `geom2fcn` not needed because units cancel out
    l1 = distance(na, nb)
    l2 = distance(nb, nc)
    l3 = distance(nc, na)
    if max(l1,l2,l3)/min(l1,l2,l3) > skinnytriangle
        avgnode = (na+nb+nc)/3
        push!(newnodes, avgnode)
    end

    return nothing
end

"""
    splittriangles!(zone1triangles, newnodes, tess, edge_idxs, params)

Add zone 2 triangles to `newnodes` and then update vector of zone 1 triangles, which require
special handling.
"""
function splittriangles!(
    zone1triangles::Vector{DelaunayTriangle{IndexablePoint2D}},
    newnodes::AbstractVector,
    tess::DelaunayTessellation2D{IndexablePoint2D},
    edge_idxs::Vector{<:Integer},
    params::GRPFParams
    )

    empty!(zone1triangles)
    for triangle in tess
        z = zone(triangle, edge_idxs)

        if z == 1
            push!(zone1triangles, triangle)
        elseif z == 2
            zone2newnode!(newnodes, triangle, params.skinnytriangle)
        end
    end
    return nothing
end

"""
    findnextnode(prevnode, refnode, tempnodes, g2f)

Find the index of the next node in the candidate region boundary process. The next one
(after the reference) is picked from the fixed set of nodes.
"""
@inline function findnextnode(
    prevnode::IndexablePoint2D,
    refnode::IndexablePoint2D,
    tempnodes::Vector{IndexablePoint2D},
    g2f::Geometry2Function
    )

    P = g2f(prevnode)
    S = g2f(refnode)

    minphi = 2π + 1  # max diff of angles is 2π, so this is guaranteed larger
    minphi_idx = 1

    for i in eachindex(tempnodes)
        @inbounds N = g2f(tempnodes[i])

        SP = P - S
        SN = N - S

        phi = mod2pi(angle(SP) - angle(SN))

        if phi < minphi
            minphi = phi
            minphi_idx = i
        end
    end

    return minphi_idx
end

"""
    contouredges(tess, edges)

Find contour edges from all candidate edges.
"""
function contouredges(
    tess::DelaunayTessellation2D{IndexablePoint2D},
    edges::Vector{DelaunayEdge{IndexablePoint2D}}
    )

    C = Vector{DelaunayEdge{IndexablePoint2D}}()
    sizehint!(C, length(edges))

    # Edges of triangles that contain at least 1 of `edges`
    for triangle in tess
        pa, pb, pc = geta(triangle), getb(triangle), getc(triangle)
        pai, pbi, pci = getindex(pa), getindex(pb), getindex(pc)

        for edge in edges
            eai, ebi = getindex(geta(edge)), getindex(getb(edge))

            if ((eai == pai) && (ebi == pbi)) | ((eai == pbi) && (ebi == pai)) |
                ((eai == pbi) && (ebi == pci)) | ((eai == pci) && (ebi == pbi)) |
                ((eai == pci) && (ebi == pai)) | ((eai == pai) && (ebi == pci))
                push!(C, DelaunayEdge(pa,pb), DelaunayEdge(pb,pc), DelaunayEdge(pc,pa))
                break  # only count each triangle once
            end
        end
    end

    # Remove duplicate edges
    sameunique!(C)

    return C
end

"""
    evaluateregions!(C, g2f)
"""
function evaluateregions!(
    C::Vector{DelaunayEdge{IndexablePoint2D}},
    g2f::Geometry2Function
    )

    # NOTE: The nodes of each region are in reverse order compared to Matlab with respect
    # to their quadrants

    # Initialize
    numregions = 1

    regions = [[geta(C[1])]]

    refnode = getb(C[1])  # type annotated to assist with boxing
    popfirst!(C)

    nextedgeidxs = Vector{Int}()
    while length(C) > 0

        # This loop is equivalent to `findall(e->geta(e)==refnode, C)`
        # but avoids closure Core.Box issue
        @inbounds for i in eachindex(C)
            if geta(C[i]) == refnode
                push!(nextedgeidxs, i)
            end
        end

        if !isempty(nextedgeidxs)
            if length(nextedgeidxs) == 1
                nextedgeidx = only(nextedgeidxs)
            else
                prevnode = regions[numregions][end]
                tempnodes = getb.(C[nextedgeidxs])
                idx = findnextnode(prevnode, refnode, tempnodes, g2f)
                nextedgeidx = nextedgeidxs[idx]
            end

            nextedge = C[nextedgeidx]
            push!(regions[numregions], geta(nextedge))
            refnode = getb(nextedge)
            deleteat!(C, nextedgeidx)
        else # isempty
            push!(regions[numregions], refnode)
            # New region
            numregions += 1
            push!(regions, [geta(C[1])])
            refnode = getb(C[1])
            popfirst!(C)
        end

        # reset `nextedgeidxs`
        empty!(nextedgeidxs)
    end

    push!(regions[numregions], refnode)

    return regions
end

"""
    rootsandpoles(regions, quadrants, g2f)

Identify roots and poles of function based on regions and quadrants.
"""
function rootsandpoles(
    regions::Vector{Vector{IndexablePoint2D}},
    quadrants::Vector{Int8},
    g2f::Geometry2Function{T}
    ) where T

    numregions = length(regions)

    complexT = complex(T)
    zroots = Vector{complexT}()
    zpoles = Vector{complexT}()
    for ii in eachindex(regions)
        quadrantsequence = [quadrants[getindex(node)] for node in regions[ii]]

        # Sign flip because `regions[ii]` are in opposite order of Matlab?
        dQ = -diff(quadrantsequence)
        @inbounds for jj in eachindex(dQ)
            if dQ[jj] == 3
                dQ[jj] = -1
            elseif dQ[jj] == -3
                dQ[jj] = 1
            elseif abs(dQ[jj]) == 2
                # ``|ΔQ| = 2`` is ambiguous; cannot tell whether phase increases or
                # decreases by two quadrants
                dQ[jj] = 0
            end
        end
        q = sum(dQ)/4
        z = sum(g2f.(regions[ii]))/length(regions[ii])

        if q > 0
            push!(zroots, convert(complexT, z))  # convert in case T isn't Float64
        elseif q < 0
            push!(zpoles, convert(complexT, z))
        end
    end

    return zroots, zpoles
end

"""
    tesselate!(tess, newnodes, f, params)

Label quadrants, identify candidate edges, and iteratively split triangles, returning
the tuple `(tess, E, quadrants)`.
"""
function tesselate!(
    tess::DelaunayTessellation2D{IndexablePoint2D},
    newnodes::Vector{IndexablePoint2D},
    f::ScaledFunction,
    params::GRPFParams
    )

    # Initialize
    numnodes = tess._total_points_added
    @assert numnodes == 0

    g2f = f.g2f

    E = Vector{DelaunayEdge{IndexablePoint2D}}()
    quadrants = Vector{Int8}()
    edge_idxs = Vector{Int}()
    zone1triangles = Vector{DelaunayTriangle{IndexablePoint2D}}()

    iteration = 0
    while (iteration < params.maxiterations) && (numnodes < params.maxnodes)
        iteration += 1

        # Determine which quadrant function value belongs at each node
        numnewnodes = length(newnodes)
        append!(quadrants, Vector{Int8}(undef, numnewnodes))
        assignquadrants!(quadrants, newnodes, f, params.multithreading)

        # Add new nodes to `tess`
        push!(tess, newnodes)
        numnodes += numnewnodes

        # Determine candidate edges that may be near a root or pole
        empty!(E)  # start with a blank E
        candidateedges!(E, tess, quadrants)
        isempty(E) && tess, E, quadrants  # no roots or poles found

        # Select candidate edges that are longer than the chosen tolerance
        selectE = filter(e -> longedge(e, params.tolerance, g2f), E)
        isempty(selectE) && return tess, E, quadrants

        # return if maximum edge length has reached tolerance
        maxElength = maximum(distance(g2f(e)) for e in selectE)
        maxElength < params.tolerance && return tess, E, quadrants

        # Get unique indices of nodes in `edges`
        uniqueindices!(edge_idxs, selectE)

        # Refine (split) triangles
        empty!(newnodes)
        splittriangles!(zone1triangles, newnodes, tess, edge_idxs, params)

        # Add new nodes in zone 1
        zone1newnodes!(newnodes, zone1triangles, g2f, params.tolerance)

        # Have to assign indexes to new nodes (which are all currently -1)
        setindex!.(newnodes, (1:length(newnodes)).+numnodes)
    end

    (iteration >= params.maxiterations) && @warn "GRPFParams.maxiterations reached"
    (numnodes >= params.maxnodes) && @warn "GRPFParams.maxnodes reached"

    return tess, E, quadrants
end

function tesselate!(
    tess::DelaunayTessellation2D{IndexablePoint2D},
    newnodes::Vector{IndexablePoint2D},
    f::ScaledFunction,
    params::GRPFParams,
    ::PlotData
    )

    # Initialize
    numnodes = tess._total_points_added
    @assert numnodes == 0

    g2f = f.g2f

    E = Vector{DelaunayEdge{IndexablePoint2D}}()
    quadrants = Vector{Int8}()
    edge_idxs = Vector{Int}()
    zone1triangles = Vector{DelaunayTriangle{IndexablePoint2D}}()

    iteration = 0
    let phasediffs
    while (iteration < params.maxiterations) && (numnodes < params.maxnodes)
        iteration += 1

        # Determine which quadrant function value belongs at each node
        numnewnodes = length(newnodes)
        append!(quadrants, Vector{Int8}(undef, numnewnodes))
        assignquadrants!(quadrants, newnodes, f, params.multithreading)

        # Add new nodes to `tess`
        push!(tess, newnodes)
        numnodes += numnewnodes

        # Determine candidate edges that may be near a root or pole
        empty!(E)  # always start with a blank E
        E, phasediffs = candidateedges!(E, tess, quadrants, PlotData())
        isempty(E) && return tess, E, quadrants, phasediffs  # no roots or poles found

        # Select candidate edges that are longer than the chosen tolerance
        selectE = filter(e -> longedge(e, params.tolerance, g2f), E)
        isempty(selectE) && return tess, E, quadrants, phasediffs

        # return if maximum edge length has reached tolerance
        maxElength = maximum(distance(g2f(e)) for e in selectE)
        maxElength < params.tolerance && return tess, E, quadrants, phasediffs

        # Get unique indices of nodes in `edges`
        uniqueindices!(edge_idxs, selectE)

        # Refine (split) triangles
        empty!(newnodes)
        splittriangles!(zone1triangles, newnodes, tess, edge_idxs, params)

        # Add new nodes in zone 1
        zone1newnodes!(newnodes, zone1triangles, g2f, params.tolerance)

        # Have to assign indexes to new nodes (which are all currently -1)
        setindex!.(newnodes, (1:length(newnodes)).+numnodes)
    end

    (iteration >= params.maxiterations) && @warn "GRPFParams.maxiterations reached"
    (numnodes >= params.maxnodes) && @warn "GRPFParams.maxnodes reached"

    return tess, E, quadrants, phasediffs
    end  # let
end

function setup(origcoords)
    rmin, rmax = minimum(real, origcoords), maximum(real, origcoords)
    imin, imax = minimum(imag, origcoords), maximum(imag, origcoords)

    # Scaling parameters
    g2f = Geometry2Function(rmin, rmax, imin, imax)
    scaledorigcoords = fcn2geom.(origcoords, rmin, rmax, imin, imax)

    @assert minimum(real, scaledorigcoords) >= min_coord &&
        minimum(imag, scaledorigcoords) >= min_coord &&
        maximum(real, scaledorigcoords) <= max_coord &&
        maximum(imag, scaledorigcoords) <= max_coord "Scaled coordinates out of bounds"

    newnodes = [IndexablePoint2D(real(coord), imag(coord), idx) for (idx, coord) in
                enumerate(scaledorigcoords)]

    return newnodes, g2f
end

"""
    grpf(fcn, origcoords, params=GRPFParams())

Return a vector `roots` and a vector `poles` of a single (complex) argument function
`fcn`.

Searches within a domain specified by the vector of complex `origcoords`.

# Examples
```jldoctest
julia> simplefcn(z) = (z - 1)*(z - im)^2*(z + 1)^3/(z + im)

julia> xb, xe = -2, 2

julia> yb, ye = -2, 2

julia> r = 0.1

julia> origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)

julia> roots, poles = grpf(simplefcn, origcoords);

julia> roots
3-element Array{Complex{Float64},1}:
    -0.9999999998508017 - 8.765385802810127e-11im
 1.3587683359414186e-10 + 1.0000000001862643im
     1.0000000002536367 + 1.0339757656912847e-26im

julia> poles
1-element Array{Complex{Float64},1}:
 -2.5363675604239815e-10 - 1.0000000002980232im
```
"""
function grpf(fcn, origcoords::AbstractArray{<:Complex}, params=GRPFParams())
    newnodes, g2f = setup(origcoords)
    f = ScaledFunction(fcn, g2f)

    tess = DelaunayTessellation2D{IndexablePoint2D}(params.tess_sizehint)

    tess, E, quadrants = tesselate!(tess, newnodes, f, params)

    complexT = complex(eltype(g2f))
    isempty(E) && return Vector{complexT}(), Vector{complexT}()

    C = contouredges(tess, E)
    regions = evaluateregions!(C, g2f)
    zroots, zpoles = rootsandpoles(regions, quadrants, g2f)

    return zroots, zpoles
end

"""
    grpf(fcn, origcoords, ::PlotData, params=GRPFParams())

Variant of `grpf` that returns `quadrants`, `phasediffs`, the `VoronoiDelauany`
tesselation `tess`, and the `Geometry2Function` struct `g2f` for converting from the
`VoronoiDelaunay` to function space, in addition to `zroots` and `zpoles`.

These additional outputs are primarily for plotting or diagnostics.

# Examples
```jldoctest
julia> simplefcn(z) = (z - 1)*(z - im)^2*(z + 1)^3/(z + im)

julia> xb, xe = -2, 2

julia> yb, ye = -2, 2

julia> r = 0.1

julia> origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)

julia> roots, poles, quadrants, phasediffs, tess, g2f = grpf(simplefcn, origcoords, PlotData());
```
"""
function grpf(fcn, origcoords::AbstractArray{<:Complex}, ::PlotData, params=GRPFParams())
    newnodes, g2f = setup(origcoords)
    f = ScaledFunction(fcn, g2f)

    tess = DelaunayTessellation2D{IndexablePoint2D}(params.tess_sizehint)

    tess, E, quadrants, phasediffs = tesselate!(tess, newnodes, f, params, PlotData())

    complexT = complex(eltype(g2f))
    isempty(E) && return (Vector{complexT}(), Vector{complexT}(), quadrants, phasediffs,
                          tess, g2f)

    C = contouredges(tess, E)
    regions = evaluateregions!(C, g2f)
    zroots, zpoles = rootsandpoles(regions, quadrants, g2f)

    return zroots, zpoles, quadrants, phasediffs, tess, g2f
end

end # module


text:
function speed_test(repeats)
	Hnat, controlHams, controlFields, controlFreqs = sim_setup(16, 2000, 4)

	# run once for JIT
	evolution_unitary(Hnat, controlHams, controlFields, controlFreqs)
	parallel_evolution_unitary(Hnat, controlHams, controlFields, controlFreqs)

	# then time each
	tserial = @elapsed for ct = 1:repeats
		evolution_unitary(Hnat, controlHams, controlFields, controlFreqs)
	end

	tparallel = @elapsed for ct = 1:repeats
		parallel_evolution_unitary(Hnat, controlHams, controlFields, controlFreqs)
	end

	println("Serial execution in $(tserial/repeats) seconds")
	println("Parallel execution in $(tparallel/repeats) seconds")

	tserial, tparallel
end

function sim_setup(dimension,
                   numTimeSteps,
                   numControls)
    #Create a random natural hamiltonian
    tmpMat = randn(dimension, dimension) + 1im*randn(dimension, dimension)
    Hnat = tmpMat+tmpMat'

    #Create random control Hamiltonians
    controlHams = Array(Matrix{Complex128}, numControls)
    for ct = 1:numControls
        tmpMat[:] = randn(dimension, dimension) + 1im*randn(dimension, dimension)
        controlHams[ct] = tmpMat+tmpMat'
    end
    #Create random controlfields
    controlFields = randn(numControls, numTimeSteps)

    #Control frequencies
    controlFreqs = randn(numControls)

    return Hnat, controlHams, controlFields, controlFreqs
end


text:
function  facenumbering(p,t,elemtype,bndexpr,prdexpr)

dim = size(p,1);
ne = size(t,2);
face = getelemface(dim,elemtype);
nvf,nfe = size(face);

t2fl = reshape(t[face[:],:], nvf, nfe*ne);
#pf = reshape(p[:,t2fl[:]], dim, nvf, nfe*ne);
#pf = reshape(sum(pf,dims=2)/nvf,dim,nfe,ne);
pf = reshape(p[:,t2fl[:]], dim, nvf, nfe, ne);

# interior faces are zero
f = zeros(Int,nfe,ne);

# face-to-element connectivity
#f2t = mkf2t(t,elemtype,dim);
f2t,t2t = mkf2e(t,elemtype,dim);

# find elements on the domain boundary
ind = findall(f2t[3,:].==0);
for i = 1:length(ind) # for each element on the domain boundary
    e = f2t[1,ind[i]]; # element e
    l = f2t[2,ind[i]]; # local face index
    for k = 1:length(bndexpr) # for each boundary expression
        #a = bndexpr[k](pf[:,l,e]); # evaluate the boundary expression
        a1 = bndexpr[k](pf[:,1,l,e]);
        a2 = bndexpr[k](pf[:,2,l,e]);
        a3 = bndexpr[k](pf[:,nvf,l,e]);
        if a1[1] && a2[1] && a3[1] # check if element e belong to this boundary
            f[l,e] = k; # then set f(l,e) to k
            break;
        end
    end
end
# # regular boundary faces are positive
# for i = 1:length(bndexpr)
#     tm = bndexpr[i](pf);
#     f[tm .== 1] .= i;
# end

pf = reshape(sum(pf,dims=2)/nvf,dim,nfe,ne);
nprd = size(prdexpr,1);
if nprd>0
    f = f[:];
    pf = reshape(pf, dim, nfe*ne);
    tprd = copy(t);
    # periodic boundary faces are negative
    for i = 1:nprd
        i1 = findall(f .== prdexpr[i,1]);
        f[i1] = -f[i1];
        p1 = prdexpr[i,2](pf[:,i1[:]]);
        e1 = Int.(ceil.(Float64.(i1)./Float64(nfe)));      # 1st elements
        l1 = i1 .- (e1.-1).*nfe;   # 1st local faces

        i2 = findall(f .== prdexpr[i,3]);
        f[i2] = -f[i2];
        p2 = prdexpr[i,4](pf[:,i2[:]]);
        e2 = Int.(ceil.(Float64.(i2)./Float64(nfe)));      # 1st elements
        l2 = i2 .- (e2.-1).*nfe;   # 1st local faces

        # update t2t to connect periodic elements
        for j=1:length(e1)
            t2t[l2[j],e2[j]] = e1[j];
            t2t[l1[j],e1[j]] = e2[j];
        end

        p1 = reshape(p1, Int64(length(p1)/length(i1)), length(i1));
        p2 = reshape(p2, Int64(length(p2)/length(i2)), length(i2));
        in = xiny(p1,p2,0);

        v1 = t2fl[:,i1];
        v1 = sort(unique(v1[:]));
        p1 = prdexpr[i,2](p[:,v1]);

        v2 = t2fl[:,i2[in]];
        v2 = sort(unique(v2[:]));
        p2 = prdexpr[i,4](p[:,v2]);

        p1 = reshape(p1, Int64(length(p1)/length(v1)), length(v1));
        p2 = reshape(p2, Int64(length(p2)/length(v2)), length(v2));
        in = xiny(p1,p2,0);
        v2 = v2[in];

        for j=1:length(v1)
            tprd[tprd .== v2[j]] .= v1[j];
        end
    end
else
    #tprd = zeros(1,1);
    tprd = t;
end

f = reshape(f, nfe, ne);

return f, tprd, t2t

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim