text:
################################################################################
#
#          EllCrv/EllCrv.jl : Elliptic curves over general fields
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2016 Tommy Hofmann
# (C) 2016 Robin Ammon
# (C) 2016 Sofia Brenner
#
################################################################################

export EllCrv, EllCrvPt

export base_field, disc, division_polynomialE, EllipticCurve, infinity,
       isinfinite, isshort, ison_curve, j_invariant, Psi_polynomial,
       psi_poly_field, short_weierstrass_model, +, *

################################################################################
#
#  Types
#
################################################################################

mutable struct EllCrv{T}
  base_field::Ring
  short::Bool
  coeff::Array{T, 1}
  a_invars::Tuple{T, T, T, T, T}
  b_invars::Tuple{T, T, T, T}
  long_c::Array{T, 1}
  disc::T
  j::T

  torsion_points#::Array{EllCrvPt, 1}
  torsion_structure#Tuple{Array{Int, 1}, Array{EllCrvPt, 1}}

  function EllCrv{T}(coeffs::Array{T, 1}, check::Bool = true) where {T}
    if length(coeffs) == 2
      if check
        d = 4*coeffs[1]^3 + 27*coeffs[2]^2
        if d != 0
          E = new{T}()
          E.short = true
          # fixed on Nemo master
          E.coeff = [ deepcopy(z) for z in coeffs]
          E.base_field = parent(coeffs[1])
          #E.disc = d
        else
          error("discriminant is zero")
        end
      else
        E = new{T}()
        E.short = true
        E.coeff = [ deepcopy(z) for z in coeffs]
        E.base_field = parent(coeffs[1])
      end
    elseif length(coeffs) == 5 # coeffs = [a1, a2, a3, a4, a6]
      if check
        a1 = coeffs[1]
        a2 = coeffs[2]
        a3 = coeffs[3]
        a4 = coeffs[4]
        a6 = coeffs[5]

        b2 = a1^2 + 4*a2
        b4 = a1*a3 + 2*a4
        b6 = a3^2 + 4*a6
        b8 = a1^2*a6 - a1*a3*a4 + 4*a2*a6 + a2*a3^2 - a4^2

        d = (-b2^2*b8 - 8*b4^3 - 27*b6^2 + 9*b2*b4*b6)

        if d != 0
          E = new{T}()
          E.coeff = [ deepcopy(z) for z in coeffs]
          E.short = false
          E.b_invars = (b2, b4, b6, b8)
          E.a_invars = (a1, a2, a3, a4, a6)
          E.disc = d
          E.base_field = parent(coeffs[1])
        else
          error("Discriminant is zero")
        end
      else
        E = new{T}()
        E.short = false
        E.coeff = [ deepcopy(z) for z in coeffs]
        E.base_field = parent(coeffs[1])
      end
    else
      error("Length of coefficient array must be 2 or 5")
    end
    return E
  end
end

mutable struct EllCrvPt{T}
  coordx::T
  coordy::T
  isinfinite::Bool
  parent::EllCrv{T}

  function EllCrvPt{T}(E::EllCrv{T}, coords::Array{T, 1}, check::Bool = true) where {T}
    if check
      if ison_curve(E, coords)
        P = new{T}(coords[1], coords[2], false, E)
        return P
      else
        error("Point is not on the curve")
      end
    else
      P = new{T}(coords[1], coords[2], false, E)
      return P
    end
  end

  function EllCrvPt{T}(E::EllCrv{T}) where {T}
    z = new{T}()
    z.parent = E
    z.isinfinite = true
    return z
  end
end

################################################################################
#
#  Constructors for users
#
################################################################################

function EllipticCurve(x::Array{T, 1}, check::Bool = true) where T
  E = EllCrv{T}(x, check)
  return E
end

#  Implicit promotion in characterstic 0
function EllipticCurve(x::Array{Int, 1}, check::Bool = true)
  return EllipticCurve(fmpq[ FlintQQ(z) for z in x], check)
end

function EllipticCurve(x::Array{fmpz, 1}, check::Bool = true)
  return EllipticCurve(fmpq[ FlintQQ(z) for z in x], check)
end

function (E::EllCrv{T})(coords::Array{S, 1}, check::Bool = true) where {S, T}
  if length(coords) != 2
    error("Need two coordinates")
  end

  if S == T
    parent(coords[1]) != base_field(E) &&
        error("Objects must be defined over same field")
    return EllCrvPt{T}(E, coords, check)
  else
    return EllCrvPt{T}(E, map(base_field(E), coords), check)
  end
end

################################################################################
#
#  Field access
#
################################################################################

function base_field(E::EllCrv{T}) where T
  return E.base_field::parent_type(T)
end

function Base.deepcopy_internal(E::EllCrv, dict::IdDict)
    return EllipticCurve(E.coeff)
end

function parent(P::EllCrvPt)
  return P.parent
end

function isfinite(P::EllCrvPt)
  return !P.isinfinite
end

function isinfinite(P::EllCrvPt)
  return P.isinfinite
end

function isshort(E::EllCrv)
  return E.short
end

function a_invars(E::EllCrv)
  if isdefined(E, :a_invars)
    return [ deepcopy(z) for z in E.a_invars ]
  else
    t = (E.coeff[1], E.coeff[2], E.coeff[3], E.coeff[4], E.coeff[5])
    E.a_invars = t
    return t
  end
end

function b_invars(E::EllCrv)
  if isdefined(E, :long_b)
    # fixed on Nemo master
    return deepcopy(E.b_invars)
  else
    a1 = E.coeff[1]
    a2 = E.coeff[2]
    a3 = E.coeff[3]
    a4 = E.coeff[4]
    a6 = E.coeff[5]

    b2 = a1^2 + 4*a2
    b4 = a1*a3 + 2*a4
    b6 = a3^2 + 4*a6
    b8 = a1^2*a6 - a1*a3*a4 + 4*a2*a6 + a2*a3^2 - a4^2

    E.b_invars = (b2, b4, b6, b8)
    return (b2, b4, b6, b8)
  end
end

################################################################################
#
# Changing the model
#
################################################################################

@doc Markdown.doc"""
    short_weierstrass_model(E::EllCrv{fmpq}) -> (EE::EllCrv, function(EllCrvPt), function(EllCrvPt))

Transforms a curve given in long Weierstrass form over QQ to short Weierstrass form
returns short form and both transformations for points on the curve; first transformation from E (long form) to EE (short form), second transformation the other way round
"""
function short_weierstrass_model(E::EllCrv)
  return _short_weierstrass_model(E)
end

function _short_weierstrass_model(E::EllCrv{T}) where T
  R = base_field(E)
  p = characteristic(R)

  if (p == 2) || (p == 3)
      error("Converting to short form not possible in characteristic 2 and 3")
  end

  a1 = E.coeff[1]
  a2 = E.coeff[2]
  a3 = E.coeff[3]
  a4 = E.coeff[4]
  a6 = E.coeff[5]

  b2, b4, b6, b8 = b_invars(E)

#  b2 = a1^2 + 4*a2
#  b4 = a1*a3 + 2*a4
#  b6 = a3^2 + 4*a6

  c4 = (b2^2 - 24*b4)
  c6 = (-b2^3 + 36* b2*b4 - 216* b6)

  Anew = -divexact(c4, 48)
  Bnew = -divexact(c6, 864)

  EE = EllipticCurve([Anew, Bnew])::EllCrv{T}

  # we are hitting https://github.com/JuliaLang/julia/issues/15276

  _b2 = deepcopy(b2)
  _a1 = deepcopy(a1)
  _a3 = deepcopy(a3)

  # transformes a point on E (long form) to a point on EE (short form)
  trafo = function(P::EllCrvPt)

    if P.isinfinite
      return infinity(EE)
    end

    xnew = P.coordx + divexact(_b2, 12)
    ynew = P.coordy + divexact(_a1*P.coordx + _a3, 2)
    Q = EE([xnew, ynew])::EllCrvPt{T}
    return Q
  end

  # transformes a point on EE (short form) back to a point on E (long form)
  ruecktrafo = function(R::EllCrvPt)
    if R.isinfinite
        return infinity(E)
    end

    xnew = R.coordx - divexact(_b2, 12)
    ynew = R.coordy - divexact(_a1*xnew + _a3, 2)
    S = E([xnew, ynew])
    return S::EllCrvPt{T}
  end

  # type annotation necessary due to #15276
  return EE::EllCrv{T}, trafo, ruecktrafo
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, E::EllCrv)
  print(io, "Elliptic curve with equation\n")
  if E.short
    print(io, "y^2 = ")
    sum = Expr(:call, :+)
    push!(sum.args, Expr(:call, :^, :x, 3))
    c = E.coeff[1]
    if !iszero(c)
      if isone(c)
        push!(sum.args, Expr(:call, :*, :x))
      else
        push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), :x))
      end
    end
    c = E.coeff[2]
    if !iszero(c)
      push!(sum.args, AbstractAlgebra.expressify(c, context = io))
    end
    print(io, AbstractAlgebra.expr_to_string(AbstractAlgebra.canonicalize(sum)))
  else
    sum = Expr(:call, :+)
    push!(sum.args, Expr(:call, :^, :y, 2))
    c = E.coeff[1]
    if !iszero(c)
      if isone(c)
        push!(sum.args, Expr(:call, :*, :x, :y))
      else
        push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), :x, :y))
      end
    end
    c = E.coeff[3]
    if !iszero(c)
      if isone(c)
        push!(sum.args, Expr(:call, :*, :y))
      else
        push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), :y))
      end
    end
    print(io, AbstractAlgebra.expr_to_string(AbstractAlgebra.canonicalize(sum)))

    print(io, " = ")
    sum = Expr(:call, :+)
    push!(sum.args, Expr(:call, :^, :x, 3))

    c = E.coeff[2]
    if !iszero(c)
      if isone(c)
        push!(sum.args, Expr(:call, :*, Expr(:call, :^, :x, 2)))
      else
        push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), Expr(:call, :^, :x, 2)))
      end
    end

    c = E.coeff[4]
    if !iszero(c)
      if isone(c)
        push!(sum.args, Expr(:call, :*, :x))
      else
        push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), :x))
      end
    end

    c = E.coeff[5]
    if !iszero(c)
      push!(sum.args, AbstractAlgebra.expressify(c, context = io))
    end

    print(io, AbstractAlgebra.expr_to_string(AbstractAlgebra.canonicalize(sum)))
  end
end

function show(io::IO, P::EllCrvPt)
    if P.isinfinite
        print(io, "Point at infinity of $(P.parent)")
    else
        print(io, "Point $(P.coordx),$(P.coordy) of $(P.parent)")
    end
end

################################################################################
#
#  Point at infinity
#
################################################################################

@doc Markdown.doc"""
    infinity(E::EllCrv) -> EllCrvPt

Creates the point at infinity.
"""
function infinity(E::EllCrv{T}) where T
  infi = EllCrvPt{T}(E)
  return infi
end

################################################################################
#
#  Test for inclusion
#
################################################################################

@doc Markdown.doc"""
    ison_curve(E::EllCrv{T}, coords::Array{T, 1}) -> Bool

Returns true if `coords` defines a point  on E and false otherwise. The array
`coords` must have length 2.
"""
function ison_curve(E::EllCrv{T}, coords::Array{T, 1}) where T
  length(coords) != 2 && error("Array must be of length 2")

  x = coords[1]
  y = coords[2]

  if E.short == true
    if y^2 == x^3 + (E.coeff[1])*x + (E.coeff[2])
      return true
    else
      return false
    end
  else
    if (y^2 + (E.coeff[1])*x*y + (E.coeff[3])*y ==
            x^3 + (E.coeff[2])*x^2+(E.coeff[4])*x + (E.coeff[5]))
      return true
    else
      return false
    end
  end
end

################################################################################
#
#  Discriminant
#
################################################################################

@doc Markdown.doc"""
    disc(E::EllCrv{T}) -> T

Computes the discriminant of $E$.
"""
function disc(E::EllCrv{T}) where T
  if isdefined(E, :disc)
    return E.disc
  end

  if E.short == true
    # fall back to the formula for the long form
    # this should be done in a more clever way
    R = base_field(E)
    F = EllipticCurve([R(0), R(0), R(0), E.coeff[1], E.coeff[2]])
    d = disc(F)
    E.disc = d
    return d::T
  else
    a1 = E.coeff[1]
    a2 = E.coeff[2]
    a3 = E.coeff[3]
    a4 = E.coeff[4]
    a6 = E.coeff[5]

    (b2, b4, b6, b8) = b_invars(E)

    d = -b2^2*b8 - 8*b4^3 - 27*b6^2 + 9*b2*b4*b6
    E.disc = d
    return d::T
  end
end

################################################################################
#
#  j-invariant
#
################################################################################

# p. 46 Washington, p. 72 Cohen
@doc Markdown.doc"""
    j(E::EllCrv{T}) -> T
Computes the j-invariant of $E$.
"""
function j_invariant(E::EllCrv{T}) where T
  if isdefined(E, :j)
    return E.j
  end

  if E.short == true

    R = base_field(E)
    F = EllipticCurve([R(0), R(0), R(0), E.coeff[1], E.coeff[2]])
    j = j_invariant(F)
    E.j = j
    return j::T
  else
    a1 = E.coeff[1]
    a2 = E.coeff[2]
    a3 = E.coeff[3]
    a4 = E.coeff[4]
    a6 = E.coeff[5]

    (b2, b4, b6, b8) = b_invars(E)
    #b2 = a1^2 + 4*a2
    #b4 = a1*a3 + 2*a4
    #b6 = a3^2 + 4*a6
    #b8 = a1^2*a6 - a1*a3*a4 + 4*a2*a6 + a2*a3^2 - a4^2

    c4 = b2^2 - 24*b4

    j = divexact(c4^3, disc(E))
    E.j = j
    return j::T
  end
end

################################################################################
#
#  Addition of Points
#
################################################################################

# washington p. 14, cohen p. 270
@doc Markdown.doc"""
    +(P::EllCrvPt, Q::EllCrvPt) -> EllCrvPt
Adds two points on an elliptic curve.
does not work in characteristic 2
"""
function +(P::EllCrvPt{T}, Q::EllCrvPt{T}) where T
  parent(P) != parent(Q) && error("Points must live on the same curve")

  characteristic(base_field(parent(P))) == 2 &&
      error("Not yet implemented in characteristic 2")

  # Is P = infinity or Q = infinity?
  if P.isinfinite
      return Q
  elseif Q.isinfinite
      return P
  end

  E = P.parent

  # Distinguish between long and short form
  if E.short == true
    if P.coordx != Q.coordx
        m = divexact(Q.coordy - P.coordy, Q.coordx - P.coordx)
        x = m^2 - P.coordx - Q.coordx
        y = m * (P.coordx - x) - P.coordy
    elseif P.coordy != Q.coordy
        return infinity(E)
    elseif P.coordy != 0
        m = divexact(3*(P.coordx)^2 + (E.coeff[1]), 2*P.coordy)
        x = m^2 - 2*P.coordx
        y = m* (P.coordx - x) - P.coordy
    else
        return infinity(E)
    end

    Erg = E([x, y], false)

  else
    a1 = E.coeff[1]
    a2 = E.coeff[2]
    a3 = E.coeff[3]
    a4 = E.coeff[4]
    a6 = E.coeff[5]

    # Use [Cohen, p. 270]
    if P.coordx == Q.coordx
      if Q.coordy == -a1*P.coordx - a3 - P.coordy # then P = -Q
        return infinity(E)
      elseif P.coordy == Q.coordy # then P = Q
        m = divexact(3*((P.coordx)^2) + 2*a2*P.coordx + a4 - a1*P.coordy, 2*P.coordy + a1*P.coordx + a3)
        x = -P.coordx - Q.coordx - a2 + a1*m + m^2
        y = -P.coordy - m*(x - P.coordx) - a1*x - a3
      else # then P != +-Q
        m = divexact(Q.coordy - P.coordy, Q.coordx - P.coordx)
        x = -P.coordx - Q.coordx - a2 + a1*m + m^2
        y = -P.coordy - m*(x - P.coordx) - a1*x - a3
      end
    else # now P != +-Q
      m = divexact(Q.coordy - P.coordy, Q.coordx - P.coordx)
      x = -P.coordx - Q.coordx - a2 + a1*m + m^2
      y = -P.coordy - m*(x - P.coordx) - a1*x - a3
    end

    Erg = E([x, y], false)

  end
  return Erg
end

################################################################################
#
#  Inverse
#
################################################################################

@doc Markdown.doc"""
    -(P::EllCrvPt) -> EllCrvPt

Computes the inverse of the point $P$ on an elliptic curve.
"""
function -(P::EllCrvPt)
  E = P.parent

  if !isfinite(P)
    return infinity(E)
  end

  if E.short == true
    Q = E([P.coordx, -P.coordy], false)
  else
    Q = E([P.coordx, -E.coeff[1]*P.coordx - E.coeff[3] - P.coordy], false)
  end

  return Q
end

@doc Markdown.doc"""
    ==(P::EllCrvPt, Q::EllCrvPt) -> Bool

Returns true if $P$ and $Q$ are equal and live over the same elliptic curve
$E$.
"""
function ==(P::EllCrvPt{T}, Q::EllCrvPt{T}) where T
  # both are infinite
  if P.isinfinite && Q.isinfinite
    return true
  end

  # one of them is infinite
  if xor(P.isinfinite, Q.isinfinite)
    return false
  end

  # Otherwise, compare coordinates
  if P.coordx == Q.coordx && P.coordy == Q.coordy
    return true
  else
    return false
  end
end

################################################################################
#
#  Scalar multiplication
#
################################################################################

# algorithm 'integer times a point', [Washington, p. 18]
@doc Markdown.doc"""
    *(n::Int, P::EllCrvPt) -> EllCrvPt

Computes the point $nP$.
"""
function *(n::Int, P::EllCrvPt)
  B = infinity(P.parent)
  C = P

  if n >= 0
      a = n
  else
      a = -n
  end

  while a != 0
    if mod(a,2) == 0
      a = div(a,2)
      C = C + C
    else
      a = a - 1
      B = B + C
    end
  end

  if n < 0
    B = -B
  end

  return B
end

################################################################################
#
#  Division polynomials
#
################################################################################
# TODO:
# - Adjust to Sage/Magma
# - Implement for long form, see Mazur, Tate, "p-adic sigma function"

# computes the n-th division polynomial psi_n in ZZ[x,y] for a given elliptic curve E over ZZ
function division_polynomialE(E::EllCrv, n::Int, x = nothing, y = nothing)

  A = numerator(E.coeff[1])
  B = numerator(E.coeff[2])

  if x === nothing
    Zx, _x = PolynomialRing(FlintZZ, "x")
    Zxy, y = PolynomialRing(Zx, "y")
    x = Zxy(_x)
  else
    Zxy = parent(x)
  end

  if n == 1
    return one(Zxy)
  elseif n == 2
    return 2*y
  elseif n == 3
    return 3*x^4 + 6*(A)*x^2 + 12*(B)*x - (A)^2
  elseif n == 4
    return 4*y*(x^6 + 5*(A)*x^4 + 20*(B)*x^3 - 5*(A)^2*x^2 - 4*(A)*(B)*x - 8*(B)^2 - (A)^3)
  elseif mod(n,2) == 0
    m = div(n,2)
    return divexact( (division_polynomialE(E,m,x,y))*(division_polynomialE(E,m+2, x, y)*division_polynomialE(E,m-1,x, y)^2 - division_polynomialE(E,m-2,x, y)*division_polynomialE(E,m+1,x,y)^2), 2*y)
  else m = div(n-1,2)
    return division_polynomialE(E,m+2,x,y)*division_polynomialE(E,m,x,y)^3 - division_polynomialE(E,m-1,x,y)*division_polynomialE(E,m+1,x,y)^3
  end
end

function replace_all_squares(f, g)
    # assumes that f is in Z[x,y^2] and g in Z[x]. Replaces y^2 with g.
    # the result will be in Z[x]
    z = zero(parent(g)) # this is the zero in Z[x]
    d = div(degree(f), 2) # degree of f in y^2 should be even
    for i in 0:d
        z = z + coeff(f, 2*i)*g^i
    end
    return z
end

function replace_all_squares_modulo(f, g, F)
    # assumes that f is in Z[x,y^2] and g in Z[x]. Replaces y^2 with g.
    # the result will be in Z[x]
    z = zero(parent(g)) # this is the zero in Z[x]
    d = div(degree(f), 2) # degree of f in y^2 should be even
    for i in 0:d
        z = z + coeff(f, 2*i)*powmod(g, i, F)
    end
    return z
end

# here is an example: f = _Zxy(_x)*y^2 + _Zxy(_x + 1)*y^4; EC.replace_all_squares(f, x^3)
# this should give _x^7+_x^6+_x^4

# computes the n-th Psi-polynomial Psi_n in ZZ[_x]
function Psi_polynomial(E::EllCrv, n::Int)

    if n < 2
        error("Psi-polynomial not defined")
    end

    Zx, _x = PolynomialRing(FlintZZ, "x")
    Zxy, y = PolynomialRing(Zx, "y")
    x = Zxy(_x)

    # g = y^2
    g = _x^3 + (numerator(E.coeff[1]))*_x + numerator(E.coeff[2])

    h = division_polynomialE(E, n, x, y)
    # make h into an element of ZZ[x]

    # in the even case, first divide by 2y and then replace y^2
    if mod(n,2) == 0
        f = divexact(h,2*y)
        f = replace_all_squares(f,g)
    else
        f = replace_all_squares(h,g)
    end

    return f
end

# Division polynomials in general for an elliptic curve over an arbitrary field

# standard divison polynomial Psi (as needed in Schoof's algorithm)
function psi_poly_field(E::EllCrv, n::Int, x, y)

    R = base_field(E)
    A = E.coeff[1]
    B = E.coeff[2]

    if n == -1
        return -y^0
    elseif n == 0
        return zero(parent(y))
    elseif n == 1
        return y^0
    elseif n == 2
        return 2*y
    elseif n == 3
        return (3*x^4 + 6*(A)*x^2 + 12*(B)*x - (A)^2)*y^0
    elseif n == 4
        return 4*y*(x^6 + 5*(A)*x^4 + 20*(B)*x^3 - 5*(A)^2*x^2 - 4*(A)*(B)*x - 8*(B)^2 - (A)^3)
    elseif mod(n,2) == 0
        m = div(n,2)
        return divexact( (psi_poly_field(E,m,x,y))*(psi_poly_field(E,m+2,x,y)*psi_poly_field(E,m-1,x,y)^2 - psi_poly_field(E,m-2,x,y)*psi_poly_field(E,m+1,x,y)^2), 2*y)
    else m = div(n-1,2)
        return psi_poly_field(E,m+2,x,y)*psi_poly_field(E,m,x,y)^3 - psi_poly_field(E,m-1,x,y)*psi_poly_field(E,m+1,x,y)^3
    end
end


text:
using BinDeps
using Compat

@BinDeps.setup
libnames = ["libCLBlast", "libclblast", "clblast"]
libCLBlast = library_dependency("libCLBlast", aliases = libnames)
version = "1.4.1"

if Compat.Sys.iswindows()
    if Sys.ARCH == :x86_64
        uri = URI("https://github.com/CNugteren/CLBlast/releases/download/" *
                  version * "/CLBlast-" * version * "-Windows-x64.zip")
        basedir = @__DIR__
        provides(
            Binaries, uri,
            libCLBlast, unpacked_dir = ".",
            installed_libpath = joinpath(basedir, "libCLBlast", "lib"), os = :Windows
        )
    else
        error("Only 64 bit windows supported with automatic build.")
    end
end

if Compat.Sys.islinux()
    #=if Sys.ARCH == :x86_64
        name, ext = splitext(splitext(basename(baseurl * "Linux-x64.tar.gz"))[1])
        uri = URI(baseurl * "Linux-x64.tar.gz")
        basedir = joinpath(@__DIR__, name)
        provides(
            Binaries, uri,
            libCLBlast, unpacked_dir = basedir,
            installed_libpath = joinpath(basedir, "lib"), os = :Linux
        )
    else
        error("Only 64 bit linux supported with automatic build.")
    end=#
    provides(Sources, URI("https://github.com/CNugteren/CLBlast/archive/" * version * ".tar.gz"),
             libCLBlast, unpacked_dir="CLBlast-" * version)

    builddir = joinpath(@__DIR__, "src", "CLBlast-" * version, "build")
    libpath = joinpath(builddir, "libclblast.so")
    provides(BuildProcess,
        (@build_steps begin
            GetSources(libCLBlast)
            CreateDirectory(builddir)
            FileRule(libpath, @build_steps begin
                ChangeDirectory(builddir)
                `cmake -DSAMPLES=OFF -DTESTS=OFF -DTUNERS=OFF -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ ..`
                `make -j 4`
            end)
        end),
        libCLBlast, installed_libpath=libpath, os=:Linux)
end

if Compat.Sys.isapple()
    using Homebrew
    provides(Homebrew.HB, "homebrew/core/clblast", libCLBlast, os = :Darwin)
end

if Compat.Sys.islinux()
    # BinDeps.jl seems to be broken, cf. https://github.com/JuliaLang/BinDeps.jl/issues/172
    wd = pwd()
    sourcedir = joinpath(@__DIR__, "CLBlast-" * version)
    builddir = joinpath(sourcedir, "build")
    libpath = joinpath(builddir, "libclblast.so")
    if !isdir(sourcedir)
        url = "https://github.com/CNugteren/CLBlast/archive/" * version * ".tar.gz"
        download_finished = false
        try
            run(pipeline(`wget -q -O - $url`, `tar xzf -`))
            global download_finished = true
        catch e
            println(e)
        end
        if download_finished == false
            try
                run(pipeline(`curl -L $url`, `tar xzf -`))
                global download_finished = true
            catch e
                println(e)
            end
        end
    end
    isdir(builddir) || mkdir(builddir)
    cd(builddir)
    run(`cmake -DSAMPLES=OFF -DTESTS=OFF -DTUNERS=OFF -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ ..`)
    run(`make -j 4`)
    cd(wd)
    open(joinpath(@__DIR__, "deps.jl"), "w") do file
        write(file,
"""
if VERSION >= v"0.7.0-DEV.3382"
    using Libdl
end
# Macro to load a library
macro checked_lib(libname, path)
    if Libdl.dlopen_e(path) == C_NULL
        error("Unable to load \n\n\$libname (\$path)\n\nPlease ",
              "re-run Pkg.build(CLBlast), and restart Julia.")
    end
    quote
        const \$(esc(libname)) = \$path
    end
end

# Load dependencies
@checked_lib libCLBlast "$libpath"
""")
    end
else
    @BinDeps.install Dict(:libCLBlast => :libCLBlast)
end

ext: jl


text:
"""
    push_integral_or_branch!(step_obj, cnode)

Add integral or branch node to step_obj
"""
function push_integral_or_branch!(step_obj, cnode)
    # check if all int vars are int
    if cnode.state == :Integral
        push!(step_obj.integral, cnode)
    else
        push!(step_obj.branch, cnode)
    end
end

"""
    set_cnode_state!(cnode, m, step_obj, disc2var_idx)

Set the state of the current node to :Integral or :Branch
"""
function set_cnode_state!(cnode, m, step_obj, disc2var_idx)
    # check if all int vars are int
    if are_type_correct(cnode.solution, m.var_type, disc2var_idx, m.options.atol)
        if only_almost_solved(cnode.relaxation_state) && m.options.allow_almost_solved_integral
            @warn "Integral leaf node only almost locally solved. Disable with `allow_almost_solved_integral=false`"
        end
        cnode.state = :Integral
    else
        cnode.state = :Branch
    end
end

"""
    new_integral!(tree, node)

Update the incumbent and add obj constr if in options
Add to solutions if list_of_solutions in options
"""
function new_integral!(tree, node)
    # node.best_bound is the objective for integral values
    tree.nsolutions += 1
    if tree.options.list_of_solutions
        push!(tree.m.solutions, Juniper.SolutionObj(node.solution, node.best_bound))
    end
    if update_incumbent!(tree,node) # returns if new
        if tree.options.incumbent_constr
            if tree.options.processors > 1
                np = nprocs()  # determine the number of processes available
                if tree.options.processors+1 < np
                    np = tree.options.processors+1
                end
                for p=2:np
                    sendto(p, is_newincumbent=true)
                end
            end
            add_incumbent_constr(tree.m, tree.incumbent)
            tree.m.ncuts += 1
        end
    end
end

"""
    push_to_branch_list!(tree, node)

Push a node to the list of branch nodes if better than incumbent
or if all solutions are wanted
"""
function push_to_branch_list!(tree, node)
    incu = isdefined(tree,:incumbent) ? tree.incumbent : false
    if tree.options.all_solutions || !isdefined(tree,:incumbent) || tree.obj_fac*node.best_bound >= tree.obj_fac*incu.objval
        push!(tree.branch_nodes, node)
    end
end

"""
    upd_integral_branch!(tree, step_obj)

Update the list of integral and branch nodes using the new step_obj
Return true if break
"""
function upd_integral_branch!(tree, step_obj)
    for integral_node in step_obj.integral
        new_integral!(tree,integral_node)
        if isbreak_new_incumbent_limits(tree)
            return true
        end
    end

    for branch_node in step_obj.branch
        push_to_branch_list!(tree, branch_node)
    end
    return false
end


text:
module VulkanCore

include("LibVulkan.jl")
using .LibVulkan

const api = LibVulkan
include(joinpath(@__DIR__, "..", "docs", "vulkan_doc.jl"))

const vk = api
export vk

end # module


text:
# Sample script for plotting fieldlines with handpicked seeds.
#
# Hongyang Zhou, hyzhou@umich.edu

# FieldTracer is not a registered pkg yet
# using Pkg
# Pkg.add(url="https://github.com/henry2004y/FieldTracer.jl.git", rev="master")

using Vlasiator, PyPlot, FieldTracer

file = "bulk.0000999.vlsv"
meta = load(file)
coordmin, coordmax, ncells = meta.coordmin, meta.coordmax, meta.ncells
pcolormesh(meta, "proton/vg_rho", colorscale=Linear)

Re = Vlasiator.Re
# regular Cartesian mesh
gridx = range(coordmin[1], coordmax[1], length=ncells[1])
gridy = range(coordmin[3], coordmax[3], length=ncells[3])

ns1 = 14 # number of seeding points on the dayside
ns2 = 10 # number of seeding points on the nightside
ns4 = 4  # number of seeding points in the polar region
seeds = Matrix{Float64}(undef, 2, ns1+2*ns2+2*ns4)
for i in 1:ns1
   seeds[1,i] = coordmax[1] / ns1 * (i - 1)
   seeds[2,i] = 0.0
end
for i in 1:ns2
   seeds[1,ns1+i] = -10Re
   seeds[2,ns1+i] = coordmin[3] + (coordmax[3] - coordmin[3]) / ns2 * (i - 1)
   seeds[1,ns1+ns2+i] = -30Re
   seeds[2,ns1+ns2+i] = coordmin[3] + (coordmax[3] - coordmin[3]) / ns2 * (i - 1)
end
for i in 1:ns4
   seeds[1,ns1+2*ns2+i] = -20Re + 10Re*(i - 1)
   seeds[2,ns1+2*ns2+i] = 30Re
   seeds[1,ns1+2*ns2+ns4+i] = -20Re + 10Re*(i - 1)
   seeds[2,ns1+2*ns2+ns4+i] = -30Re
end

b = meta["vg_b_vol"]
bx = reshape(b[1,:], ncells[1], ncells[3])
bz = reshape(b[3,:], ncells[1], ncells[3])

for i = axes(seeds,2)
   startx, starty = seeds[:,i]
   x1, y1 = trace2d(bx, bz, startx, starty, gridx, gridy;
      ds=0.5, maxstep=3000, gridType="ndgrid")
   x1 ./= Re
   y1 ./= Re
   if length(x1) < 5; continue; end
   line = plot(x1, y1, color="w")
   add_arrow(line[1])
end


text:
const PERFORM_CHECKS = false

using MAT
using DataFrames
using CSV

using ForwardDiff: jacobian
using SparseArrays: sparse
using LinearAlgebra: eigvals

function map_complex_2D(array)
    n = length(array)
    out = zeros(2n)
    for i = 1:n
        out[2i-1] = real(array[i])
        out[2i] = imag(array[i])
    end
    return out
end

function map_2D_complex(array)
    return array[1:2:end] .+ 1im .* array[2:2:end]
end

function stable_power_flow(pg, ic_guess::AbstractArray; sparse_jac=false, project=false, eig=true, p=nothing, t0=0.0)
    pgr = rhs(pg)
    M = Array(pgr.mass_matrix)
    n = systemsize(pg)

    f!(dx, x) = pgr(dx, x, p, t0)
    f(x) = (dx = similar(x); f!(dx, x); dx)
    j!(J, x) = (J .= jacobian(f, x))
    j(x) = (dx = similar(x); jacobian(f!, dx, x))

    if project
        mpm = M * pinv(M)
        proj_rhs!(dx, x) = (f!(dx, x); dx .-= mpm * dx)
        proj_guess = nlsolve(proj_rhs!, ic_guess; xtol=1E-10)
    end

    if sparse_jac
        # we can specify the Jacobian to be sparse
        dx0 = similar(ic_guess)
        F0 = similar(ic_guess)
        J0 = sparse(zeros(n, n))
        df = OnceDifferentiable(f!, j!, dx0, F0, J0)
        if project
            fixed_point = nlsolve(df, proj_guess.zero; xtol = 1E-10)
        else
            fixed_point = nlsolve(df, ic_guess; xtol = 1E-10)
        end
    else
        if project
            fixed_point = nlsolve(f!, j!, proj_guess.zero; xtol = 1E-10)
        else
            fixed_point = nlsolve(f!, j!, ic_guess; xtol = 1E-10)
        end
    end

    if eig
        λ = eigvals(Array(M), j(fixed_point.zero)) .|> real |> extrema
        println("Jacobian spectrum \nmin : ", first(λ), "\nmax : ",last(λ), "\nstable : ", last(λ)==0)
    end

    return State(pg, fixed_point.zero)
end

stable_power_flow(
    pg,
    ic_guess::State;
    sparse_jac = false,
    project = false,
    eig = true,
    p = nothing,
    t0 = 0.0,
) = stable_power_flow(
    pg,
    ic_guess.vec;
    sparse_jac = sparse_jac,
    project = project,
    eig = eig,
    p = p,
    t0 = t0,
)

function pf_sol(pg, ic_guess, p)
    #### solve power flow
    root_rhs = (dx, x) -> rhs(pg)(dx, x, p, 0.)
    pf = nlsolve(root_rhs, ic_guess, autodiff = :forward, method = :newton, xtol = 1E-10).zero #
    return State(pg, pf)
end

function pf_sol(busses_static, lines; matfile="$dir/LFCigre.mat", outfile="$dir/power_flow_sol.csv")

    #### read from Matlab
    matlab_data = matread(matfile)
    uK = Array(matlab_data["uK"])[:]
    uK[1] *= 1e3/base_voltage_HV
    uK[2:end] .*= 1e3/base_voltage
    sK = Array(matlab_data["sK"])[:]

    #### solve power flow

    pg = PowerGrid(busses_static, lines)
    rpg! = rhs(pg)

    ic_guess = map_complex_2D(uK)
    root_rhs = RootRhs_ic(rhs(pg))
    pf = nlsolve(root_rhs, ic_guess, autodiff = :forward,  method = :newton, xtol = 1E-10).zero #

    line_currents = [PowerDynamics.get_current(State(pg, pf), k) for k = 1:length(busses_static)]
    sK_julia = map_2D_complex(pf) .* conj.(line_currents)

    power_flow_sol = DataFrame(
        V = abs.(map_2D_complex(pf)),
        V_matlab = abs.(uK),
        A = rad2deg.(angle.(map_2D_complex(pf))),
        A_matlab = rad2deg.(angle.(uK)),
        s_real = real.(sK_julia),
        s_real_matlab = real.(sK),
        s_imag = imag.(sK_julia),
        s_imag_matlab = imag.(sK),
        real = pf[1:2:end],
        imag = pf[2:2:end],)

    if outfile != nothing
        CSV.write(outfile, power_flow_sol)
    end

    return power_flow_sol
end

if PERFORM_CHECKS
    using GraphPlot
    using Plots
    using BenchmarkTools

    #### read from Matlab
    Q_marcel = matlab_data["LeistungsflussQuelle"]["Q0las"]
    P_marcel = matlab_data["LeistungsflussQuelle"]["P0las"]

    #### check output

    pf = map_complex_2D(complex.(power_flow_sol.real, power_flow_sol.imag))
    niceprint(pf)

    du = similar(pf)
    rpg!(du, ic_guess, 0, 0)
    println(sum(du))

    #### plot result

    # nodesize: Spannungsabfall [%]
    val = 100abs.(maximum(power_flow_sol.V[2:end]) .- power_flow_sol.V) ./ 20
    val[1] = val[2]
    gplot(pg.graph, nodesize = val, nodelabel = v_angle(pf))

    #### compare with Matlab

    @assert any(Q_marcel[2:end] .- LoadQ .== 0)
    @assert any(P_marcel[2:end] .- LoadP .== 0)

    # Trafo
    # need to multiply with .* [-1 -1; 1 1] to take the sign convention in PowerDynamics into account
    matlab_data["Y"]["TTVsoll"] .- PiModel(lines[1]) .* [-1 -1; 1 1] |> sum |> println

    # line admittances Vierpoldarstellung
    # need to multiply with .* [-1 -1; 1 1] to take the sign convention in PowerDynamics into account
    for (i, l) in enumerate(lines[2:end])
        PiModel(l) .* [-1 -1; 1 1] .- matlab_data["Y"]["LLsoll"][2i-1:2i, 2i-1:2i] |> sum |> println
    end

    @show abs.(map_2D_complex(pf)) ./ abs.(uK)
    @show angle.(map_2D_complex(pf)) - angle.(uK) .|> rad2deg

    currents = [PowerDynamics.get_current(State(pg, pf), k) for k = 1:length(busses)]

    @show sK_julia = map_2D_complex(pf) .* conj.(currents)

    @show real.(sK_julia) ./ real.(sK)
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
############################################################################################
"""
$(SIGNATURES)
Contains information about log-likelihood, expected sample size and proposal trajectory.

# Examples
```julia
```

"""
struct ParticleFilterDiagnostics{P,G} <: AbstractDiagnostics
    "Diagnostics used for all Baytes kernels"
    base::BaytesCore.BaseDiagnostics{P}
    "Incremental log objective at current iteration."
    ℓincrement::Float64
    "Number of particles used in proposal steps."
    Nparticles::Int64
    "Average number of resampling steps."
    resampled::Float64
    "Generated quantities specified for objective"
    generated::G
    function ParticleFilterDiagnostics(
        base::BaytesCore.BaseDiagnostics{P},
        ℓincrement::Float64,
        Nparticles::Int64,
        resampled::Float64,
        generated::G,
    ) where {P, G}
        return new{P, G}(base, ℓincrement, Nparticles, resampled, generated)
    end
end

############################################################################################
"""
$(SIGNATURES)
Show relevant diagnostic results.

# Examples
```julia
```

"""
function generate_showvalues(diagnostics::D) where {D<:ParticleFilterDiagnostics}
    return function showvalues()
        return (:pf, "diagnostics"),
        (:loglik_estimate, diagnostics.base.ℓobjective),
        (:Temperature, diagnostics.base.temperature),
        (:Nparticles, diagnostics.Nparticles),
        (:resampled, diagnostics.resampled)
    end
end

############################################################################################
#export
export ParticleFilterDiagnostics, generate_showvalues


text:
using Folds
using FoldsPreview
using UnicodePlots

function online_extrema(;
    N = 100_000_000_000,
    interval = 10,
    ex = ThreadedEx(),
)
    vmin_trace = Float64[]
    vmax_trace = Float64[]
    t_trace = Float64[]
    t0 = time()
    function on_preview((vmin, vmax))
        push!(vmin_trace, -vmin)
        push!(vmax_trace, vmax)
        push!(t_trace, time() - t0)
        rows, columns = displaysize(stdout)
        opts = (height = max(5, rows ÷ 2 - 5), width = max(5, columns - 40))
        xs = log10.(t_trace)
        ylim = (
            floor(min(vmin_trace[1], vmax_trace[1]); digits = 3),
            ceil(max(-vmin, vmax); digits = 3),
        )
        plt1 = lineplot(xs, vmax_trace; ylabel = "extrema", ylim = ylim, opts...)
        lineplot!(plt1, xs, vmin_trace)
        plt2 = lineplot(
            xs,
            log10.(abs.(vmax_trace .- vmin_trace));
            ylabel = "log10(diff)",
            opts...,
        )
        display(plt1)
        display(plt2)
    end
    with_preview(on_preview, ex, interval) do previewer
        Folds.extrema(previewer(1:N), ex) do _
            randn()
        end
    end
end


text:
# This file is a part of TextClassification.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using MLDataUtils, LinearAlgebra
import TextSearch: vectorize, vectorize_corpus
import StatsBase: predict
import Base: hash, isequal
export filtered_power_set, predict, predict_corpus, vectorize, vectorize_corpus,
        MicroTC, AngleDistance, CosineDistance, NormalizedAngleDistance, NormalizedCosineDistance
import Base: hash, isequal
using SparseArrays

struct MicroTC{C_<:MicroTC_Config, CLS_<:Any, TextModel_<:TextModel, LabelType<:Any}
    config::C_
    cls::CLS_
    textmodel::TextModel_
    tok::Tokenizer
    levels::Vector{LabelType}
end

function Base.show(io::IO, model::MicroTC)
    print(io, "{MicroTC")
    show(io, model.config)
    print(io, " ")
    show(io, typeof(model.cls))
    print(io, " ")
    show(io, model.textmodel)
    print(io, " ")
    show(io, model.tok)
    print(io, " ")
    show(io, model.levels)
    print(io, "}")
end

Base.copy(c::MicroTC; config=c.config, cls=c.cls, textmodel=c.textmodel) = MicroTC(config, cls, textmodel)
Base.broadcastable(tc::MicroTC) = (tc,)

"""
    create(config, train_X, train_y) # config describes a text model
    create(config, train_X, train_y, tok, dim) # config describes a classifier

Creates a new object from a configuration and a train / test datasets.

- train_X is an array of BOW objects for text models
- train_X is an array of SVEC objects for classifiers
"""
create(config::KncConfig, train_X, train_y, dim) = Knc(config, train_X, train_y)
create(config::KncProtoConfig, train_X, train_y, dim) = KncProto(config, train_X, train_y)

"""
    MicroTC(
        config::MicroTC_Config,
        train_corpus::AbstractVector,
        train_y::CategoricalArray;
        tok=Tokenizer(config.textconfig),
        verbose=true)
    MicroTC(config::MicroTC_Config, textmodel::TextModel, train_X::AbstractVector{S}, train_y::CategoricalArray; verbose=true) where {S<:SVEC}

Creates a MicroTC model on the given dataset and configuration
"""
function MicroTC(
        config::MicroTC_Config,
        train_corpus::AbstractVector,
        train_y::CategoricalArray;
        tok=Tokenizer(config.textconfig, invmap=nothing),
        verbose=true)
    train_corpus_bow = compute_bow_corpus(tok, train_corpus)
    textmodel = create(config.textmodel, train_corpus_bow, train_y)
    X = SVEC[]
    mask = Bool[]
    for bow in train_corpus_bow
        x = vectorize(textmodel, bow)
        push!(mask, false)
        length(x) == 1 && haskey(x, 0) && continue
        push!(X, x)
        mask[end] = true
    end

    m = sum(mask)
    if m != length(mask)
        @info "WARNING considering $(m) of $(length(mask)) examples after vectorization using $(config)"
    end
    MicroTC(config, textmodel, X, train_y[mask], tok=Tokenizer(tok))  # tok is a copy with isconstruction=false
end

function MicroTC(
        config::MicroTC_Config,
        textmodel::TextModel,
        train_X::AbstractVector{S},
        train_y::CategoricalArray;
        tok=Tokenizer(config.textconfig, invmap=nothing),
        verbose=true) where {S<:SVEC}
    cls = create(config.cls, train_X, train_y, textmodel.m)
    MicroTC(config, cls, textmodel, tok, copy(levels(train_y)))
end

"""
    vectorize(tc::MicroTC, text;
        bow=BOW(),
        tok=Tokenizer(tc.config.textconfig, invmap=nothing),
        normalize=true)
    vectorize(tc::MicroTC, bow::BOW)

Creates a weighted vector using the model. The input `text` can be a string or an array of strings;
it also can be an already computed bag of words.

"""
function vectorize(
        tc::MicroTC,
        text;
        bow=BOW(),
        tok=tc.tok,
        normalize=true
    )::SVEC
    vectorize(tc.textmodel, compute_bow(tok, text, bow); normalize)
end

function vectorize(
        tc::MicroTC,
        bow::BOW;
        normalize=true
    )::SVEC
    vectorize(tc.textmodel, bow; normalize=normalize)
end

function vectorize_corpus(
        tc::MicroTC,
        corpus;
        bow=BOW(),
        tok=tok=tc.tok,
        normalize=true
    )
    V = Vector{SVEC}(undef, length(corpus))

    for i in eachindex(corpus)
        empty!(bow)
        V[i] = vectorize(tc, corpus[i]; bow, tok, normalize)
    end

    V
end

"""
    predict(tc::MicroTC, text)
    predict(tc::MicroTC, vec::SVEC)

Predicts the label of the given input
"""
predict(tc::MicroTC, text) = predict(tc.cls, vectorize(tc, text))
predict(tc::MicroTC, vec::SVEC) = tc.levels[predict(tc.cls, vec)]

function predict_corpus(tc::MicroTC, corpus;
    bow=BOW(),
    tok=tc.tok,
    normalize=true)
    V = Vector{eltype(tc.levels)}(undef, length(corpus))

    for i in eachindex(corpus)
        empty!(bow)
        empty!(tok)
        v = vectorize(tc, corpus[i]; bow, tok, normalize)
        V[i] = predict(tc, v)
    end

    V
end


text:
# phase-retrieval.jl
#
# Code for phase retrieval problems

# (A, b, x) = GeneratePhaseRetrievalData(m::Int64 = 200, n::Int64 = 20;
#                                        condition_number::Float64 = 1.0)
#
# Generates an m-by-n matrix A of measurement vectors a_i = A[i, :], where each
# a_i is uniform on the unit sphere. Chooses x ~ N(0, I), and then sets each
# observation b_i = (a_i' * x)^2.
function GeneratePhaseRetrievalData(m::Int64 = 200, n::Int64 = 20;
                                    condition_number::Float64 = 1.0)
  A = randn(m, n);
  A = A ./ sqrt.(repeat(sum(A.^2, dims=2), outer = (1, n)));
  if (condition_number > 1.0)
    # Make A have average norm 1 columns
    condition_vec = Vector(LinRange(1.0, condition_number, n));
    condition_vec = condition_vec / norm(condition_vec);
    A = A .* repeat(condition_vec', outer = (m, 1));
  end
  x = randn(n);
  b = (A * x).^2;
  return (A, b, x);
end

# PhaseRetrievalObj(A::Matrix{Float64}, b::Vector{Float64}, x::Vector{Float64})
#
# Computes the (average) absolute loss for the phase retrieval problem
#
#  p(b | lambda) = exp(-lambda) * lambda^b / b!
#
# For a data vector a, the parameter lambda = exp(a' * x), so we have log loss
# (ignoring the factorial term)
#
# l(x; (a, b)) = -b * dot(a, x) + exp(dot(a, x)).
function PhaseRetrievalObj(A::Matrix{Float64}, b::Vector{Float64},
                           x::Vector{Float64})
  mm = size(A, 1);
  objective = sum(abs.((A * x).^2 - b)) / mm;
  return objective;
end

# x = ProximalUpdatePhaseRetrieval(a::Vector{Float64}, b::Float64,
#                              x_init::Vector, stepsize::Float64; accuracy)
#
# Assuming b >= 0, sets x to be the minimizer of
#
#   |(a'*x)^2 - b| + ||x - x_init||^2 / (2 * alpha).
function ProximalUpdatePhaseRetrieval(a::Vector{Float64}, b::Float64,
                                      x_init::Vector{Float64},
                                      stepsize::Float64)
  # Replace this code to return the correct update.
  return x_init;
end

# x = TruncatedUpdatePhaseRetrieval(a, b, x_init, stepsize)
#
# Let F(x, a, b) = abs((a'*x)^2 - b) and
#
#  F_approx(x) = max{ F(x_init) + F'(x_init) * (x - x_init)), 0}
#
# Sets x to minimize
#
#  F_approx(x) + norm(x - x_init)^2 / (2 * stepsize).
#
# Returns the minimizing x.
function TruncatedUpdatePhaseRetrieval(a::Vector{Float64}, b::Float64,
                                       x_init::Vector{Float64},
                                       stepsize::Float64)
  # Replace this code to return the correct update.
  return x_init;
end

# x = SGUpdatePhaseRetrieval(a::Vector{Float64}, b::Float64,
#                     x_init::Vector{Float64}, stepsize::Float64)
#
# Computes the stochastic gradient update for PhaseRetrieval regression.
function SGUpdatePhaseRetrieval(a::Vector{Float64}, b::Float64,
                                x_init::Vector{Float64}, stepsize::Float64)
  # Replace this code to return the correct update.
  return x_init;
end

max_stars_repo_path: Machine Learning Summer School 2019 (London, UK)/tutorials/optimization/phase-retrieval.jl


text:
    function print_slice3(xst, xsp, yst, ysp, zst, zsp)
        print("[", xst, ":", xsp, ", ", yst, ":", ysp, ", ", zst, ":", zsp, "]")
    end

    function debug_slice3(xst, xsp, yst, ysp, zst, zsp)
        @debug "[", xst, ":", xsp, ", ", yst, ":", ysp, ", ", zst, ":", zsp, "]"
    end

    function print_array_stats(data3d)
        print(
            "size: ", size(data3d),", type: ", eltype(data3d),
            ", min, max: ", minimum(data3d), ", ", maximum(data3d), "\n"
        )
    end

    function debug_array_stats(data3d)
        @debug "size: ", size(data3d),", type: ", eltype(data3d),
            ", min, max: ", minimum(data3d), ", ", maximum(data3d), "\n"

    end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

include("RegistersTest.jl")
include("CircuitTest.jl")
include("InstructionsTest.jl")
include("CircuitSimplifyTest.jl")
include("CircuitUnitaryTest.jl")
include("GatesTest.jl")
include("OtherGatesTest.jl")


text:
function Binary_Search(x,key)     #declared function
mid = Int(floor((1 + length(x))/2))  #calculating the mid of array
           try
           if x[mid] == key
           println("Key found")        # if key is found
           elseif x[mid]< key           #comparision to see if right side of key is to iterated
           Binary_Search(x[mid+1:end],key) #recursive call
           else x[mid]>key            #comparision to see if left side of key is to iterated
           Binary_Search(x[1:mid-1],key)  #recursive call
           end
           catch
           println("Key: $key not found") # else key is not found
           end
end

x = [1,4,6,8,9,19,56]
key = 10
Binary_Search(x,key)

max_stars_repo_name: rajatenzyme/Coding-Journey-

max_issues_repo_name: rajatenzyme/Coding-Journey-

max_forks_repo_name: rajatenzyme/Coding-Journey-


text:
"""
Country fileformat

Create the link between the country code and the country name

Example:
#;COUNTRY;2;446;20181206;20190102;307;EAR_P
AG;Solomon Islands;N
AN;Nauru;N
AY;Papua New Guinea;N
BG;Greenland (Denmark);N
BI;Iceland;Y
"""
module Country

export read

struct CountryDc
    name::String
    member::Bool
    CountryDc(name, member) = new(name, member=="Y")
    # Country(name::AbstractString, member::Bool) = new(name, member)
end

function read(file)
    countries = Dict{String, CountryDc}()
    for line in eachline(file)
        line_elements = split(line, ';')
        if length(line_elements) == 3
            countries[line_elements[1]] = CountryDc(line_elements[2],
            line_elements[3])
        end
    end
    return countries
end

end # module


text:
"""
    ECMA_NUMBER_FORMAT

reference:[ECMA-376 standard](http://www.ecma-international.org/publications/standards/Ecma-376.htm)
           5th Edition Part1 - 18.8.30 numFmt (Number Format)

"""
const ECMA_NUMBER_FORMAT = Dict{UInt16, Tuple}(
   0 => (GeneralFormat,   "General"),
   1 => (Int,             "0"),
   2 => (Float64,         "0.00"),
   3 => (Float64,         "#,##0"),
   4 => (Float64,         "#,##0.00"),
   # not specified by ECMA-376 standard
   5 => (Float64,         "0.00"),
   6 => (Float64,         "0.00"),
   7 => (Float64,         "0.00"),
   8 => (Float64,         "0.00"),
   9 => (Float64,         "0%"),
  10 => (Float64,         "0.00%"),
  11 => (ScientificFormat,"0.00E+00"),
  12 => (Float64,         "# ?/?"),
  13 => (Float64,         "# ??/??"),
  14 => (ExcelDateTime,   "m/d/yy"),
  15 => (ExcelDateTime,   "d-mmm-yy"),
  16 => (ExcelDateTime,   "d-mmm"),
  17 => (ExcelDateTime,   "mmm-yy"),
  18 => (ExcelDateTime,   "h:mm AM/PM"),
  19 => (ExcelDateTime,   "h:mm:ss AM/PM"),
  20 => (ExcelDateTime,   "h:mm"),
  21 => (ExcelDateTime,   "h:mm:ss"),
  22 => (ExcelDateTime,   "m/d/yy h:mm"),
  27 => (ExcelDateTime,   "yyyy\"年\" mm\"月\" dd\"日\""),
  28 => (ExcelDateTime,   "mm-dd"),
  29 => (ExcelDateTime,   "mm-dd"),
  30 => (ExcelDateTime,   "mm-dd-yy"),
  31 => (ExcelDateTime,   "mm\"월\" dd\"일\""),
  32 => (ExcelDateTime,   "h\"시\" mm\"분\""),
  33 => (ExcelDateTime,   "h\"시\" mm\"분\" ss\"초\""),
  34 => (ExcelDateTime,   "yyyy-mm-dd"),
  35 => (ExcelDateTime,   "yyyy-mm-dd"),
  36 => (ExcelDateTime,   "mm\"月\" dd\"日\""),
  37 => (Float64,         "#,##0_);(#,##0)"),
  38 => (Float64,         "#,##0_);[Red](#,##0)"),
  39 => (Float64,         "#,##0.00_);(#,##0.00)"),
  40 => (Float64,         "#,##0.00_);[Red](#,##0.00)"),
  41 => (Float64,         "_(* #,##0_);_(* (#,##0);_(* \"-\"_);_(@_)"),
  42 => (Float64,         "_(u0024* #,##0_);_(u0024* (#,##0);_(u0024* \"-\"_);_(@_)"),
  43 => (Float64,         "_(* #,##0.00_);_(* (#,##0.00);_(* \"-\"??_);_(@_)"),
  44 => (Float64,         "_(u0024* #,##0.00_);_(u0024* (#,##0.00);_(u0024* \"-\"??_);_(@_)"),
  45 => (ExcelDateTime,   "mm:ss"),
  46 => (ExcelDateTime,   "[h]:mm:ss"),
  47 => (ExcelDateTime,   "mm:ss.0"),
  48 => (ScientificFormat,"##0.0E+0"),
  49 => (String,          "@"),
  50 => (ExcelDateTime,   "mm\"月\" dd\"日\""),
  51 => (ExcelDateTime,   "mm-dd"),
  52 => (ExcelDateTime,   "yyyy-mm-dd"),
  53 => (ExcelDateTime,   "yyyy-mm-dd"),
  54 => (ExcelDateTime,   "mm-dd"),
  55 => (ExcelDateTime,   "yyyy-mm-dd"),
  56 => (ExcelDateTime,   "yyyy-mm-dd"),
  57 => (ExcelDateTime,   "mm\"月\" dd\"日\""),
  58 => (ExcelDateTime,   "mm-dd"))

"""
    ECMA_CELLTYPE
e: treats all error type of 'ECMA_ERROR_VALUES' as 'missing'
n: excel treats all number as 'Float'

reference:[ECMA-376 standard](http://www.ecma-international.org/publications/standards/Ecma-376.htm)
           5th Edition Part1 - 18.18.11 ST_ST_CellType (Cell Type)
"""
const ECMA_CELLTYPE = Dict(
    "b"         => Bool,            # a Boolean
    "d"         => DateTime,        # a date in the ISO 8601 format
    "e"         => Missing,         # an error
    "inlineStr" => String,          # (inline) rich string
    "n"         => Float64,         # a Number
    "str"       => String,
    "s"         => SharedString)

"""
    ECMA_ERROR_VALUES
currently not being used in 'EzExcel'

reference:[ECMA-376 standard](http://www.ecma-international.org/publications/standards/Ecma-376.htm)
        5th Edition Part1 - L.2.16.8 Error values
"""
 const ECMA_ERROR_VALUES = Dict(
     "#NULL!" => "Intersection of two cell ranges is empty",
     "#DIV/0!"=> "Division by zero",
     "#VALUE!"=> "Wrong type of operand",
     "#REF!"  => "Illegal or deleted cell reference",
     "#NAME?" => "Wrong function or range name",
     "#NUM!"  => "Value range overflow",
     "#N/A"   => "Argument or function not available")

"""
    guess_datatype(fmtcode::AbstractString)
guess `DataType` from Excel format code
"""
function guess_datatype(fmtcode::AbstractString)
    T = String
    # seperating 'print format' and 'data format'
    print_fmt = search(fmtcode, r"\[(.*?)\]")
    data_fmt = first(print_fmt) == 0 ? fmtcode : fmtcode[last(print_fmt)+1:end]

    for reg in ((EzExcel.ExcelDateTime,    r"[ymdhs]"),
                (EzExcel.ScientificFormat, r"E\+|E\-"),
                (Float64,          r"0\.|#\.|0_|#_|\%|\?\/"),
                (Int,          r"0|#"))
        if ismatch(reg[2], data_fmt)
            T = reg[1]
            break
        end
    end
    # T == String && warn("failed to find 'DataType' for '$fmtcode'")
    return T
end


text:
#EasyPlot color tools & color schemes
#-------------------------------------------------------------------------------

#Alias in case underlying code changes:
getcolor(v::Symbol) = parse(Colorant, v)

#==Types
===============================================================================#
#ColorRef that can also reference another
#Int: Pick specific color from theme/ColorScheme
#Nothing: Pick appropriate color from theme/ColorScheme (Varies with sweep value)
const ColorRef = Union{Nothing, Colorant, Int}

#TODO: Possible to use structure from Colors module?
mutable struct ColorScheme
	colors::Vector{Colorant} #Colorant: might support alpha in the future
end

#==Color scheme definitions (not to be accessed directly)
===============================================================================#
const COLORSCHEME_DEFAULT = ColorScheme([getcolor(c) for c in
	[:blue, :yellow3, :red, :cyan, :green, :magenta, :orange, :violet, :grey50]
])

#==External interface for users
===============================================================================#
const COLORSCHEME = Dict{Symbol, ColorScheme}(
	:default => COLORSCHEME_DEFAULT,
)

getcolor(s::ColorScheme, ::Nothing) = colorant"black" #TODO: something better?

#Assumes value "v" starts at 1 (not 0):
function getcolor(s::ColorScheme, v::Integer)
	idx = mod(v-1, length(s.colors)) + 1
	return s.colors[idx]
end

#Default behaviour: ignore color scheme, and get color
#Alternative: Find closest color in scheme?
function getcolor(s::ColorScheme, v::Symbol)
	dfltcolor = getcolor(s, nothing)
	if :default == v
		return dfltcolor
	end

	result = dfltcolor
	try
		result = getcolor(v)
	catch
		@info("Color not supported: $v")
	end
	return result
end

#Default behaviour: ignore color scheme, and get color
#Alternative: Find closest color in scheme?
getcolor(s::ColorScheme, v::Colorant) = v

#Last line


text:
using JuMP
using GLPK
using Crayons
using OVERT
using Requires

"""
----------------------------------------------
main structure
----------------------------------------------
"""

"""
This structure is a mixed-integer-representation of Overt. The relu and max
operations are turned into mix integer programs following ideas from
MIPVerify paper https://arxiv.org/abs/1711.07356
"""
mutable struct OvertMIP
    overt_app::OverApproximation  # overt OverApproximation object
    model::JuMP.Model             # a JuMP model for MIP
    vars_dict::Dict{Symbol, JuMP.VariableRef}  # dictionary of Overt symbols and their associated variable in mip
end

gurobi_model(threads) = error("Gurobi not loaded")

DEFAULT_MODEL = "glpk"

function __init__()
    @require Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" begin
        DEFAULT_MODEL = "gurobi"

        gurobi_model(threads) = Model(optimizer_with_optimizer(Gurobi.Optimizer, "OutputFlag" => 0, "Threads" => threads))
    end
end

# default constructor
function OvertMIP(overt_app::OverApproximation; threads=0, model=DEFAULT_MODEL)
    if model == "glpk" || model == "GLPK"
        model = Model(GLPK.Optimizer)
        set_optimizer_attribute(model, "msg_lev", GLPK.MSG_OFF)
    elseif model == "gurobi" || model == "Gurobi"
        model = gurobi_model(threads)
    else
        error("Model not supported")
    end
    overt_mip_model = OvertMIP(overt_app,
                         model,
                         Dict{Symbol, JuMP.VariableRef}())
    overt_2_mip(overt_mip_model)
    return overt_mip_model
end

#TODO: secondary constructor with specified solver

Base.print(overt_mip_model::OvertMIP) = println(overt_mip_model.model)
Base.display(overt_mip_model::OvertMIP) = println(overt_mip_model.model)

"""
----------------------------------------------
some utility functions
----------------------------------------------
"""

"""
This function returns mip variable associated with a Symbol in Overt.
if the Symbol does not have an mip variable yet, one is created.
With a new variable created, the range of that variable will be added to the
mip model. The ranges are obtained from range attribute of OverApproximation object.
"""
function get_mip_var(var::Symbol, overt_mip_model::OvertMIP)
    if var in keys(overt_mip_model.vars_dict)
        mip_var = overt_mip_model.vars_dict[var]
    else
        l, u = overt_mip_model.overt_app.ranges[var]
        mip_var = @variable(overt_mip_model.model, lower_bound=l, upper_bound=u, base_name="$var")
        overt_mip_model.vars_dict[var] = mip_var
    end
    return mip_var
end

"""
This function creates an auxillary mip variable. These are the intermediate defined
    when relu or max are turned into mip.
"""
aux_var_counter = 1
function get_mip_aux_var(overt_mip_model::OvertMIP; binary::Bool = false)
    global aux_var_counter;
    base_name = "aux_var_$aux_var_counter"
    mip_aux_var = @variable(overt_mip_model.model, binary=binary, base_name=base_name)
    aux_var_counter += 1
    return mip_aux_var
end

"""
----------------------------------------------
main loop
----------------------------------------------
"""

"""
This function reads an OverApproximation object and translate all constraints
    in terms of a Mixed Integer Program.
    input: overt_mip_model, an object for mip model of overt.

    Equality constraints in Overt are of two types:
        - linear equalities are an affine relation like: like z == 2x+3y
        - max equalities contain a max operation. They are of two kinds themselves:
                -- max-alone: like z = max(x, 0) which is basically a relu
                -- max-min:   like z = max(min(x,y), 0)
        linear equalities are processed in linear_2_mip as an additional constraints
        added to the model.
        max-alone equalities are processed in max_alone_2_mip. they are turned to a
        mip representation of a relu (following MIPVerify algorithm) and then
        added to the model
        max-min equalities are processed in max_min_2_mip. they are turned to a
        mip representation of a max-min, which is similar to the what is derived
        in MIPVerify algorithms for relu and max.
"""
function overt_2_mip(overt_mip_model::OvertMIP)
    # parsing equalities
    for eq in overt_mip_model.overt_app.approx_eq
        eq_2_mip(eq, overt_mip_model)
    end

    # parsing inequalities
    for ineq in overt_mip_model.overt_app.approx_ineq
        ineq_2_mip(ineq, overt_mip_model)
    end
end

"""
turn equality constraints to mip
"""
function eq_2_mip(expr::Expr, overt_mip_model::OvertMIP)
    #println(expr)
    fix_negation!(expr)
    if is_max_expr(expr)
        max_2_mip(expr, overt_mip_model)
    else
        affine_to_mip(expr, overt_mip_model)
    end
end

"""
turn inequality constraints to mip
"""
function ineq_2_mip(expr::Expr, overt_mip_model::OvertMIP)
    leftvar = get_mip_var(expr.args[2], overt_mip_model)
    ritevar = get_mip_var(expr.args[3], overt_mip_model)
    @debug "ineq_2_mip: Constraint is: $leftvar <= $ritevar"
    @constraint(overt_mip_model.model, leftvar <= ritevar )
end

"""
----------------------------------------------
low level functions: affine equalities
----------------------------------------------
"""

""" parse affine equalities z == 2*x + 3*y -1 """
function affine_to_mip(expr::Expr, overt_mip_model::OvertMIP)
    lhs = get_mip_var(expr.args[2], overt_mip_model)
    var_list, coeff_list, scalar = get_linear_coeffs(expr.args[3])
    for (v, c) in zip(var_list, coeff_list)
        mip_var = get_mip_var(v, overt_mip_model)
        lhs -= mip_var*c
    end
    lhs -= scalar
    @debug "affine_to_mip: Constraint is: $lhs == 0"
    @constraint(overt_mip_model.model, lhs == 0)
end

"""
----------------------------------------------
low level functions: max equalities
----------------------------------------------
"""

""" parse a big max into some of max-only or max-min expressions """
function max_2_mip(expr::Expr, overt_mip_model::OvertMIP)
    lhs = get_mip_var(expr.args[2], overt_mip_model)
    max_list = expr.args[3].args[2:end]
    for e in max_list
        @assert e.args[1] == :*
        c = e.args[2]
        tmp_var = get_mip_aux_var(overt_mip_model)
        lhs -= c*tmp_var
        if is_min_expr(e)
            max_min_2_mip(e.args[3], tmp_var, overt_mip_model)
        else
            max_alone_2_mip(e.args[3], tmp_var, overt_mip_model)
        end
    end
    @debug "max_2_mip: Constraint is: $lhs == 0"
    @constraint(overt_mip_model.model, lhs == 0)
end

""" parse expr like y == max(0, x) """
function max_alone_2_mip(expr::Expr, outvar::JuMP.VariableRef, overt_mip_model::OvertMIP)
    @assert expr.args[1] == :max
    @assert expr.args[2] == 0.0
    var, coef, scalar = get_linear_coeffs(expr.args[3])
    var, coef = var[1], coef[1]
    l, u = overt_mip_model.overt_app.ranges[var]
    if coef > 0
        L = l * coef + scalar
        U = u * coef + scalar
    else
        L = u * coef + scalar
        U = l * coef + scalar
    end
    var_mip = get_mip_var(var, overt_mip_model)
    a = get_mip_aux_var(overt_mip_model; binary=true)

    @constraint(overt_mip_model.model, outvar >= 0)
    @constraint(overt_mip_model.model, outvar >= var_mip*coef + scalar)
    @constraint(overt_mip_model.model, outvar <= a*U)
    @constraint(overt_mip_model.model, outvar <= var_mip*coef + scalar - L*(1-a))

    @debug "max_alone_2_mip: $outvar >= 0"
    @debug "max_alone_2_mip: $outvar >= $var_mip*$coef + $scalar"
    @debug "max_alone_2_mip: $outvar <= $a*$U"
    @debug "max_alone_2_mip: $outvar <= $var_mip*$coef + $scalar - $L*(1-$a)"
end

""" parse expr like z == max(0, min(x, y)) """
function max_min_2_mip(expr::Expr, outvar::JuMP.VariableRef, overt_mip_model::OvertMIP)
    @assert expr.args[1] == :max
    @assert expr.args[2] == 0
    min_arg = expr.args[3]
    @assert min_arg.args[1] == :min

    var1, coef1, scalar1 = get_linear_coeffs(min_arg.args[2])
    var1, coef1 = var1[1], coef1[1]
    l1, u1 = overt_mip_model.overt_app.ranges[var1]
    if coef1 > 0
        L1 = l1 * coef1 + scalar1
        U1 = u1 * coef1 + scalar1
    else
        L1 = u1 * coef1 + scalar1
        U1 = l1 * coef1 + scalar1
    end

    var2, coef2, scalar2 = get_linear_coeffs(min_arg.args[3])
    var2, coef2 = var2[1], coef2[1]
    l2, u2 = overt_mip_model.overt_app.ranges[var2]
    if coef2 > 0
        L2 = l2 * coef2 + scalar2
        U2 = u2 * coef2 + scalar2
    else
        L2 = u2 * coef2 + scalar2
        U2 = l2 * coef2 + scalar2
    end

    L3 = min(L1, L2)
    U3 = max(U1, U2)

    var1_mip = get_mip_var(var1, overt_mip_model)
    var2_mip = get_mip_var(var2, overt_mip_model)

    a  = get_mip_aux_var(overt_mip_model; binary=true)
    b  = get_mip_aux_var(overt_mip_model; binary=true)
    x3 = get_mip_aux_var(overt_mip_model)

    @constraint(overt_mip_model.model, x3 >= coef1*var1_mip + scalar1 - (1 - b) * (U1 - L2))
    @constraint(overt_mip_model.model, x3 >= coef2*var2_mip + scalar2 - b * (U2 - L1))
    @constraint(overt_mip_model.model, x3 <= coef1*var1_mip + scalar1)
    @constraint(overt_mip_model.model, x3 <= coef2*var2_mip + scalar2)

    @debug "max_min_2_mip: $x3 >= $coef1*$var1_mip + $scalar1 - (1 - $b) * ($U1 - $L2)"
    @debug "max_min_2_mip: $x3 >= $coef2*$var2_mip + $scalar2 - $b * ($U2 - $L1)"
    @debug "max_min_2_mip: $x3 <= $coef1*$var1_mip + $scalar1"
    @debug "max_min_2_mip: $x3 <= $coef2*$var2_mip + $scalar2"

    @constraint(overt_mip_model.model, outvar >= 0)
    @constraint(overt_mip_model.model, outvar >= x3)
    @constraint(overt_mip_model.model, outvar <= a*U3)
    @constraint(overt_mip_model.model, outvar <= x3 - L3*(1-a))

    @debug "max_min_2_mip: $outvar >= 0"
    @debug "max_min_2_mip: $outvar >= $x3"
    @debug "max_min_2_mip: $outvar <= $a*$U3"
    @debug "max_min_2_mip: $outvar <= $x3 - $L3*(1-$a)"
end

function mip_summary(model)
    MathOptInterface = MOI
    const_types = list_of_constraint_types(model)
    l_lin = 0
    l_bin = 0
	println(Crayon(foreground = :yellow), "="^50)
	println(Crayon(foreground = :yellow), "="^18 * " mip summary " * "="^19)
    for i = 1:length(const_types)
        var = const_types[i][1]
        const_type = const_types[i][2]
        l = length(all_constraints(model, var, const_type))
        #println("there are $l constraints of type $const_type with variables type $var.")
        if const_type != MathOptInterface.ZeroOne
            l_lin += l
        else
            l_bin += l
        end
    end
	println(Crayon(foreground = :yellow), "there are $l_lin linear constraints and $l_bin binary constraints.")
	println(Crayon(foreground = :yellow), "="^50)
    println(Crayon(foreground = :white), " ")
end


text:
import Base.issubset

"""
    issubset(X::LazySet, Y::LazySet, [witness]::Bool=false, args...)

Alias for `⊆` (inclusion check).

### Input

- `X`       -- set
- `Y`       -- set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X ⊆ Y``
* If `witness` option is activated:
  * `(true, [])` iff ``X ⊆ Y``
  * `(false, v)` iff ``X ⊈ Y`` and ``v ∈ X \\setminus Y``

### Notes

For more documentation see `⊆`.
"""
function issubset end

# this operation is forbidden, but it is a common error so we give a detailed error message
function ⊆(::AbstractVector, ::LazySet)
    throw(ArgumentError("cannot make an inclusion check if the left-hand side " *
          "is a vector; either wrap it as a set with one element, as in " *
          "`Singleton(v) ⊆ X`, or check for set membership, as in `v ∈ X` " *
          "(they behave equivalently although the implementations may differ)"))
end

# conversion for IA types
⊆(X::LazySet, Y::IA.Interval) = ⊆(X, Interval(Y))
⊆(X::IA.Interval, Y::LazySet) = ⊆(Interval(X), Y)

⊆(X::LazySet, Y::IA.IntervalBox) = ⊆(X, convert(Hyperrectangle, Y))
⊆(X::IA.IntervalBox, Y::LazySet) = ⊆(convert(Hyperrectangle, X), Y)

# --- fall-back with constraints_list of rhs ---

"""
    ⊆(X::LazySet, P::LazySet, [witness]::Bool=false)

Check whether a set is contained in a polyhedral set, and if not,
optionally compute a witness.

### Input

- `X`       -- inner set
- `Y`       -- outer polyhedral set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X ⊆ P``
* If `witness` option is activated:
  * `(true, [])` iff ``X ⊆ P``
  * `(false, v)` iff ``X ⊈ P`` and ``v ∈ X \\setminus P``

### Notes

We require that `constraints_list(P)` is available.

### Algorithm

We check inclusion of `X` in every constraint of `P`.
"""
function ⊆(X::LazySet, P::LazySet, witness::Bool=false)
    if applicable(constraints_list, P)
        return _issubset_constraints_list(X, P, witness)
    else
        error("an inclusion check for the given combination of set types is " *
              "not available")
    end
end

# --- AbstractHyperrectangle ---

"""
    ⊆(S::LazySet, H::AbstractHyperrectangle, [witness]::Bool=false)

Check whether a convex set is contained in a hyperrectangular set, and if not,
optionally compute a witness.

### Input

- `S` -- inner convex set
- `H` -- outer hyperrectangular set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``S ⊆ H``
* If `witness` option is activated:
  * `(true, [])` iff ``S ⊆ H``
  * `(false, v)` iff ``S ⊈ H`` and ``v ∈ S \\setminus H``

### Algorithm

``S ⊆ H`` iff ``\\operatorname{ihull}(S) ⊆ H``, where  ``\\operatorname{ihull}``
is the interval hull operator.
"""
function ⊆(S::LazySet, H::AbstractHyperrectangle, witness::Bool=false)
    return _issubset_in_hyperrectangle(S, H, witness)
end

function _issubset_in_hyperrectangle(S, H, witness)
    Shull = Approximations.interval_hull(S)
    return ⊆(Shull, H, witness)
end

# fix ambiguity
⊆(P::AbstractPolytope, H::AbstractHyperrectangle, witness::Bool=false) = _issubset_in_hyperrectangle(P, H, witness)

"""
    ⊆(H1::AbstractHyperrectangle, H2::AbstractHyperrectangle, [witness]::Bool=false)

Check whether a given hyperrectangular set is contained in another
hyperrectangular set, and if not, optionally compute a witness.

### Input

- `H1` -- inner hyperrectangular set
- `H2` -- outer hyperrectangular set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``H1 ⊆ H2``
* If `witness` option is activated:
  * `(true, [])` iff ``H1 ⊆ H2``
  * `(false, v)` iff ``H1 ⊈ H2`` and ``v ∈ H1 \\setminus H2``

### Algorithm

``H1 ⊆ H2`` iff ``c_1 + r_1 ≤ c_2 + r_2 ∧ c_1 - r_1 ≥ c_2 - r_2`` iff
``r_1 - r_2 ≤ c_1 - c_2 ≤ -(r_1 - r_2)``, where ``≤`` is taken component-wise.
"""
function ⊆(H1::AbstractHyperrectangle, H2::AbstractHyperrectangle, witness::Bool=false)
    @assert dim(H1) == dim(H2)
    N = promote_type(eltype(H1), eltype(H2))

    @inbounds for i in 1:dim(H1)
        c_dist = center(H1, i) - center(H2, i)
        r_dist = radius_hyperrectangle(H1, i) - radius_hyperrectangle(H2, i)
        # check if c_dist is not in the interval [r_dist, -r_dist]
        if !_leq(r_dist, c_dist) || !_leq(c_dist, -r_dist)
            if witness
                # compute a witness 'p' in the difference
                p = copy(center(H1))
                if c_dist >= zero(N)
                    p[i] += radius_hyperrectangle(H1, i)
                else
                    p[i] -= radius_hyperrectangle(H1, i)
                end
                return (false, p)
            else
                return false
            end
        end
    end

    if witness
        return (true, N[])
    else
        return true
    end
end

# --- AbstractPolytope ---

"""
    ⊆(P::AbstractPolytope, S::LazySet, [witness]::Bool=false;
      algorithm=_default_issubset(P, S))

Check whether a polytope is contained in a convex set, and if not, optionally
compute a witness.

### Input

- `P` -- inner polytope
- `S` -- outer convex set
- `witness`   -- (optional, default: `false`) compute a witness if activated
- `algorithm` -- (optional, default: `"constraints"` if the constraints list of `S`
                 is available, otherwise `"vertices"`) algorithm for the inclusion
                 check; available options are:

    * `"constraints"`, using the list of constraints of `P` and support function
      evaluations of `S`

    * `"vertices"`, using the list of vertices of `P` and membership evaluations
      of `S`

### Output

* If `witness` option is deactivated: `true` iff ``P ⊆ S``
* If `witness` option is activated:
  * `(true, [])` iff ``P ⊆ S``
  * `(false, v)` iff ``P ⊈ S`` and ``v ∈ P \\setminus S``

### Algorithm

Since ``S`` is convex, ``P ⊆ S`` iff ``v_i ∈ S`` for all vertices ``v_i`` of
``P``.
"""
function ⊆(P::AbstractPolytope, S::LazySet, witness::Bool=false;
           algorithm=_default_issubset(P, S))
    @assert dim(P) == dim(S)

    if algorithm == "constraints"
        return _issubset_constraints_list(P, S, witness)
    elseif algorithm == "vertices"
        return _issubset_vertices_list(P, S, witness)
    else
        error("algorithm $algorithm unknown")
    end
end

@inline function _default_issubset(P, S)
    if applicable(constraints_list, S)
        return "constraints"
    else
        return "vertices"
    end
end

# check whether P ⊆ S by testing if each vertex of P belongs to S
function _issubset_vertices_list(P, S, witness)
    for v in vertices_list(P)
        if v ∉ S
            if witness
                return (false, v)
            else
                return false
            end
        end
    end
    if witness
        return (true, N[])
    else
        return true
    end
end

"""
    ⊆(X::LazySet, P::AbstractPolyhedron, [witness]::Bool=false)

Check whether a convex set is contained in a polyhedron, and if not, optionally
compute a witness.

### Input

- `X` -- inner convex set
- `P` -- outer polyhedron (including a half-space)
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X ⊆ P``
* If `witness` option is activated:
  * `(true, [])` iff ``X ⊆ P``
  * `(false, v)` iff ``X ⊈ P`` and ``v ∈ P \\setminus X``

### Algorithm

Since ``X`` is convex, we can compare the support function of ``X`` and ``P`` in
each direction of the constraints of ``P``.

For witness generation, we use the support vector in the first direction where
the above check fails.
"""
function ⊆(X::LazySet, P::AbstractPolyhedron, witness::Bool=false)
    return _issubset_constraints_list(X, P, witness)
end

# for documentation see
# ⊆(X::LazySet, P::AbstractPolyhedron, witness::Bool=false)
function _issubset_constraints_list(S::LazySet, P::LazySet, witness::Bool=false)
    @assert dim(S) == dim(P)

    @inbounds for H in constraints_list(P)
        if !_leq(ρ(H.a, S), H.b)
            if witness
                return (false, σ(H.a, S))
            else
                return false
            end
        end
    end
    N = promote_type(eltype(S), eltype(P))
    return witness ? (true, N[]) : true
end

# disambiguations
for ST in [AbstractPolytope, AbstractHyperrectangle, LineSegment]
    @eval ⊆(X::($ST), P::AbstractPolyhedron, witness::Bool=false) = _issubset_constraints_list(X, P, witness)
end

# --- AbstractSingleton ---

"""
    ⊆(S::AbstractSingleton, X::LazySet, [witness]::Bool=false)

Check whether a given set with a single value is contained in a convex set, and
if not, optionally compute a witness.

### Input

- `S`       -- inner set with a single value
- `X`       -- outer convex set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``S ⊆ X``
* If `witness` option is activated:
  * `(true, [])` iff ``S ⊆ X``
  * `(false, v)` iff ``S ⊈ X`` and ``v ∈ S \\setminus X``
"""
function ⊆(S::AbstractSingleton, X::LazySet, witness::Bool=false)
    return _issubset_singleton(S, X, witness)
end

function _issubset_singleton(S, X, witness)
    result = element(S) ∈ X

    if witness
        N = promote_type(eltype(S), eltype(X))
        return (result, result ? N[] : element(S))
    else
        return result
    end
end

# disambiguations
for ST in [AbstractHyperrectangle, AbstractPolyhedron]
    @eval ⊆(X::AbstractSingleton, Y::($ST), witness::Bool=false) = _issubset_singleton(X, Y, witness)
end

"""
    ⊆(S1::AbstractSingleton, S2::AbstractSingleton, witness::Bool=false)

Check whether a given set with a single value is contained in another set with a
single value, and if not, optionally compute a witness.

### Input

- `S1` -- inner set with a single value
- `S2` -- outer set with a single value
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``S1 ⊆ S2`` iff ``S1 == S2``
* If `witness` option is activated:
  * `(true, [])` iff ``S1 ⊆ S2``
  * `(false, v)` iff ``S1 ⊈ S2`` and ``v ∈ S1 \\setminus S2``
"""
function ⊆(S1::AbstractSingleton, S2::AbstractSingleton, witness::Bool=false)
    result = element(S1) == element(S2) # TODO use _isapprox or isequivalent
    if witness
        N = promote_type(eltype(S1), eltype(S2))
        return (result, result ? N[] : element(S1))
    else
        return result
    end
end

# --- Ball2 ---

"""
    ⊆(B1::Ball2, B2::Ball2{N}, [witness]::Bool=false
     ) where {N<:AbstractFloat}

Check whether a ball in the 2-norm is contained in another ball in the 2-norm,
and if not, optionally compute a witness.

### Input

- `B1` -- inner ball in the 2-norm
- `B2` -- outer ball in the 2-norm
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``B1 ⊆ B2``
* If `witness` option is activated:
  * `(true, [])` iff ``B1 ⊆ B2``
  * `(false, v)` iff ``B1 ⊈ B2`` and ``v ∈ B1 \\setminus B2``

### Algorithm

``B1 ⊆ B2`` iff ``‖ c_1 - c_2 ‖_2 + r_1 ≤ r_2``
"""
function ⊆(B1::Ball2, B2::Ball2, witness::Bool=false)
    result = norm(B1.center - B2.center, 2) + B1.radius <= B2.radius
    if witness
        if result
            N = promote_type(eltype(B1), eltype(B2))
            return (result, N[])
        end
    else
        return result
    end

    # compute a witness 'v'
    v = B1.center .+ B1.radius * (B1.center .- B2.center)
    return (false, v)
end

# --- Ball2/Ballp ---

"""
    ⊆(B::Union{Ball2, Ballp}, S::AbstractSingleton, witness::Bool=false)

Check whether a ball in the 2-norm or p-norm is contained in a set with a single
value, and if not, optionally compute a witness.

### Input

- `B` -- inner ball in the 2-norm or p-norm
- `S` -- outer set with a single value
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``B ⊆ S``
* If `witness` option is activated:
  * `(true, [])` iff ``B ⊆ S``
  * `(false, v)` iff ``B ⊈ S`` and ``v ∈ B \\setminus S``
"""
function ⊆(B::Union{Ball2, Ballp}, S::AbstractSingleton, witness::Bool=false)
    result = B.center == element(S) && B.radius == 0 # TODO use _isapprox
    if witness
        if result
            N = promote_type(eltype(B), eltype(S))
            return (result, N[])
        end
    else
        return result
    end

    # compute a witness 'p' in the difference
    if B.center != element(S)
        p = B.center
    else
        p = copy(B.center)
        p[1] += B.radius
    end
    return (false, p)
end

# --- LineSegment ---

"""
    ⊆(L::LineSegment, S::LazySet, witness::Bool=false)

Check whether a line segment is contained in a convex set, and if not,
optionally compute a witness.

### Input

- `L` -- inner line segment
- `S` -- outer convex set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``L ⊆ S``
* If `witness` option is activated:
  * `(true, [])` iff ``L ⊆ S``
  * `(false, v)` iff ``L ⊈ S`` and ``v ∈ L \\setminus S``

### Algorithm

Since ``S`` is convex, ``L ⊆ S`` iff ``p ∈ S`` and ``q ∈ S``, where ``p, q`` are
the end points of ``L``.
"""
function ⊆(L::LineSegment, S::LazySet, witness::Bool=false)
    return _issubset_line_segment(L, S)
end

function _issubset_line_segment(L, S, witness)
    p_in_S = L.p ∈ S
    result = p_in_S && L.q ∈ S
    if !witness
        return result
    elseif result
        N = promote_type(eltype(L), eltype(S))
        return (result, N[])
    else
        return (result, p_in_S ? L.q : L.p)
    end
end

# fix ambiguity
⊆(L::LineSegment, H::AbstractHyperrectangle, witness::Bool=false) = _issubset_line_segment(L, H, witness)

# --- Interval ---

"""
    ⊆(x::Interval{N}, y::Interval, [witness]::Bool=false) where {N}

Check whether an interval is contained in another interval.

### Input

- `x`       -- interval
- `y`       -- interval
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

`true` iff ``x ⊆ y``.
"""
function ⊆(x::Interval{N}, y::Interval, witness::Bool=false) where {N}
    if min(y) > min(x)
        witness && return (false, low(x))
        return false
    elseif max(x) > max(y)
        witness && return (false, high(x))
        return false
    end
    witness && return (true, N[])
    return true
end

"""
    ⊆(x::Interval, U::UnionSet, [witness]::Bool=false)

Check whether an interval is contained in the union of two sets.

### Input

- `x` -- interval
- `U` -- union of two sets

### Output

`true` iff `x ⊆ U`.

### Notes

This implementation assumes that `U` is a union of one-dimensional convex sets.
Since these are equivalent to intervals, we convert to `Interval`s.

### Algorithm

Let ``U = a ∪ b `` where ``a`` and ``b`` are intervals and assume that the lower
bound of ``a`` is to the left of ``b``.
If the lower bound of ``x`` is to the left of ``a``, we have a counterexample.
Otherwise we compute the set difference ``y = x \\ a`` and check whether
``y ⊆ b`` holds.
"""
function ⊆(x::Interval, U::UnionSet, witness::Bool=false)
    @assert dim(U) == 1 "an interval is incompatible with a set of dimension " *
        "$(dim(U))"
    return _issubset_interval(x, convert(Interval, U.X), convert(Interval, U.Y),
                              witness)
end

function ⊆(x::Interval, U::UnionSet{N, <:Interval, <:Interval},
           witness::Bool=false) where {N}
    return _issubset_interval(x, U.X, U.Y, witness)
end

function _issubset_interval(x::Interval{N}, a::Interval, b::Interval,
                            witness) where {N}
    if min(a) > min(b)
        c = b
        b = a
        a = c
    end
    # a is on the left of b
    if min(x) < min(a)
        witness && return (false, low(x))
        return false
    end
    y = difference(x, a)
    if y ⊆ b
        witness && return (true, N[])
        return true
    elseif witness
        if min(b) > min(y)
            w = [(min(y) + min(b)) / 2]
        else
            w = high(y)
        end
        return (false, w)
    else
        return false
    end
end

function ⊆(x::Interval, U::UnionSetArray, witness::Bool=false)
    @assert dim(U) == 1 "an interval is incompatible with a set of dimension " *
        "$(dim(U))"
    V = _get_interval_array_copy(U)
    return _issubset_interval!(x, V, witness)
end

function _get_interval_array_copy(U::UnionSetArray)
    return [convert(Interval, X) for X in array(U)]
end

function _get_interval_array_copy(U::UnionSetArray{N, <:AbstractVector{<:Interval}}) where {N}
    return copy(array(U))
end

function _issubset_interval!(x::Interval{N}, intervals, witness) where {N}
    # sort intervals by lower bound
    sort!(intervals, lt=(x, y) -> min(x) <= min(y))

    # subtract intervals from x
    for y in intervals
        if min(y) > min(x)
            # lowest point of x is not contained
            witness && return (false, center(Interval(min(x), min(y))))
            return false
        end
        x = difference(x, y)
        if isempty(x)
            witness && return (true, N[])
            return true
        end
    end

    witness && return (false, low(x))
    return false
end

# --- EmptySet ---

"""
    ⊆(∅::EmptySet, X::LazySet, witness::Bool=false)

Check whether an empty set is contained in another set.

### Input

- `∅`       -- empty set
- `X`       -- another set
- `witness` -- (optional, default: `false`) compute a witness if activated
               (ignored, just kept for interface reasons)

### Output

`true`.
"""
function ⊆(∅::EmptySet, X::LazySet, witness::Bool=false)
    return _issubset_emptyset(∅, X, witness)
end

function _issubset_emptyset(∅, X, witness)
    N = promote_type(eltype(∅), eltype(X))
    return witness ? (true, N[]) : true
end

# disambiguations
⊆(∅::EmptySet, P::AbstractPolyhedron, witness::Bool=false) = _issubset_emptyset(∅, P, witness)
⊆(∅::EmptySet, H::AbstractHyperrectangle, witness::Bool=false) = _issubset_emptyset(∅, H, witness)

"""
    ⊆(X::LazySet, ∅::EmptySet, [witness]::Bool=false)

Check whether a set is contained in an empty set.

### Input

- `X`       -- another set
- `∅`       -- empty set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

`true` iff `X` is empty.

### Algorithm

We rely on `isempty(X)` for the emptiness check and on `an_element(X)` for
witness production.
"""
function ⊆(X::LazySet, ∅::EmptySet, witness::Bool=false)
    return _issubset_in_emptyset(X, ∅, witness)
end

function _issubset_in_emptyset(X, ∅, witness)
    if isempty(X)
        N = promote_type(eltype(∅), eltype(X))
        return witness ? (true, N[]) : true
    else
        return witness ? (false, an_element(X)) : false
    end
end

# disambiguations
⊆(X::AbstractPolytope, ∅::EmptySet, witness::Bool=false) = _issubset_in_emptyset(X, ∅, witness)

function ⊆(X::AbstractSingleton, ::EmptySet, witness::Bool=false)
    return witness ? (false, an_element(X)) : false
end

function ⊆(X::LineSegment, ::EmptySet, witness::Bool=false)
    return witness ? (false, an_element(X)) : false
end

function ⊆(∅₁::EmptySet, ∅₂::EmptySet, witness::Bool=false)
    N = promote_type(eltype(∅₁), eltype(∅₂))
    return witness ? (true, N[]) : true
end

# --- UnionSet ---

"""
    ⊆(cup::UnionSet, X::LazySet, [witness]::Bool=false)

Check whether a union of two convex sets is contained in another set.

### Input

- `cup`     -- union of two convex sets
- `X`       -- another set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``\\text{cup} ⊆ X``
* If `witness` option is activated:
  * `(true, [])` iff ``\\text{cup} ⊆ X``
  * `(false, v)` iff ``\\text{cup} \\not\\subseteq X`` and
    ``v ∈ \\text{cup} \\setminus X``
"""
function ⊆(cup::UnionSet, X::LazySet, witness::Bool=false)
    return ⊆(UnionSetArray([cup.X, cup.Y]), X, witness) # TODO implement here
end

"""
    ⊆(cup::UnionSetArray, X::LazySet, [witness]::Bool=false)

Check whether a union of a finite number of convex sets is contained in another
set.

### Input

- `cup`     -- union of a finite number of convex sets
- `X`       -- another set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``\\text{cup} ⊆ X``
* If `witness` option is activated:
  * `(true, [])` iff ``\\text{cup} ⊆ X``
  * `(false, v)` iff ``\\text{cup} \\not\\subseteq X`` and
    ``v ∈ \\text{cup} \\setminus X``
"""
function ⊆(cup::UnionSetArray, X::LazySet, witness::Bool=false)
    result = true
    N = promote_type(eltype(cup), eltype(X))
    w = N[]
    for Y in array(cup)
        if witness
            result, w = ⊆(Y, X, witness)
        else
            result = ⊆(Y, X, witness)
        end
        if !result
            break
        end
    end
    return witness ? (result, w) : result
end

# --- Universe ---

"""
    ⊆(X::LazySet, U::Universe, [witness]::Bool=false)

Check whether a convex set is contained in a universe.

### Input

- `U`       -- universe
- `X`       -- convex set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true`
* If `witness` option is activated: `(true, [])`
"""
function ⊆(X::LazySet, U::Universe, witness::Bool=false)
    return _issubset_universe(X, U, witness)
end

function _issubset_universe(X, U, witness)
    N = promote_type(eltype(X), eltype(U))
    witness ? (true, N[]) : true
end

# disambiguations
for ST in [AbstractPolytope, AbstractHyperrectangle, AbstractSingleton, LineSegment, EmptySet]
    @eval ⊆(X::($ST), U::Universe, witness::Bool=false) = _issubset_universe(X, U, witness)
end

"""
    ⊆(U::Universe, X::LazySet, [witness]::Bool=false)

Check whether a universe is contained in another convex set, and otherwise
optionally compute a witness.

### Input

- `U`       -- universe
- `X`       -- convex set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``U ⊆ X``
* If `witness` option is activated:
  * `(true, [])` iff ``U ⊆ X``
  * `(false, v)` iff ``U \\not\\subseteq X`` and
    ``v ∈ U \\setminus X``

### Algorithm

We fall back to `isuniversal(X)`.
"""
function ⊆(U::Universe, X::LazySet, witness::Bool=false)
    return isuniversal(X, witness)
end

# disambiguations
function ⊆(U1::Universe, U2::Universe, witness::Bool=false)
    N = promote_type(eltype(U1), eltype(U2))
    return witness ? (true, N[]) : true
end

# disambiguations
for ST in [AbstractPolyhedron, AbstractPolytope, AbstractHyperrectangle, AbstractSingleton, EmptySet]
    @eval ⊆(U::Universe, X::($ST), witness::Bool=false) = isuniversal(X, witness)
end

# --- Complement ---

"""
    ⊆(X::LazySet, C::Complement, [witness]::Bool=false)

Check whether a convex set is contained in the complement of another convex set,
and otherwise optionally compute a witness.

### Input

- `X`       -- convex set
- `C`       -- complement of a convex set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X ⊆ C``
* If `witness` option is activated:
  * `(true, [])` iff ``X ⊆ C``
  * `(false, v)` iff ``X \\not\\subseteq C`` and
    ``v ∈ X \\setminus C``

### Algorithm

We fall back to `isdisjoint(X, C.X)`, which can be justified as follows.

```math
    X ⊆ Y^C ⟺ X ∩ Y = ∅
```
"""
function ⊆(X::LazySet, C::Complement, witness::Bool=false)
    return isdisjoint(X, C.X, witness)
end

# --- CartesianProduct ---

"""
    ⊆(X::CartesianProduct, Y::CartesianProduct, [witness]::Bool=false;
      check_block_equality::Bool=true)

Check whether a Cartesian product of two convex sets is contained in another
Cartesian product of two convex sets, and otherwise optionally compute a
witness.

### Input

- `X`       -- Cartesian product of two convex sets
- `Y`       -- Cartesian product of two convex sets
- `witness` -- (optional, default: `false`) compute a witness if activated
- `check_block_equality` -- (optional, default: `true`) flag for checking that
                            the block structure of the two sets is identical

### Output

* If `witness` option is deactivated: `true` iff ``X ⊆ Y``
* If `witness` option is activated:
  * `(true, [])` iff ``X ⊆ Y``
  * `(false, v)` iff ``X \\not\\subseteq Y`` and
    ``v ∈ X \\setminus Y``

### Notes

This algorithm requires that the two Cartesian products share the same block
structure.
If `check_block_equality` is activated, we check this property and, if it does
not hold, we use a fallback implementation based on conversion to constraint
representation (assuming that the sets are polyhedral).

### Algorithm

We check for inclusion for each block of the Cartesian products.

For witness production, we obtain a witness in one of the blocks.
We then construct a high-dimensional witness by obtaining any point in the other
blocks (using `an_element`) and concatenating these points.
"""
function ⊆(X::CartesianProduct, Y::CartesianProduct, witness::Bool=false;
           check_block_equality::Bool=true)
    n1 = dim(X.X)
    n2 = dim(X.Y)
    if check_block_equality && (n1 != dim(Y.X) || n2 != dim(Y.Y))
        return _issubset_constraints_list(X, Y, witness)
    end

    # check first block
    result = ⊆(X.X, Y.X, witness)
    if !witness && !result
        return false
    elseif witness && !result[1]
        # construct a witness
        w = vcat(result[2], an_element(X.Y))
        return (false, w)
    end

    # check second block
    result = ⊆(X.Y, Y.Y, witness)
    if !witness && !result
        return false
    elseif witness && !result[1]
        # construct a witness
        w = vcat(an_element(X.X), result[2])
        return (false, w)
    end

    N = promote_type(eltype(X), eltype(Y))
    return witness ? (true, N[]) : true
end

"""
    ⊆(X::CartesianProductArray, Y::CartesianProductArray, [witness]::Bool=false;
      check_block_equality::Bool=true)

Check whether a Cartesian product of finitely many convex sets is contained in
another Cartesian product of finitely many convex sets, and otherwise optionally
compute a witness.

### Input

- `X`       -- Cartesian product of finitely many convex sets
- `Y`       -- Cartesian product of finitely many convex sets
- `witness` -- (optional, default: `false`) compute a witness if activated
- `check_block_equality` -- (optional, default: `true`) flag for checking that
                             the block structure of the two sets is identical

### Output

* If `witness` option is deactivated: `true` iff ``X ⊆ Y``
* If `witness` option is activated:
  * `(true, [])` iff ``X ⊆ Y``
  * `(false, v)` iff ``X \\not\\subseteq Y`` and
    ``v ∈ X \\setminus Y``

### Notes

This algorithm requires that the two Cartesian products share the same block
structure.
If `check_block_equality` is activated, we check this property and, if it does
not hold, we use a fallback implementation based on conversion to constraint
representation (assuming that the sets are polyhedral).

### Algorithm

We check for inclusion for each block of the Cartesian products.

For witness production, we obtain a witness in one of the blocks.
We then construct a high-dimensional witness by obtaining any point in the other
blocks (using `an_element`) and concatenating these points.
"""
function ⊆(X::CartesianProductArray, Y::CartesianProductArray, witness::Bool=false;
           check_block_equality::Bool=true)
    aX = array(X)
    aY = array(Y)
    if check_block_equality && !same_block_structure(aX, aY)
        return _issubset_constraints_list(X, Y, witness)
    end

    N = promote_type(eltype(X), eltype(Y))
    for i in 1:length(aX)
        result = ⊆(aX[i], aY[i], witness)
        if !witness && !result
            return false
        elseif witness && !result[1]
            # construct a witness
            w = Vector{N}(undef, dim(X))
            k = 1
            for j in 1:length(aX)
                Xj = aX[j]
                l = k + dim(Xj)
                w[k:l-1] = j == i ? result[2] : an_element(Xj)
                k = l
            end
            return (false, w)
        end
    end

    return witness ? (true, N[]) : true
end

# --- AbstractZonotope ---

"""
    ⊆(Z::AbstractZonotope, H::AbstractHyperrectangle, [witness]::Bool=false)

Check whether a zonotopic set is contained in a hyperrectangular set.

### Input

- `Z`       -- inner zonotopic set
- `H`       -- outer hyperrectangular set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

`true` iff ``Z ⊆ H`` otherwise `false`

### Algorithm

Algorithm based on Lemma 3.1 of [1]

[1] Mitchell, I. M., Budzis, J., & Bolyachevets, A. (2019, April). Invariant,
 viability and discriminating kernel under-approximation via zonotope scaling.
 In Proceedings of the 22nd ACM International Conference on Hybrid Systems:
 Computation and Control (pp. 268-269).
"""
function ⊆(Z::AbstractZonotope, H::AbstractHyperrectangle, witness::Bool=false)
    witness && raise(ValueError("witness production is not supported yet"))

    c = center(Z)
    G = genmat(Z)
    n, m = size(G)
    N = promote_type(eltype(Z), eltype(H))
    @inbounds for i = 1:n
        aux = zero(N)
        for j in 1:m
            aux += abs(G[i, j])
        end
        ubound = c[i] + aux
        lbound = c[i] - aux
        if !_leq(ubound, high(H, i)) || !_geq(lbound, low(H, i))
            return false
        end
    end
    return true
end


text:
import LinearAlgebra

function relative_error_l2(a, b)
  return LinearAlgebra.norm(a - b)^2 / LinearAlgebra.norm(b)^2
end

# @testset "$(format("SGD Optimizer: Bernoulli Loss [Small][Forgiving][Native]"))" begin
#   let
#     input_size = 1600* 400
#     y = rand(Distributions.Gamma(10, 2), input_size) * 1.0
#     mle_x = MatrixCompletion.Losses.sgd_train(Loss{Gamma}(),
#                                               fx         = y,
#                                               y          = y,
#                                               c          = zeros(input_size),
#                                               ρ          = 0.2,
#                                               α          = 0.2,
#                                               ρ₁         = 0.9,
#                                               ρ₂         = 0.999,
#                                               batch_size = 1600 * 400,
#                                               epoch      = 20);
#     prediction = predict(Gamma(),forward_map(Gamma(),mle_x))
#     # @show(relative_error_l2(prediction, y))
#     tc = provide(Diagnostics{Poisson()}(),
#                  input_data=prediction, reference=y)
#     display(tc)
#   end
# end

@testset "$(format("GD Optimizer: Gamma Loss [Small][Forgiving][Native]"))" begin
  let
    input_size = 200 * 200
    y = rand(Distributions.Gamma(10, 2), input_size) * 1.0
    mle_x = MatrixCompletion.Losses.train(Loss{Gamma}(),
                                          fx      = rand(input_size),
                                          y       = y,
                                          c       = zeros(input_size),
                                          ρ       = 0,
                                          γ       = 0.2,
                                          iter    = 500,
                                          verbose = true)
    @show(mle_x[1:20])
    prediction = predict(Gamma(),forward_map(Gamma(),mle_x))
    tc = provide(Diagnostics{Poisson()}(),
                 input_data=prediction, reference=y)
    display(tc)
  end
end