text:
################################################################################
#
#  NfRel/NfRel.jl : Relative number field extensions
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: 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.
#
#
#  Copyright (C) 2017 Claus Fieker, Tommy Hofmann
#
################################################################################

export absolute_field

################################################################################
#
#  Copy
#
################################################################################

function Base.deepcopy_internal(a::NfRelElem{T}, dict::IdDict) where T
  z = NfRelElem{T}(Base.deepcopy_internal(data(a), dict))
  z.parent = parent(a)
  return z
end

################################################################################
#
#  Comply with Nemo ring interface
#
################################################################################

Nemo.elem_type(::Type{NfRel{T}}) where {T} = NfRelElem{T}

Nemo.elem_type(::NfRel{T}) where {T} = NfRelElem{T}

Nemo.parent_type(::Type{NfRelElem{T}}) where {T} = NfRel{T}

order_type(K::NfRel{T}) where {T} = NfRelOrd{T, frac_ideal_type(order_type(base_field(K)))}

order_type(::Type{NfRel{T}}) where {T} = NfRelOrd{T, frac_ideal_type(order_type(parent_type(T)))}

Nemo.needs_parentheses(::NfRelElem) = true

Nemo.isnegative(x::NfRelElem) = Nemo.isnegative(data(x))

Nemo.show_minus_one(::Type{NfRelElem{T}}) where {T} = true

function Nemo.iszero(a::NfRelElem)
  reduce!(a)
  return iszero(data(a))
end

function Nemo.isone(a::NfRelElem)
  reduce!(a)
  return isone(data(a))
end

Nemo.zero(K::NfRel) = K(Nemo.zero(parent(K.pol)))

Nemo.one(K::NfRel) = K(Nemo.one(parent(K.pol)))

function Nemo.zero!(a::NfRelElem)
  Nemo.zero!(a.data)
  return a
end

################################################################################
#
#  Promotion
#
################################################################################

Nemo.promote_rule(::Type{NfRelElem{S}}, ::Type{T}) where {T <: Integer, S} = NfRelElem{S}

Nemo.promote_rule(::Type{NfRelElem{T}}, ::Type{fmpz}) where {T <: Nemo.RingElement} = NfRelElem{T}

Nemo.promote_rule(::Type{NfRelElem{T}}, ::Type{fmpq}) where {T <: Nemo.RingElement} = NfRelElem{T}

Nemo.promote_rule(::Type{NfRelElem{T}}, ::Type{T}) where {T} = NfRelElem{T}

Nemo.promote_rule(::Type{NfRelElem{T}}, ::Type{NfRelElem{T}}) where T <: Nemo.RingElement = NfRelElem{T}

function Nemo.promote_rule(::Type{NfRelElem{T}}, ::Type{U}) where {T <: Nemo.RingElement, U <: Nemo.RingElement}
  Nemo.promote_rule(T, U) == T ? NfRelElem{T} : Union{}
end

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

@inline Nemo.base_ring(a::NfRel{T}) where {T} = a.base_ring::parent_type(T)

@inline base_field(a::NfRel{T}) where {T} = a.base_ring::parent_type(T)

@inline Nemo.data(a::NfRelElem) = a.data

@inline Nemo.parent(a::NfRelElem{T}) where {T} = a.parent::NfRel{T}

issimple(a::NfRel) = true

################################################################################
#
#  Reduction
#
################################################################################

function reduce!(a::NfRelElem)
  a.data = mod(a.data, parent(a).pol)
  return a
end

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

function Base.show(io::IO, a::NfRel)
  print(io, "Relative number field over\n")
  print(io, a.base_ring, "\n")
  print(io, " with defining polynomial ", a.pol)
end

function _show(io::IO, x::PolyElem, S::String)
   len = length(x)
   if len == 0
      print(io, base_ring(x)(0))
   else
      for i = 1:len - 1
         c = coeff(x, len - i)
         bracket = needs_parentheses(c)
         if !iszero(c)
            if i != 1 && !isnegative(c)
               print(io, "+")
            end
            if !isone(c) && (c != -1 || show_minus_one(typeof(c)))
               if bracket
                  print(io, "(")
               end
               show(io, c)
               if bracket
                  print(io, ")")
               end
               print(io, "*")
            end
            if c == -1 && !show_minus_one(typeof(c))
               print(io, "-")
            end
            print(io, string(S))
            if len - i != 1
               print(io, "^")
               print(io, len - i)
            end
         end
      end
      c = coeff(x, 0)
      bracket = needs_parentheses(c)
      if !iszero(c)
         if len != 1 && !isnegative(c)
            print(io, "+")
         end
         if bracket
            print(io, "(")
         end
         show(io, c)
         if bracket
            print(io, ")")
         end
      end
   end
end

function Base.show(io::IO, a::NfRelElem)
  f = data(a)
  _show(io, f, string(parent(a).S))
end

################################################################################
#
#  Constructors and parent object overloading
#
################################################################################

@doc Markdown.doc"""
    NumberField(f::Generic.Poly{T}, s::String; cached::Bool = false, check::Bool = false) where T
Given an irreducible polynomial $f$ over some number field $K$,
create the field $K[t]/f$.
$f$ must be irreducible - although this is not tested.
"""
function NumberField(f::Generic.Poly{T}, s::String; cached::Bool = false, check::Bool = false) where T
  S = Symbol(s)
  check && !isirreducible(f) && error("Polynomial must be irreducible")
  K = NfRel{T}(f, S, cached)
  return K, K(gen(parent(f)))
end

@doc Markdown.doc"""
    NumberField(f::Generic.Poly{T}; cached::Bool = false, check::Bool = false) where T
Given an irreducible polynomial $f$ over some number field $K$,
create the field $K[t]/f$.
$f$ must be irreducible - although this is not tested.
"""
function NumberField(f::Generic.Poly{T}; cached::Bool = false, check::Bool = false) where T
  return NumberField(f, "_\$", cached = cached, check = check)
end

function (K::NfRel{T})(a::Generic.Poly{T}) where T
  z = NfRelElem{T}(mod(a, K.pol))
  z.parent = K
  return z
end

function (K::NfRel{T})(a::T) where T
  parent(a) != base_ring(parent(K.pol)) == error("Cannot coerce")
  z = NfRelElem{T}(parent(K.pol)(a))
  z.parent = K
  return z
end

function (K::NfRel{T})(a::Vector{T}) where T
  @assert length(a) <= degree(K)
  z = NfRelElem{T}(parent(K.pol)(a))
  z.parent = K
  return z
end

(K::NfRel)(a::Integer) = K(parent(K.pol)(a))

(K::NfRel)(a::Rational{T}) where {T <: Integer} = K(parent(K.pol)(a))

(K::NfRel)(a::fmpz) = K(parent(K.pol)(a))

(K::NfRel)(a::fmpq) = K(parent(K.pol)(a))

(K::NfRel)() = zero(K)

Nemo.gen(K::NfRel) = K(Nemo.gen(parent(K.pol)))

################################################################################
#
#  Unary operators
#
################################################################################

function Base.:(-)(a::NfRelElem)
  return parent(a)(-data(a))
end

################################################################################
#
#  Binary operators
#
################################################################################

function Base.:(+)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  return parent(a)(data(a) + data(b))
end

function Base.:(-)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  return parent(a)(data(a) - data(b))
end

function Base.:(*)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  return parent(a)(data(a) * data(b))
end

function Nemo.divexact(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  b == 0 && error("Element not invertible")
  return a*inv(b)
end

Base.:(//)(a::NfRelElem{T}, b::NfRelElem{T}) where {T} = divexact(a, b)

################################################################################
#
#  Inversion
#
################################################################################

function Base.inv(a::NfRelElem)
  a == 0 && error("Element not invertible")
  g, s, _ = gcdx(data(a), parent(a).pol)
  @assert g == 1
  return parent(a)(s)
end

################################################################################
#
#  Powering
#
################################################################################

Base.:(^)(a::NfRelElem, n::UInt) = a^Int(n)

function Base.:(^)(a::NfRelElem, n::Int)
  K = parent(a)
  if iszero(a)
    return zero(K)
  end

  if n == 0
    return one(K)
  end

  if n < 0 && iszero(a)
    error("Element is not invertible")
  end

  return K(powmod(data(a), n, K.pol))
end

function Base.:(^)(a::NfRelElem, b::fmpz)
  if b < 0
    return inv(a)^(-b)
  elseif b == 0
    return parent(a)(1)
  elseif b == 1
    return deepcopy(a)
  elseif mod(b, 2) == 0
    c = a^(div(b, 2))
    return c*c
  elseif mod(b, 2) == 1
    return a^(b - 1)*a
  end
end

################################################################################
#
#  Comparison
#
################################################################################

function Base.:(==)(a::NfRelElem{T}, b::NfRelElem{T}) where T
  reduce!(a)
  reduce!(b)
  return data(a) == data(b)
end

if !isdefined(Nemo, :promote_rule1)
  function Base.:(==)(a::NfRelElem{T}, b::Union{Integer, Rational}) where T
    return a == parent(a)(b)
  end
end

################################################################################
#
#  Unsafe operations
#
################################################################################

function Nemo.mul!(c::NfRelElem{T}, a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  mul!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

function mul!(c::NfRelElem{T}, a::NfRelElem{T}, b::T) where {T}
  mul!(c.data, a.data, b)
  c = reduce!(c)
  return c
end

function Nemo.addeq!(b::NfRelElem{T}, a::NfRelElem{T}) where {T}
  addeq!(b.data, a.data)
  b = reduce!(b)
  return b
end

function Nemo.add!(c::NfRelElem{T}, a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  add!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

################################################################################
#
#  Isomorphism to absolute field (AnticNumberField)
#
################################################################################

@doc Markdown.doc"""
    absolute_field(K::NfRel{nf_elem}, cached::Bool = false) -> AnticNumberField, Map, Map
Given an extension $K/k/Q$, find an isomorphic extension of $Q$.
"""
function absolute_field(K::NfRel{nf_elem}, cached::Bool = false)
  Ka, a, b, c = _absolute_field(K, cached)
  return Ka, NfRelToNf(K, Ka, a, b, c), hom(base_ring(K), Ka, a, check = false)
end

@doc Markdown.doc"""
    absolute_field(K::NfRel{NfRelElem}, cached::Bool = false) -> NfRel, Map, Map
Given an extension $E/K/k$, find an isomorphic extension of $k$.
In a tower, only the top-most steps are collapsed.
"""
function absolute_field(K::NfRel{NfRelElem{T}}, cached::Bool = false) where T
  Ka, a, b, c = _absolute_field(K)
  return Ka, NfRelRelToNfRel(K, Ka, a, b, c), hom(base_ring(K), Ka, a, check = false)
end

#Trager: p4, Algebraic Factoring and Rational Function Integration
function _absolute_field(K::NfRel, cached::Bool = false)
  f = K.pol
  kx = parent(f)
  k = base_ring(kx)
  Qx = parent(k.pol)

  l = 0
  g = f
  N = norm(g)

  while true
    @assert degree(N) == degree(g) * degree(k)

    if !isconstant(N) && issquarefree(N)
      break
    end

    l += 1

    g = compose(f, gen(kx) - l*gen(k))
    N = norm(g)
  end

  Ka = NumberField(N, "x", cached = cached, check = false)[1]

  KaT, T = PolynomialRing(Ka, "T", cached = false)

  # map Ka -> K: gen(Ka) -> gen(K)+ k gen(k)

  # gen(k) -> Root(gcd(g, poly(k)))  #gcd should be linear:
  # g in kx = (Q[a])[x]. Want to map x -> gen(Ka), a -> T

  gg = zero(KaT)
  for i=degree(g):-1:0
    gg = gg*gen(Ka) + Qx(coeff(g, i))(gen(KaT))
  end

  q = gcd(gg, k.pol(gen(KaT)))
  @assert degree(q) == 1
  al = -trailing_coefficient(q)//lead(q)
  be = gen(Ka) - l*al
  ga = gen(K) + l*gen(k)

  #al -> gen(k) in Ka
  #be -> gen(K) in Ka
  #ga -> gen(Ka) in K
  return Ka, al, be, ga
end

function Nemo.check_parent(a, b)
  return a==b
end

function Nemo.content(a::Generic.Poly{T}) where T <: Nemo.Field
  return base_ring(a)(1)
end

Nemo.canonical_unit(a::NfRelElem) = a

#function +(a::NfRelElem{NfRelElem{T}}, b::NfRelElem{T}) where T
#  c = deepcopy(a)
#  setcoeff!(c.data, 0, coeff(c.data, 0)+b)
#  return c
#end
#
#+(a::NfRelElem{T}, b::NfRelElem{NfRelElem{T}}) where T = b+a
#
#function *(a::NfRelElem{NfRelElem{T}}, b::NfRelElem{T}) where T
#  c = deepcopy(a)
#  setcoeff!(c.data, 0, coeff(c.data, 0)*b)
#  return c
#end
#
#*(a::NfRelElem{T}, b::NfRelElem{NfRelElem{T}}) where T = b*a

@inline coeff(a::NfRelElem{T}, i::Int) where {T} = coeff(a.data, i)

@inline setcoeff!(a::NfRelElem{T}, i::Int, c::T) where {T} = setcoeff!(a.data, i, c)

@inline degree(L::Hecke.NfRel) = degree(L.pol)

#######################################################################
# convenience constructions
#######################################################################
@doc Markdown.doc"""
    ispure_extension(K::NfRel) -> Bool
Tests if $K$ is pure, ie. if the defining polynomial is $x^n-g$.
"""
function ispure_extension(K::NfRel)
  if !ismonic(K.pol)
    return false
  end
  return all(i->iszero(coeff(K.pol, i)), 1:degree(K)-1)
end

@doc Markdown.doc"""
    iskummer_extension(K::Hecke.NfRel{nf_elem}) -> Bool
Tests if $K$ is Kummer, ie. if the defining polynomial is $x^n-g$ and
if the coefficient field contains the $n$-th roots of unity.
"""
function iskummer_extension(K::Hecke.NfRel{nf_elem})
  if !ispure_extension(K)
    return false
  end

  k = base_ring(K)
  Zk = maximal_order(k)
  _, o = Hecke.torsion_units_gen_order(Zk)
  if o % degree(K) != 0
    return false
  end
  return true
end

@doc Markdown.doc"""
    pure_extension(n::Int, gen::FacElem{nf_elem, AnticNumberField}) -> NfRel{nf_elem}, NfRelElem
    pure_extension(n::Int, gen::nf_elem) -> NfRel{nf_elem}, NfRelElem
Create the field extension with the defining polynomial $x^n-gen$.
"""
function pure_extension(n::Int, gen::FacElem{nf_elem, AnticNumberField})
  return pure_extension(n, evaluate(gen))
end

function pure_extension(n::Int, gen::nf_elem)
  k = parent(gen)
  kx, x = PolynomialRing(k, cached = false)
  return number_field(x^n-gen)
end

function hash(a::Hecke.NfRelElem{nf_elem}, b::UInt)
  return hash(a.data, b)
end

################################################################################
#
#  Representation Matrix
#
################################################################################

function basis_mat(v::Vector{<: NfRelElem})
  K = parent(v[1])
  k = base_field(K)
  z = zero_matrix(k, length(v), degree(K))
  for i in 1:length(v)
    for j in 0:(degree(K) - 1)
      z[i, j + 1] = coeff(v[i], j)
    end
  end
  return z
end

function elem_to_mat_row!(M::Generic.Mat{T}, i::Int, a::NfRelElem{T}) where T
  for c = 1:ncols(M)
    M[i, c] = deepcopy(coeff(a, c - 1))
  end
  return nothing
end

function elem_from_mat_row(L::NfRel{T}, M::Generic.Mat{T}, i::Int) where T
  t = L(1)
  a = L()
  for c = 1:ncols(M)
    a += M[i, c]*t
    mul!(t, t, gen(L))
  end
  return a
end

function representation_matrix(a::NfRelElem)
  L = a.parent
  n = degree(L)
  M = zero_matrix(base_ring(L), n, n)
  t = gen(L)
  b = deepcopy(a)
  for i = 1:n-1
    elem_to_mat_row!(M, i, b)
    mul!(b, b, t)
  end
  elem_to_mat_row!(M, n, b)
  return M
end

function norm(a::NfRelElem{nf_elem}, new::Bool = !true)
  if new && ismonic(parent(a).pol) #should be much faster - eventually
    return resultant_mod(a.data, parent(a).pol)
  end
  M = representation_matrix(a)
  return det(M)
end

function norm(a::NfRelElem, new::Bool = true)
  if new && ismonic(parent(a).pol)
    return resultant(a.data, parent(a).pol)
  end
  M = representation_matrix(a)
  return det(M)
end

function tr(a::NfRelElem)
  M = representation_matrix(a)
  return tr(M)
end

################################################################################
#
#  Random elements from arrays
#
################################################################################

function rand!(c::NfRelElem, b::Vector{<: NfRelElem}, r::UnitRange)
  K = parent(b[1])
  length(b) == 0 && error("Array must not be empty")

  # TODO: Avoid promotion to K
  mul!(c, b[1], K(rand(r)))
  t = zero(K)

  for i = 2:length(b)
    mul!(t, b[i], K(rand(r)))
    add!(c, t, c)
  end

  return c
end

######################################################################
# fun in towers..
######################################################################

isunit(a::NfRelElem) = !iszero(a)

absolute_degree(A::AnticNumberField) = degree(A)

function absolute_degree(A::NfRel)
  return absolute_degree(base_ring(A))*degree(A)
end

function tr(a::NfRelElem, k::Union{NfRel, AnticNumberField, FlintRationalField})
  b = tr(a)
  while parent(b) != k
    b = tr(b)
  end
  return b
end

function norm(a::NfRelElem, k::Union{NfRel, AnticNumberField, FlintRationalField})
  b = norm(a)
  while parent(b) != k
    b = norm(b)
  end
  return b
end

function absolute_tr(a::NfRelElem)
  return tr(a, FlintQQ)
end

function absolute_norm(a::NfRelElem)
  return norm(a, FlintQQ)
end

#TODO: investigate charpoly/ minpoly from power_sums, aka tr(a^i) and
#      Newton identities
#TODO: cache traces of powers of the generator on the field, then
#      the trace does not need the matrix

@doc Markdown.doc"""
    charpoly(a::NfRelElem) -> PolyElem

Given an element $a$ in an extension $L/K$, this function returns the
characteristic polynomial of $a$ over $K$.
"""
function charpoly(a::NfRelElem)
  M = representation_matrix(a)
  R = PolynomialRing(base_ring(parent(a)), cached = false)[1]
  return minpoly(R, M, true)
end

@doc Markdown.doc"""
    minpoly(a::NfRelElem) -> PolyElem

Given an element $a$ in an extension $L/K$, this function returns the minimal
polynomial of $a$ of $K$.
"""
function minpoly(a::NfRelElem)
  M = representation_matrix(a)
  R = PolynomialRing(base_ring(parent(a)), cached = false)[1]
  return minpoly(R, M, false)
end

function charpoly(a::NfRelElem, k::Union{NfRel, AnticNumberField, FlintRationalField})
  f = charpoly(a)
  while base_ring(f) != k
    f = norm(f)
  end
  return f
end

function absolute_charpoly(a::NfRelElem)
  return charpoly(a, FlintQQ)
end

function minpoly(a::NfRelElem, k::Union{NfRel, AnticNumberField, FlintRationalField})

  if parent(a) == k
    R, y  = PolynomialRing(k, cached = false)
    return y - a
  end

  f = minpoly(a)
    while base_ring(f) != k
    @show typeof(f)
    f = norm(f)
    g = gcd(f, derivative(f))
    if !isone(g)
      f = divexact(f, g)
    end
  end
  return f
end

function absolute_minpoly(a::NfRelElem)
  return minpoly(a, FlintQQ)
end

#

(K::NfRel)(x::NfRelElem) = K(base_ring(K)(x))

(K::NfRel)(x::nf_elem) = K(base_ring(K)(x))

(K::NfRel{T})(x::NfRelElem{T}) where {T} = K(x.data)

(K::NfRel{NfRelElem{T}})(x::NfRelElem{T}) where {T} = K(parent(K.pol)(x))

(K::NfRel{nf_elem})(x::nf_elem) = K(parent(K.pol)(x))

#

function (R::Generic.PolyRing{nf_elem})(a::NfRelElem{nf_elem})
  if base_ring(R)==base_ring(parent(a))
    return a.data
  end
  error("wrong ring")
end

function (R::Generic.PolyRing{NfRelElem{T}})(a::NfRelElem{NfRelElem{T}}) where T
  if base_ring(R)==base_ring(parent(a))
    return a.data
  end
  error("wrong ring")
end

function factor(f::Generic.Poly{NfRelElem{T}}) where T
  K = base_ring(f)
  Ka, rel_abs, _ = absolute_field(K)
  function map_poly(P::Ring, mp::Map, f::Generic.Poly)
    return P([mp(coeff(f, i)) for i=0:degree(f)])
  end

  fa = map_poly(PolynomialRing(Ka, "T", cached=false)[1], rel_abs, f)
  lf = factor(fa)
  res = Fac(map_poly(parent(f), pseudo_inv(rel_abs), lf.unit), Dict(map_poly(parent(f), pseudo_inv(rel_abs), k)=>v for (k,v) = lf.fac))

  return res
end

function isirreducible(f::Generic.Poly{NfRelElem{T}}) where T
  lf = factor(f)
  return sum(values(lf.fac)) == 1
end

function roots(f::Generic.Poly{NfRelElem{T}}) where T
  lf = factor(f)
  @assert degree(lf.unit) == 0
  scale = inv(coeff(lf.unit, 0))
  return [-constant_coefficient(x)*scale for x = keys(lf.fac) if degree(x)==1]
end

################################################################################
#
#  issubfield and isisomorphic
#
################################################################################

@doc Markdown.doc"""
    issubfield(K::NfRel, L::NfRel) -> Bool, NfRelToNfRelMor

Returns "true" and an injection from $K$ to $L$ if $K$ is a subfield of $L$.
Otherwise the function returns "false" and a morphism mapping everything to
$0$.
"""
function issubfield(K::NfRel, L::NfRel)
  @assert base_ring(K) == base_ring(L)
  f = K.pol
  g = L.pol
  if mod(degree(g), degree(f)) != 0
    return false, hom(K, L, zero(L), check = false)
  end
  Lx, x = PolynomialRing(L, "x", cached = false)
  fL = Lx()
  for i = 0:degree(f)
    setcoeff!(fL, i, L(coeff(f, i)))
  end
  fac = factor(fL)
  for (a, b) in fac
    if degree(a) == 1
      c1 = coeff(a, 0)
      c2 = coeff(a, 1)
      h = parent(K.pol)(-c1*inv(c2))
      return true, hom(K, L, h(gen(L)), check = false)
    end
  end
  return false, hom(K, L, zero(L), check = false)
end

@doc Markdown.doc"""
    isisomorphic(K::NfRel, L::NfRel) -> Bool, NfRelToNfRelMor

Returns "true" and an isomorphism from $K$ to $L$ if $K$ and $L$ are isomorphic.
Otherwise the function returns "false" and a morphism mapping everything to $0$.
"""
function isisomorphic(K::NfRel, L::NfRel)
  @assert base_ring(K) == base_ring(L)
  f = K.pol
  g = L.pol
  if degree(f) != degree(g)
    return false, hom(K, L, zero(L), check = false)
  end
  return issubfield(K, L)
end

@doc Markdown.doc"""
    discriminant(K::AnticNumberField) -> fmpq
    discriminant(K::NfRel) ->
The discriminant of the defining polynomial of $K$ {\bf not} the discriminant
of the maximal order.
"""
function Nemo.discriminant(K::AnticNumberField)
  return discriminant(K.pol)
end

function Nemo.discriminant(K::NfRel)
  d = discriminant(K.pol)
  return d
end

@doc Markdown.doc"""
    discriminant(K::AnticNumberField, FlintQQ) -> fmpq
    discriminant(K::NfRel, FlintQQ) ->
The absolute discriminant of the defining polynomial of $K$ {\bf not} the discriminant
of the maximal order. Ie the norm of the discriminant time the power of the discriminant
of the base field.
"""
function Nemo.discriminant(K::AnticNumberField, ::FlintRationalField)
  return discriminant(K)
end

function Nemo.discriminant(K::NfRel, ::FlintRationalField)
  d = norm(discriminant(K)) * discriminant(base_field(K))^degree(K)
  return d
end


text:
module ParamSetmod

"""
Liste les paramètres utilisés dans la résolution du MPCC
"""

type ParamSet
 #paramètres pour la résolution du MPCC
 precmpcc::Float64

 #paramètres algorithmiques pour la relaxation
 # le t_bar
 tb::Function #dépend des trois paramètres (r,s,t)
 # l'oracle de choix de epsilon
 prec_oracle::Function #dépend des trois paramètres (r,s,t)
 # la mise à jour des paramètres
 rho_restart::Function #
 paramin::Float64 #valeur minimal pour les paramères (r,s,t)

 #paramètres algorithmiques pour l'activation de contraintes
 ite_max_alas::Int64 #entier >=0
 ite_max_viol::Int64 #entier >=0
 rho_init::Vector #nombre >= 0 ou un vecteur ?
 rho_update::Float64 #nombre >=1
 rho_max::Float64
 goal_viol::Float64 #nombre entre (0,1)

 #paramètres pour la minimisation sans contrainte
 ite_max_armijo::Int64 #nombre maximum d'itération pour la recherche linéaire >= 0
 tau_armijo::Float64 #paramètre pour critère d'Armijo doit être entre (0,0.5)
 armijo_update::Float64 #step=step_0*(1/2)^m
 ite_max_wolfe::Int64
 tau_wolfe::Float64 #entre (tau_armijo,1)
 wolfe_update::Float64

 #paramètres de sortie:
 verbose::Int64
end

#Constructeur par défaut
#nbc est le nombre de contraintes du problème
function ParamSet(nbc::Int64)
 precmpcc=1e-3

 tb=(r,s,t)->-r
 prec_oracle=(r,s,t,prec)->max(r,s,t,prec)
 rho_restart(r,s,t,prec,rho)=rho
 paramin=sqrt(eps(Float64))
 #paramin=1.0

 ite_max_alas=100
 ite_max_viol=10
 rho_init=1*ones(nbc)
 rho_update=1.5 #2.0 bien pour Newton, trop grand sinon.
 rho_max=1/precmpcc #le grand max est 1/eps(Float64)
 goal_viol=0.5

 ite_max_armijo=300
 tau_armijo=0.1 #0.4
 armijo_update=0.9 #0.9
 ite_max_wolfe=10
 tau_wolfe=0.6
 wolfe_update=5.0

 verbose=0 #0 quiet, 1 relaxation, 2 relaxation+activation, 3 relaxation+activation+linesearch

 return ParamSet(precmpcc,tb,prec_oracle,rho_restart,paramin,ite_max_alas,ite_max_viol,rho_init,rho_update,rho_max,goal_viol,ite_max_armijo,tau_armijo,armijo_update,ite_max_wolfe,tau_wolfe,wolfe_update,verbose)
end

#end of module
end

ext: jl


text:
#!/usr/bin/env julia

using Eigenbroetler

function starry(x::Number, y::Number)
    r = sqrt(x * x + y * y)
    ϕ = atan(y, x)
    return sin(r * (3 + cos(7 * ϕ)) / 4)
end

eb = Eigenbrot(starry, 512, 512)
save("starry.png", eb, ImageSetting(RealPart, LinearScale))
save("starry_dft.png", fft(eb), ImageSetting(Magn, LogScale))


text:
using Documenter
using IscaTools

makedocs(
    modules  = [IscaTools],
    doctest  = false,
    format   = Documenter.HTML(
        collapselevel = 1,
        prettyurls    = false
    ),
    authors  = "Nathanael Wong",
    sitename = "IscaTools.jl",
    pages    = [
        "Home"     => "index.md",
        # "Tutorials" => [
        #     "Initilization"        => "tutorials/initialize.md",
        #     "Preliminary Analysis" => "tutorials/analysis.md"
        # ]
    ]
)

deploydocs(
    repo = "github.com/natgeo-wong/IscaTools.jl.git",
)


text:
@testset "struct: VectorPEPOTensor" begin
    S1 = Z2Space(0 => 1, 1 => 1)
    A1 = PEPOTensor(TensorMap(randn, ComplexF64, S1 * S1 * S1, S1 * S1 * S1))
    @test A1 ≈ PEPOKit.unvec(vec(A1))
    S2 = U1Space(-1 => 1, 0 => 3, 1 => 1)
    A2 = PEPOTensor(TensorMap(randn, ComplexF64, S2 * S2 * S2, S2 * S2 * S2))
    @test A2 ≈ PEPOKit.unvec(vec(A2))
    S3 = U1Space(-1 => 1, 1 => 1)
    A3 = PEPOTensor(TensorMap(randn, ComplexF64, S3 * S3 * S2, S3 * S3 * S2))
    @test A3 ≈ PEPOKit.unvec(vec(A3))
end

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
using NearestNeighbors
using Distances
using LinearAlgebra
data = rand(3, 10^4);
point = rand(3);
#println("Our point")
#println(point)
tree = BallTree(data);

t = @elapsed minimum.(data)
t = @elapsed minimum.(data)
println("Reference time: ",t)
@elapsed tree = BallTree(data);
t = @elapsed tree = BallTree(data);
println("BallTree time: ",t)

@elapsed kdtree = KDTree(data);

k = 10;
t = @elapsed idxs, dists = knn(kdtree, point, k, true);
t = @elapsed idxs, dists = knn(kdtree, point, k, true);
#println("time:",t)

# display(idxs)
# println("")
# display(dists)
# println("")

@elapsed idxs, dists = knn(tree, point, k, true)
t = @elapsed idxs, dists = knn(tree, point, k, true)
#println("time:",t)

# display(idxs)
# println("")
# display(dists)
# println("")

r = 0.02;
@elapsed idxs = inrange(tree, point, r)
t = @elapsed idxs = inrange(tree, point, r)
# println("time:",t)
# display(data[:,idxs])
#println("**************************************************************************************")
A = randn(3,3);
A = A'*A + 0.1*Matrix(1.0I,3,3);
t = @elapsed tree = BallTree(data,Mahalanobis{Float64}(A));
t = @elapsed tree = BallTree(data,Mahalanobis{Float64}(A));
println("time:",t)

t = @elapsed idxs, dists = knn(tree, point, k, true)

t = @elapsed idxs, dists = knn(tree, point, k, true)
# println("time:",t)

# display(idxs)
# println("")
# display(dists)
# println("")

r = 0.02;
@elapsed idxs = inrange(tree, point, r)
t = @elapsed idxs = inrange(tree, point, r)
# println("time:",t)
# display(data[:,idxs])


text:
using Test
using CVChannel

"""
This script verifies that the channels formed from axisymmetric states are multiplicative
for qubits and qutrits.
"""

println("Verifying qubit multiplicativity of axisymmetric channels")
@testset "qubit multiplicativity of axisymmetric channels" begin
    y_step = 0.1
    x_step = 0.1

    d = 2
    y_bounds = CVChannel._axisymmetric_y_bounds(d)

    it = 0
    y_range = y_bounds[1]:y_step:y_bounds[2]
    for y in y_range
        it = it + 1
        println("verifying y-slice $it of $(length(y_range))")
        x_constraints = CVChannel._axisymmetric_x_constraints(d,y)
        for x in x_constraints[1]:x_step:x_constraints[2]
            @testset "(x,y) = ($x,$y)" begin
                ρ_axi = axisymmetricState(d,x,y)

                J_N = d*ρ_axi
                (opt_cv_N, opt_σAB_N) = pptCVDual(J_N, d, d)

                J_NN = permuteSubsystems(kron(J_N,J_N), [1,3,2,4], [d,d,d,d])
                (opt_cv_NN, opt_σAB_NN) = pptCVPrimal(J_NN, d^2, d^2)

                multiplicativity = opt_cv_NN - opt_cv_N^2

                @test multiplicativity ≈ 0 atol=2e-5
            end
        end
    end
end

println("Verifying qutrit multiplicativity of axisymmetric channels")
@time @testset "qutrit multiplicativity of axisymmetric channels" begin
    y_step = 0.1
    x_step = 0.1

    d = 3

    y_bounds = CVChannel._axisymmetric_y_bounds(d)

    y_range = y_bounds[1]:y_step:y_bounds[2]
    it = 0
    for y in y_range
        it = it + 1
        println("verifying y-slice $it of $(length(y_range))")
        x_constraints = CVChannel._axisymmetric_x_constraints(d,y)
        for x in x_constraints[1]:x_step:x_constraints[2]
            @testset "(x,y) = ($x,$y)" begin
                ρ_axi = axisymmetricState(d,x,y)

                J_N = d*ρ_axi
                (opt_cv_N, opt_σAB_N) = pptCVDual(J_N, d, d)

                J_NN = permuteSubsystems(kron(J_N,J_N), [1,3,2,4], [d,d,d,d])
                (opt_cv_NN, opt_σAB_NN) = pptCVPrimal(J_NN, d^2, d^2)

                multiplicativity = opt_cv_NN - opt_cv_N^2

                @test multiplicativity ≈ 0 atol=2e-5
            end
        end
    end
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
## Define Terms type that manages formula parsing and extension.

abstract AbstractTerm

type Term{H} <: AbstractTerm
    children::Vector{Union{Term,Symbol}}

    Term() = new(Term[])
    Term(children::Vector) = add_children!(new(Term[]), children)
    Term(child::Symbol) = new([child])
end

typealias InterceptTerm Union{Term{0}, Term{-1}, Term{1}}

## equality of Terms
import Base.==
=={G,H}(::Term{G}, ::Term{H}) = false
=={H}(a::Term{H}, b::Term{H}) = a.children == b.children

## display of terms
function Base.show{H}(io::IO, t::Term{H})
    print(io, string(H))
    if length(t.children) > 0
        print(io, "(")
        print(io, join(map(string, t.children), ", "))
        print(io, ")")
    end
end
Base.show(io::IO, t::Term{:eval}) = print(io, string(t.children[1]))
## show ranef term:
Base.show(io::IO, t::Term{:|}) = print(io, "(", t.children[1], " | ", t.children[2], ")")

## Converting to Term:

## Symbols are converted to :eval terms (leaf nodes)
Base.convert(::Type{Term}, s::Symbol) = Term{:eval}(s)

## Integers to intercept terms
function Term(i::Integer)
    i == 0 || i == -1 || i == 1 || error("Can't construct term from Integer $i")
    Term{i}()
end

## no-op constructor
Term{H}(t::Term{H}) = t

## convert from one head type to another
(::Type{Term{H}}){H,J}(t::Term{J}) = add_children!(Term{H}(), t, [])

## Expressions are recursively converted to Terms, depth-first, and then added
## as children.  Specific `add_children!` methods handle special cases like
## associative and distributive operators.
function Base.convert(::Type{Term}, ex::Expr)
    ex.head == :call || error("non-call expression detected: '$(ex.head)'")
    add_children!(Term{ex.args[1]}(), [Term(a) for a in ex.args[2:end]])
end

## Adding children to a Term

## General strategy: add one at a time to allow for dispatching on special
## cases, but also keep track of the rest of the children being added because at
## least the distributive rule requires that context.
add_children!(t::Term, new_children::Vector) =
    isempty(new_children) ?
    t :
    add_children!(t::Term, new_children[1], new_children[2:end])

function add_children!(t::Term, c::Term, others::Vector)
    push!(t.children, c)
    add_children!(t, others)
end

## special cases:
## Associative rule
add_children!(t::Term{:+}, new_child::Term{:+}, others::Vector) =
    add_children!(t, cat(1, new_child.children, others))
add_children!(t::Term{:&}, new_child::Term{:&}, others::Vector) =
    add_children!(t, cat(1, new_child.children, others))

## Distributive property
## &(a..., +(b...), c...) -> +(&(a..., b_i, c...)_i...)
add_children!(t::Term{:&}, new_child::Term{:+}, others::Vector) =
    Term{:+}([add_children!(deepcopy(t), c, others) for c in new_child.children])

## Expansion of a*b -> a + b + a&b
expand_star(a::Term,b::Term) = Term{:+}([a, b, Term{:&}([a,b])])
add_children!(t::Term, new_child::Term{:*}, others::Vector) =
    add_children!(Term{:+}(reduce(expand_star, new_child.children)), others)

## Handle - for intercept term -1, and throw error otherwise
add_children!(t::Term{:-}, children::Vector) =
    isa(children[2], Term{1}) ?
    Term{:+}([children[1], Term{-1}()]) :
    error("invalid subtraction of $(children[2]); subtraction only supported for -1")

## sorting term by the degree of its children: order is 1 for everything except
## interaction Term{:&} where order is number of children
degree(t::Term{:&}) = length(t.children)
degree(::Term) = 1
degree(::InterceptTerm) = 0

function Base.sort!(t::Term)
    sort!(t.children, by=degree)
    return t
end

## extract evaluation terms: children of Term{:+} and Term{:&}, nothing for
## ranef Term{:|} and intercept terms, and Term itself for everything else.
evt(t::Term) = Term[t]
evt(t::Term{:eval}) = t.children
evt(t::Term{:&}) = mapreduce(evt, vcat, t.children)
evt(t::Term{:+}) = mapreduce(evt, vcat, t.children)
evt(t::Term{:|}) = Term[]
evt(t::InterceptTerm) = Term[]

## whether a Term is for fixed effects or not
isfe(t::Term{:|}) = false
isfe(t::Term) = true

################################################################################
## Constructing a DataFrames.Terms object

type Terms
    terms::Vector
    eterms::Vector        # evaluation terms
    factors::Matrix{Bool} # maps terms to evaluation terms
    ## An eterms x terms matrix which is true for terms that need to be "promoted"
    ## to full rank in constructing a model matrx
    is_non_redundant::Matrix{Bool}
# order can probably be dropped.  It is vec(sum(factors, 1))
    order::Vector{Int}    # orders of rhs terms
    response::Bool        # indicator of a response, which is eterms[1] if present
    intercept::Bool       # is there an intercept column in the model matrix?
end

Base.:(==)(t1::Terms, t2::Terms) = all(getfield(t1, f)==getfield(t2, f) for f in fieldnames(t1))

function Terms(f::Formula)
    ## start by raising everything on the right-hand side by converting
    rhs = sort!(Term{:+}(Term(f.rhs)))
    terms = rhs.children

    ## detect intercept
    is_intercept = [isa(t, InterceptTerm) for t in terms]
    hasintercept = mapreduce(t -> isa(t, Term{1}),
                             &,
                             true, # default is to have intercept
                             terms[is_intercept])

    terms = terms[!is_intercept]
    degrees = map(degree, terms)

    evalterms = map(evt, terms)

    haslhs = f.lhs != nothing
    if haslhs
        lhs = Term(f.lhs)
        unshift!(evalterms, evt(lhs))
        unshift!(degrees, degree(lhs))
    end

    evalterm_sets = [Set(x) for x in evalterms]
    evalterms = unique(vcat(evalterms...))

    factors = Int8[t in s for t in evalterms, s in evalterm_sets]
    non_redundants = fill(false, size(factors)) # initialize to false

    Terms(terms, evalterms, factors, non_redundants, degrees, haslhs, hasintercept)

end


text:
# @recipe function f(e::Event;coloring=[], labeling=[],show_pulses=true)
#     if show_pulses==true
#         n_wps = size(e.weighting_potentials,1)
#         width = 800
#         length = width+n_wps*width/2
#         size  -->  (width,length)
#         myheights = [width/length]
#         for i in 1:n_wps
#             push!(myheights,(width/2)/length)
#         end

#         layout  -->  (n_wps+1,1) #grid(n_wps+1,1,heights=myheights)
#         @series begin
#             subplot := 1
#             coloring --> coloring
#             labeling --> labeling
#             e.detector
#         end
#         for itr in 1:e.n_sites
#             @series begin
#                 itr == 1 ? showlabel --> true : showlabel --> false
#                 subplot := 1
#                 myscale := 1/e.detector.geometry_unit_factor
#                 e.trajectories_e[itr]
#             end
#             @series begin
#                 itr == 1 ? showlabel --> true : showlabel --> false
#                 subplot := 1
#                 myscale := 1/e.detector.geometry_unit_factor
#                 e.trajectories_h[itr]
#             end
#         end
#         for iwp in 1:n_wps
#             @series begin
#                 subplot := 1+iwp
#                 e.pulses[iwp]
#             end
#         end
#     else
#         width,length = 800, 800
#         size  -->  (width,length)
#         @series begin
#             subplot := 1
#             coloring --> coloring
#             labeling --> labeling
#             e.detector
#         end

#         for itr in 1:e.n_sites
#             @series begin
#                 itr == 1 ? showlabel --> true : showlabel --> false
#                 subplot := 1
#                 myscale := 1/e.detector.geometry_unit_factor
#                 e.trajectories_e[itr]
#             end
#             @series begin
#                 itr == 1 ? showlabel --> true : showlabel --> false
#                 subplot := 1
#                 myscale := 1/e.detector.geometry_unit_factor
#                 e.trajectories_h[itr]
#             end
#         end
#     end
# end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RandomThinning(p)

Random thining with retention probability `p`.
"""
struct RandomThinning{P<:Union{Real,Function}} <: AbstractThinning
  p::P
end

# -----------------------
# thinning point process
# -----------------------
thin(p::PoissonProcess{<:Real}, t::RandomThinning{<:Real}) =
  PoissonProcess(t.p * p.λ)

thin(p::PoissonProcess{<:Function}, t::RandomThinning{<:Function}) =
  PoissonProcess(x -> t.p(x) * p.λ(x))

thin(p::PoissonProcess{<:Real}, t::RandomThinning{<:Function}) =
  PoissonProcess(x -> t.p(x) * p.λ)

thin(p::PoissonProcess{<:Function}, t::RandomThinning{<:Real}) =
  PoissonProcess(x -> t.p * p.λ(x))

# -----------------------
# thinning point pattern
# -----------------------
function thin(pp::PointSet, t::RandomThinning{<:Real})
  draws = rand(Bernoulli(t.p), nelements(pp))
  inds  = findall(isequal(1), draws)
  PointSet(coordinates(pp, inds))
end

function thin(pp::PointSet, t::RandomThinning{<:Function})
  inds = Vector{Int}()
  for j in 1:nelements(pp)
    x = coordinates(pp, j)
    if rand(Bernoulli(t.p(x)))
      push!(inds, j)
    end
  end
  PointSet(coordinates(pp, inds))
end


text:
module MolecularBoxes

export MolecularBox, Box
export isperiodic, separation, center_of_mass

using StaticArrays

"""
MolecularBox has three type parameters:
    * V: An immutable vector such as those implemented in the StaticArrays
        package
    * N: The number of dimensions
    * P: A tuple of booleans indicating which dimensions are periodic

Concrete types inheriting from MolecularBox
must implement 2 fields: vectors and lengths
 si(::Box)
"""
abstract type MolecularBox{V,N,P} end

include("center_of_mass.jl")

@inline isperiodic(b::MolecularBox{V,N,P}) where {V,N,P} = P
@inline isperiodic(b::MolecularBox{V,N,P}, dim) where {V,N,P} = P[dim]

@inline Base.eltype(::MolecularBox{V}) where V = V

struct Box{V,N,P} <: MolecularBox{V,N,P}
    vectors::NTuple{N,V}
    lengths::V
    function Box{V,N,P}(vectors::NTuple{N,V}) where {N,V,P}
        if !isa(N,Integer)
            error("N parameter of a MolecularBox must be an Integer")
        end
        if !isa(P,NTuple{N,Bool})
            error(string("P=$P parameter of a `MolecularBox` must be a tuple ",
                         "of `Bool`s, with one `Bool` for each dimension"))
        end
        for i in 1:N
            if vectors[i][i] <=0
                error("Diagonal box matrix elements must be >= 0")
            end
            for j in 1:i-1
                if vectors[j][i] != 0
                    error("box vector element [$(j)][$(i)] must be zero")
                end
            end
            for j in i+1:N
                if vectors[j][i] != 0
                    warn("Triclinic box support not fully tested ($j,$i)")
                end
            end
        end
        lengths = convert(V, collect(vectors[i][i] for i in 1:N) )
        new{V,N,P}(vectors, lengths)
    end
end

function Box(
    lengths::V;
    periodic=( (true for x in 1:length(lengths))..., )
) where V
    N = length(lengths)
    vectors = map( (i,L)-> begin
        v = zeros(eltype(V), N)
        v[i] = L
        convert(V,v)
    end, 1:N, lengths)
    Box{V,N,(periodic...,)}((vectors...,))
end

@inline _separation(
    x1,
    x2,
    length,
    periodic::Bool,
) = _separation(x1, x2, length, Val{periodic})

@inline function _separation(
    x1,
    x2,
    length,
    periodic::Type{Val{P}},
) where {P}
    r = x1-x2
    if P
        hlen = 0.5length
        while r >= hlen
            r-=length
        end
        while r < -hlen
            r+=length
        end
        # faster but unsafe version:
#       r + ifelse(
#           r >= hlen,
#           -length,
#           ifelse(r < -hlen, length, zero(length)),
#       )
        r
    else
        r
    end
end

@inline function separation(
    x1,
    x2,
    box::MolecularBox{V,N,P},
) where {V,N,P}
    _separation.(x1, x2, box.lengths, SVector(P))
end

separation

end #module


text:
# Test of Depth-first visit

using Graphs
using Base.Test

type GraphTest
    graph_edges::Array{(Int,Int),1}
    dfs_path::Array{Int,1}
    is_cyclic::Bool
    topo_sort::Array{Int,1}

    GraphTest(gedges, dfspath, iscyclic, toposort) =
        new(gedges, dfspath, iscyclic, toposort)
end

dir_acyclic = GraphTest(
    [(1,2), (1,3), (1,6), (2,4), (2,5), (3,5), (3,6)],
    [1, 2, 4, 5, 3, 6],
    false,
    [1,3,6,2,5,4])
undir_acyclic = GraphTest(
    [(1,2), (1,3), (1,6), (2,4), (2,5)],
    [],
    false,
    [1,6,3,2,5,4])
cyclic  = GraphTest(
    [(1,2), (1,3), (1,6), (2,4), (2,5), (3,5), (3,6), (5,1)],
    [1, 2, 4, 5, 3, 6],
    true,
    [])

testsets = [
    (true, [dir_acyclic, cyclic]),
    (false, [undir_acyclic, cyclic])]

for tset in testsets
    is_dir, graphtests = tset

    for gtest in graphtests

        g = simple_inclist(6, is_directed = is_dir)
        map((edg) -> add_edge!(g, edg[1], edg[2]), gtest.graph_edges)

        gEx = graph(ExVertex[], ExEdge{ExVertex}[], is_directed = is_dir)
        map((x) -> add_vertex!(gEx, "edge:" * string(x)), 1:6)
        V = vertices(gEx)
        map((edg) -> add_edge!(gEx, V[edg[1]], V[edg[2]]), gtest.graph_edges)

        # DFS traversal
        if !isempty(gtest.dfs_path)
            vs1 = visited_vertices(g, DepthFirst(), 1)
            @assert vs1 == gtest.dfs_path

            vs2 = visited_vertices(gEx, DepthFirst(), V[1])
            @assert vs2 == collect(map((x) -> gEx.vertices[x], gtest.dfs_path))
        end

        # Cyclic test
        @assert test_cyclic_by_dfs(g) == gtest.is_cyclic
        @assert test_cyclic_by_dfs(gEx) == gtest.is_cyclic

        # Topological sort
        if gtest.is_cyclic
            @test_throws ArgumentError topological_sort_by_dfs(g)  # g2 contains a loop
            @test_throws ArgumentError topological_sort_by_dfs(gEx)  # g2 contains a loop

        elseif !isempty(gtest.topo_sort)
            ts = topological_sort_by_dfs(g)
            @assert ts == gtest.topo_sort

            ts = topological_sort_by_dfs(gEx)
            @assert [vertex_index(e,gEx) for e in ts] == gtest.topo_sort
        end
    end
end


text:
"""
    kernelmatrix!(K::AbstractMatrix, κ::Kernel, X; obsdim::Integer = 2)
    kernelmatrix!(K::AbstractMatrix, κ::Kernel, X, Y; obsdim::Integer = 2)

In-place version of [`kernelmatrix`](@ref) where pre-allocated matrix `K` will be
overwritten with the kernel matrix.
"""
kernelmatrix!

"""
    kernelmatrix(κ::Kernel, X; obsdim::Int = 2)
    kernelmatrix(κ::Kernel, X, Y; obsdim::Int = 2)

Calculate the kernel matrix of `X` (and `Y`) with respect to kernel `κ`.
`obsdim = 1` means the matrix `X` (and `Y`) has size #samples x #dimension
`obsdim = 2` means the matrix `X` (and `Y`) has size #dimension x #samples
"""
kernelmatrix

"""
    kerneldiagmatrix!(K::AbstractVector, κ::Kernel, X; obsdim::Int = 2)
    kerneldiagmatrix!(K::AbstractVector, κ::Kernel, X, Y; obsdim::Int = 2)

In place version of [`kerneldiagmatrix`](@ref)
"""
kerneldiagmatrix!

"""
    kerneldiagmatrix(κ::Kernel, X; obsdim::Int = 2)

Calculate the diagonal matrix of `X` with respect to kernel `κ`
`obsdim = 1` means the matrix `X` has size #samples x #dimension
`obsdim = 2` means the matrix `X` has size #dimension x #samples

    kerneldiagmatrix(κ::Kernel, X, Y; obsdim::Int = 2)

Calculate the diagonal of `kernelmatrix(κ, X, Y; obsdim)` efficiently. Requires that `X` and
`Y` are the same length.
"""
kerneldiagmatrix

#
# Kernel implementations. Generic fallbacks that depend only on kernel evaluation.
#

kernelmatrix!(K::AbstractMatrix, κ::Kernel, x::AbstractVector) = kernelmatrix!(K, κ, x, x)

function kernelmatrix!(K::AbstractMatrix, κ::Kernel, x::AbstractVector, y::AbstractVector)
    validate_inplace_dims(K, x, y)
    K .= κ.(x, permutedims(y))
    return K
end

kernelmatrix(κ::Kernel, x::AbstractVector) = kernelmatrix(κ, x, x)

function kernelmatrix(κ::Kernel, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    return κ.(x, permutedims(y))
end

function kerneldiagmatrix!(K::AbstractVector, κ::Kernel, x::AbstractVector)
    validate_inplace_dims(K, x)
    return map!(x -> κ(x, x), K, x)
end

function kerneldiagmatrix!(
    K::AbstractVector, κ::Kernel, x::AbstractVector, y::AbstractVector,
)
    return map!(κ, x, y)
end

kerneldiagmatrix(κ::Kernel, x::AbstractVector) = map(x -> κ(x, x), x)

kerneldiagmatrix(κ::Kernel, x::AbstractVector, y::AbstractVector) = map(κ, x, y)

#
# SimpleKernel optimisations.
#

function kernelmatrix!(K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector)
    validate_inplace_dims(K, x)
    pairwise!(K, metric(κ), x)
    return map!(d -> kappa(κ, d), K, K)
end

function kernelmatrix!(
    K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector, y::AbstractVector,
)
    validate_inplace_dims(K, x, y)
    pairwise!(K, metric(κ), x, y)
    return map!(d -> kappa(κ, d), K, K)
end

function kernelmatrix(κ::SimpleKernel, x::AbstractVector)
    return map(d -> kappa(κ, d), pairwise(metric(κ), x))
end

function kernelmatrix(κ::SimpleKernel, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    return map(d -> kappa(κ, d), pairwise(metric(κ), x, y))
end

#
# Wrapper methods for AbstractMatrix inputs to maintain obsdim interface.
#

const defaultobs = 2

function kernelmatrix!(
    K::AbstractMatrix, κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs,
)
    return kernelmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim))
end

function kernelmatrix!(
    K::AbstractMatrix, κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix;
    obsdim::Int=defaultobs,
)
    return kernelmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim), vec_of_vecs(Y; obsdim=obsdim))
end

function kernelmatrix(κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs)
    return kernelmatrix(κ, vec_of_vecs(X; obsdim=obsdim))
end

function kernelmatrix(κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim=defaultobs)
    return kernelmatrix(κ, vec_of_vecs(X; obsdim=obsdim), vec_of_vecs(Y; obsdim=obsdim))
end

function kerneldiagmatrix!(
    K::AbstractVector, κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs
)
    return kerneldiagmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim))
end

function kerneldiagmatrix!(
    K::AbstractVector, κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix;
    obsdim::Int = defaultobs,
)
    return kerneldiagmatrix!(
        K, κ, vec_of_vecs(X; obsdim=obsdim), vec_of_vecs(Y; obsdim=obsdim),
    )
end

function kerneldiagmatrix(κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs)
    return kerneldiagmatrix(κ, vec_of_vecs(X; obsdim=obsdim))
end

function kerneldiagmatrix(
    κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int=defaultobs,
)
    return kerneldiagmatrix(κ, vec_of_vecs(X; obsdim=obsdim), vec_of_vecs(Y; obsdim=obsdim))
end


text:
#=
# 4. Multiple stationary bodies
Adding multiple bodies to a problem is easy, using the concepts of a `BodyList`
and `RigidTransformList`.
=#
using ViscousFlow
#-
using Plots

#=
In this example, we will set up a problem with three cylinders arranged in a
formation in a free stream.
=#
my_params = Dict()
my_params["Re"] = 200
my_params["freestream speed"] = 1.0
my_params["freestream angle"] = 0.0

#=
Set up the domain and surface point spacing
=#
xlim = (-2.0,4.0)
ylim = (-2.0,2.0)
my_params["grid Re"] = 4.0
g = setup_grid(xlim,ylim,my_params)
Δs = surface_point_spacing(g,my_params)

#=
### Set up bodies
We start by initializing a `BodyList` and an associated `RigidTransformList`.
Each member of the `RigidTransformList` will be used to place the respective body in
the correct position and orientation.
=#
bl = BodyList()
tl = RigidTransformList()

# Place the first body at (-1,0)
push!(bl,Circle(0.5,Δs))
push!(tl,RigidTransform((-1,0),0));

# Place the second body at (1,-1)
push!(bl,Circle(0.5,Δs))
push!(tl,RigidTransform((1,-1),0));

# and place the third body at (1,1)
push!(bl,Circle(0.5,Δs))
push!(tl,RigidTransform((1,1),0))

# Perform the actual transformation. Note that this operation works `in-place`:
tl(bl)

# #### Plot the initial configuration of the bodies
# Just to check they are in the right places
plot(bl,xlim=xlim,ylim=ylim)

#=
### Construct the system structure
We construct the system with the same syntax as for a single body:
=#
sys = viscousflow_system(g,bl,phys_params=my_params);

#-
u0 = init_sol(sys)
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)
#=
### Solve
Here, we run it for a little while, just to demonstrate:
=#
@time step!(integrator,2.0)
#=
### Examine
Let's inspect the results
=#
sol = integrator.sol
plt = plot(layout = (1,3), size = (800, 200), legend=:false)
tsnap = 0.0:1.0:2.0
for (i, t) in enumerate(tsnap)
    plot!(plt[i],vorticity(sol,sys,t),sys,clim=(-10,10),levels=range(-10,10,length=30), color = :RdBu)
end
plt

# Now we will examine the force on each body
fx1, fy1 = force(sol,sys,1)
fx2, fy2 = force(sol,sys,2)
fx3, fy3 = force(sol,sys,3);
#-
plt = plot(layout = (2,1), size = (600, 400))
plot!(plt[1],sol.t,2*fx1,xlim=(0,Inf),ylim=(0,4),xlabel="Convective time",ylabel="\$C_D\$",label="Lead body",title="Drag force")
plot!(plt[2],sol.t,2*fy1,xlim=(0,Inf),ylim=(-2,2),xlabel="Convective time",ylabel="\$C_L\$",label="Lead body",title="Side force")
plot!(plt[1],sol.t,2*fx2,xlim=(0,Inf),ylim=(0,4),xlabel="Convective time",ylabel="\$C_D\$",label="Trailing body",title="Drag force")
plot!(plt[2],sol.t,2*fy2,xlim=(0,Inf),ylim=(-2,2),xlabel="Convective time",ylabel="\$C_L\$",label="Trailing body",title="Side force")
#-
println("Mean drag coefficient on lead body = ", GridUtilities.mean(2*fx1[3:end]))
#-
println("Mean drag coefficient on trailing body = ", GridUtilities.mean(2*fx2[3:end]))


text:
using FactCheck, Images, Colors, FixedPointNumbers

macro chk(a, b)
    :(@fact ($a == $b && typeof($a) == typeof($b)) --> true)
end

macro chk_approx(a, b)
    :(@fact (abs($a - $b) < 2*(eps($a)+eps($b)) && typeof($a) == typeof($b)) --> true)
end

facts("Map") do
    context("MapNone") do
        mapi = MapNone{Int}()
        @chk map(mapi, 7) 7
        @chk map(mapi, 0x07) 7
        a7 = Int[7]
        @chk map(mapi, [0x07]) a7
        @fact map(mapi, a7) --> exactly(a7)
        mapi = MapNone{RGB24}()
        g = UFixed8(0.1)
        @chk Images.map1(mapi, 0.1) g
        @chk map(mapi, 0.1) RGB24(g,g,g)
        @chk map(mapi, Gray(0.1)) RGB24(g,g,g)
        @chk map(mapi, g) RGB24(g,g,g)
        @chk map(mapi, true) RGB24(1,1,1)
        @chk map(mapi, false) RGB24(0x00,0x00,0x00)
        mapi = MapNone{RGB{Float32}}()
        g = 0.1f0
        @chk Images.map1(mapi, 0.1) g
        @chk map(mapi, 0.1) RGB(g,g,g)

        c = RGB(0.1,0.2,0.3)
        mapi = MapNone{HSV{Float64}}()
        @chk map(mapi, c) convert(HSV, c)

        # issue #200
        c = RGBA{UFixed8}(1,0.5,0.25,0.8)
        mapi = MapNone{Images.ColorTypes.BGRA{UFixed8}}()
        @chk map(mapi, c) convert(Images.ColorTypes.BGRA{UFixed8}, c)
    end

    context("BitShift") do
        mapi = BitShift{UInt8,7}()
        @chk map(mapi, 0xff) 0x01
        @chk map(mapi, 0xf0ff) 0xff
        @chk map(mapi, [0xff]) UInt8[0x01]
        mapi = BitShift{UFixed8,7}()
        @chk map(mapi, 0xffuf8) 0x01uf8
        @chk map(mapi, 0xf0ffuf16) 0xffuf8
        mapi = BitShift{ARGB32,4}()
        @chk map(mapi, 0xffuf8) ARGB32(0xff0f0f0f)
        mapi = BitShift{RGB24,2}()
        @chk map(mapi, Gray(0xffuf8)) RGB24(0x003f3f3f)
        mapi = BitShift{ARGB32,2}()
        @chk map(mapi, Gray(0xffuf8)) ARGB32(0xff3f3f3f)
        @chk map(mapi, GrayA{UFixed8}(Gray(0xffuf8),0x3fuf8)) ARGB32(0x0f3f3f3f)
        mapi = BitShift{RGB{UFixed8},2}()
        @chk map(mapi, Gray(0xffuf8)) RGB(0x3fuf8, 0x3fuf8, 0x3fuf8)
        mapi = BitShift{ARGB{UFixed8},2}()
        @chk map(mapi, Gray(0xffuf8)) ARGB{UFixed8}(0x3fuf8, 0x3fuf8, 0x3fuf8, 0xffuf8)
        @chk map(mapi, GrayA{UFixed8}(Gray(0xffuf8),0x3fuf8)) ARGB{UFixed8}(0x3fuf8, 0x3fuf8, 0x3fuf8, 0x0fuf8)
        mapi = BitShift{RGBA{UFixed8},2}()
        @chk map(mapi, Gray(0xffuf8)) RGBA{UFixed8}(0x3fuf8, 0x3fuf8, 0x3fuf8, 0xffuf8)
        @chk map(mapi, GrayA{UFixed8}(Gray(0xffuf8),0x3fuf8)) RGBA{UFixed8}(0x3fuf8, 0x3fuf8, 0x3fuf8, 0x0fuf8)
        mapi = BitShift(ARGB{UFixed8}, 8)
        @chk map(mapi, RGB{UFixed16}(1.0,0.8,0.6)) ARGB{UFixed8}(1.0,0.8,0.6,1.0)
        mapi = BitShift(RGBA{UFixed8}, 8)
        @chk map(mapi, RGB{UFixed16}(1.0,0.8,0.6)) RGBA{UFixed8}(1.0,0.8,0.6,1.0)
        # Issue, #269, IJulia issue #294
        bs = BitShift(Gray{UFixed8}, 8)
        v = Gray(ufixed16(0.8))
        @chk map(bs, v) Gray{UFixed8}(0.8)
    end

    context("Clamp") do
        mapi = ClampMin(Float32, 0.0)
        @chk map(mapi,  1.2) 1.2f0
        @chk map(mapi, -1.2) 0.0f0
        mapi = ClampMin(RGB24, 0.0f0)
        @chk map(mapi, RGB{Float32}(-5.3,0.4,0.8)) RGB24(0x000066cc)
        mapi = ClampMax(Float32, 1.0)
        @chk map(mapi,  1.2)  1.0f0
        @chk map(mapi, -1.2) -1.2f0
        mapi = ClampMax(RGB24, 1.0f0)
        @chk map(mapi, RGB{Float32}(0.2,1.3,0.8)) RGB24(0x0033ffcc)
        mapi = ClampMinMax(Float32, 0.0, 1.0)
        @chk map(mapi,  1.2) 1.0f0
        @chk map(mapi, -1.2) 0.0f0
        mapi = ClampMinMax(ARGB32, 0.0f0, 1.0f0)
        @chk map(mapi, RGB{Float32}(-0.2,1.3,0.8)) ARGB32(0xff00ffcc)
        @chk map(mapi, ARGB{Float32}(-0.2,1.3,0.8,0.6)) ARGB32(0x9900ffcc)
        mapi = Clamp(Float32)
        @chk map(mapi,  1.2) 1.0f0
        @chk map(mapi, -1.2) 0.0f0
        mapi = Clamp(UFixed12)
        @chk map(mapi, UFixed12(1.2)) one(UFixed12)
        mapi = Clamp(Gray{UFixed12})
        @chk map(mapi, Gray(UFixed12(1.2))) Gray(one(UFixed12))
        mapi = ClampMinMax(RGB24, 0.0, 1.0)
        @chk map(mapi, 1.2) RGB24(0x00ffffff)
        @chk map(mapi, 0.5) RGB24(0x00808080)
        @chk map(mapi, -.3) RGB24(0x00000000)
        mapi = ClampMinMax(RGB{UFixed8}, 0.0, 1.0)
        @chk map(mapi, 1.2) RGB{UFixed8}(1,1,1)
        @chk map(mapi, 0.5) RGB{UFixed8}(0.5,0.5,0.5)
        @chk map(mapi, -.3) RGB{UFixed8}(0,0,0)
        mapi = Clamp(RGB{UFixed8})
        @chk map(mapi, RGB(1.2,0.5,-.3)) RGB{UFixed8}(1,0.5,0)
        mapi = Clamp(ARGB{UFixed8})
        @chk map(mapi, ARGB{Float64}(1.2,0.5,-.3,0.2)) ARGB{UFixed8}(1.0,0.5,0.0,0.2)
        @chk map(mapi, RGBA{Float64}(1.2,0.5,-.3,0.2)) ARGB{UFixed8}(1.0,0.5,0.0,0.2)
        @chk map(mapi, 0.2) ARGB{UFixed8}(0.2,0.2,0.2,1.0)
        @chk map(mapi, GrayA{Float32}(Gray(0.2),1.2)) ARGB{UFixed8}(0.2,0.2,0.2,1.0)
        @chk map(mapi, GrayA{Float32}(Gray(-.4),0.8)) ARGB{UFixed8}(0.0,0.0,0.0,0.8)
        mapi = Clamp(RGBA{UFixed8})
        @chk map(mapi, ARGB{Float64}(1.2,0.5,-.3,0.2)) RGBA{UFixed8}(1.0,0.5,0.0,0.2)
        @chk map(mapi, RGBA{Float64}(1.2,0.5,-.3,0.2)) RGBA{UFixed8}(1.0,0.5,0.0,0.2)
        @chk map(mapi, 0.2) RGBA{UFixed8}(0.2,0.2,0.2,1.0)
        @chk map(mapi, GrayA{Float32}(Gray(0.2),1.2)) RGBA{UFixed8}(0.2,0.2,0.2,1.0)
        @chk map(mapi, GrayA{Float32}(Gray(-.4),0.8)) RGBA{UFixed8}(0.0,0.0,0.0,0.8)
        # Issue #253
        mapi = Clamp(BGRA{UFixed8})
        @chk map(mapi, RGBA{Float32}(1.2,0.5,-.3,0.2)) BGRA{UFixed8}(1.0,0.5,0.0,0.2)

        @chk clamp(RGB{Float32}(-0.2,0.5,1.8)) RGB{Float32}(0.0,0.5,1.0)
        @chk clamp(ARGB{Float64}(1.2,0.5,-.3,0.2)) ARGB{Float64}(1.0,0.5,0.0,0.2)
        @chk clamp(RGBA{Float64}(1.2,0.5,-.3,0.2)) RGBA{Float64}(1.0,0.5,0.0,0.2)
    end

    context("Issue #285") do
        a = [Gray(0xd0uf8)]
        a1 = 10*a
        mapi = mapinfo(Gray{UFixed8}, a1)
        @chk map(mapi, a1[1]) Gray(0xffuf8)
    end

    context("ScaleMinMax") do
        mapi = ScaleMinMax(UFixed8, 100, 1000)
        @chk map(mapi, 100) UFixed8(0.0)
        @chk map(mapi, 1000) UFixed8(1.0)
        @chk map(mapi, 10) UFixed8(0.0)
        @chk map(mapi, 2000) UFixed8(1.0)
        @chk map(mapi, 550) UFixed8(0.5)
        mapinew = ScaleMinMax(UFixed8, [100,500,1000])
        @fact mapinew --> mapi
        mapinew = ScaleMinMax(UFixed8, [0,500,2000], convert(UInt16, 100), convert(UInt16, 1000))
        @fact mapinew --> mapi
        mapi = ScaleMinMax(ARGB32, 100, 1000)
        @chk map(mapi, 100) ARGB32(0,0,0,1)
        @chk map(mapi, 550) ARGB32(0x80uf8,0x80uf8,0x80uf8,0xffuf8)
        @chk map(mapi,2000) ARGB32(1,1,1,1)
        mapi = ScaleMinMax(RGB{Float32}, 100, 1000)
        @chk map(mapi,  50) RGB(0.0f0, 0.0f0, 0.0f0)
        @chk map(mapi, 550) RGB{Float32}(0.5, 0.5, 0.5)
        @chk map(mapi,2000) RGB(1.0f0, 1.0f0, 1.0f0)
        A = Gray{UFixed8}[UFixed8(0.1), UFixed8(0.9)]
        @fact mapinfo(RGB24, A) --> MapNone{RGB24}()
        mapi = ScaleMinMax(RGB24, A, zero(Gray{UFixed8}), one(Gray{UFixed8}))
        @fact map(mapi, A) --> map(mapinfo(RGB24, A), A)
        mapi = ScaleMinMax(Float32, [Gray(one(UFixed8))], 0, 1) # issue #180
        @chk map(mapi, Gray(UFixed8(0.6))) 0.6f0
        @fact_throws ErrorException ScaleMinMax(Float32, 0, 0, 1.0) # issue #245
        A = [Gray{Float64}(0.2)]
        mapi = ScaleMinMax(RGB{UFixed8}, A, 0.0, 0.2)
        @fact map(mapi, A) --> [RGB{UFixed8}(1,1,1)]
        mapi = ScaleMinMax(Gray{U8}, Gray{U8}(0.2), Gray{U8}(0.4))
        @fact Gray{U8}(0.49) <= map(mapi, Gray{U8}(0.3)) <= Gray{U8}(0.5) --> true
        @fact Gray{U8}(0.49) <= map(mapi, 0.3) <= Gray{U8}(0.5) --> true
        mapi = ScaleMinMax(Gray{U8}, 0.2, 0.4)
        @fact Gray{U8}(0.49) <= map(mapi, Gray{U8}(0.3)) <= Gray{U8}(0.5) --> true
        @fact Gray{U8}(0.49) <= map(mapi, 0.3) <= Gray{U8}(0.5) --> true

        A = [-0.5 0.5; Inf NaN]
        @fact map(Clamp01NaN(A), A) --> [0.0 0.5; 1.0 0.0]
        B = colorim(repeat(reshape(A, (1,2,2)), outer=[3,1,1]))
        @fact map(Clamp01NaN(B), B) --> [RGB(0.0,0,0) RGB(0.5,0.5,0.5); RGB(1.0,1,1) RGB(0.0,0,0)]
    end

    context("ScaleSigned") do
        mapi = ScaleSigned(Float32, 1/5)
        @chk map(mapi, 7) 1.0f0
        @chk map(mapi, 5) 1.0f0
        @chk map(mapi, 3) convert(Float32, 3/5)
        @chk map(mapi, -3) convert(Float32, -3/5)
        @chk map(mapi, -6) -1.0f0
        mapi = ScaleSigned(RGB24, 1.0f0/10)
        @chk map(mapi, 12) RGB24(0x00ff00ff)
        @chk map(mapi, -10.0) RGB24(0x0000ff00)
        @chk map(mapi, 0) RGB24(0x00000000)
    end

    context("ScaleAutoMinMax") do
        mapi = ScaleAutoMinMax()
        A = [100,550,1000]
        @chk map(mapi, A) ufixed8([0.0,0.5,1.0])
        mapi = ScaleAutoMinMax(RGB24)
        @chk map(mapi, A) RGB24[0x00000000, 0x00808080, 0x00ffffff]

        # Issue #304
        A = rand(UInt16, 3, 2, 2)
        imgr = colorim(A)
        mi1 = ScaleAutoMinMax(RGB{UFixed16})
        res1 = raw(map(mi1, imgr))
        mi2 = ScaleAutoMinMax(UFixed16)
        res2 = raw(map(mi2, raw(imgr)))
        # @fact res1 --> res2
        # Note: this fails occassionally. Reproduce it with
        #    s = 1.1269798f0
        #    val = 0xdeb5
        #    UFixed16(s*UFixed16(val,0)) == UFixed16((s/typemax(UInt16))*val)
        @fact maxabs(convert(Array{Int32}, res1) - convert(Array{Int32}, res2)) --> less_than_or_equal(1)
    end

    context("Scaling and ssd") do
        img = Images.grayim(fill(typemax(UInt16), 3, 3))
        mapi = Images.mapinfo(UFixed8, img)
        img8 = map(mapi, img)
        @fact all(img8 .== typemax(UFixed8)) --> true
        A = 0
        mnA, mxA = 1.0, -1.0
        while mnA > 0 || mxA < 0
            A = randn(3,3)
            mnA, mxA = extrema(A)
        end
        offset = 30.0
        img = convert(Images.Image, A .+ offset)
        mapi = Images.ScaleMinMax(UFixed8, offset, offset+mxA, 1/mxA)
        imgs = map(mapi, img)
        @fact minimum(imgs) --> 0
        @fact maximum(imgs) --> 1
        @fact eltype(imgs) --> UFixed8
        imgs = Images.imadjustintensity(img)
        @fact_throws MethodError Images.imadjustintensity(img, [1])
        mnA = minimum(A)
        @fact Images.ssd(imgs, (A.-mnA)/(mxA-mnA)) --> less_than(eps(UFixed16))
        A = reshape(1:9, 3, 3)
        B = map(Images.ClampMin(Float32, 3), A)
        @fact (eltype(B) == Float32 && B == [3 4 7; 3 5 8; 3 6 9]) --> true
        B = map(Images.ClampMax(UInt8, 7), A)
        @fact (eltype(B) == UInt8 && B == [1 4 7; 2 5 7; 3 6 7]) --> true

        A = reinterpret(UFixed8, [convert(UInt8,1):convert(UInt8,24);], (3, 2, 4))
        img = reinterpret(RGB{UFixed8}, A, (2,4))
        @fact separate(img) --> permutedims(A, (2,3,1))
    end

    context("sc") do
        arr = zeros(4,4)
        arr[2,2] = 0.5
        @fact sc(arr)[2,2] --> 0xffuf8
        @fact sc(arr, 0.0, 0.75)[2,2] --> 0xaauf8
    end

    context("Color conversion") do
        gray = collect(linspace(0.0,1.0,5)) # a 1-dimensional image
        gray8 = round(UInt8, 255*gray)
        gray32 = UInt32[convert(UInt32, g)<<16 | convert(UInt32, g)<<8 | convert(UInt32, g) for g in gray8]
        imgray = Images.Image(gray, Dict{ASCIIString,Any}([("colordim",0), ("colorspace","Gray")]))
        buf = map(Images.mapinfo(UInt32, imgray), imgray) # Images.uint32color(imgray)
        @fact buf --> gray32
        rgb = RGB{Float64}[RGB(g, g, g) for g in gray]
        buf = map(Images.mapinfo(UInt32, rgb), rgb) # Images.uint32color(rgb)
        @fact buf --> gray32
        r = red(rgb)
        @fact r --> gray
        img = Images.Image(reinterpret(RGB24, gray32)) # , ["colordim"-->0, "colorspace"=>"RGB24"])
        buf = map(Images.mapinfo(UInt32, img), img) # Images.uint32color(img)
        @fact buf --> gray32
        rgb = repeat(gray, outer=[1,3])
        img = Images.Image(rgb, Dict{ASCIIString,Any}([("colordim",2), ("colorspace","RGB"), ("spatialorder",["x"])]))
        buf = map(Images.mapinfo(UInt32, img), img) # Images.uint32color(img)
        @fact buf --> gray32
        g = green(img)
        @fact g --> gray
        rgb = repeat(gray', outer=[3,1])
        img = Images.Image(rgb, Dict{ASCIIString,Any}([("colordim",1), ("colorspace","RGB"), ("spatialorder",["x"])]))
        buf = map(Images.mapinfo(UInt32, img), img) # Images.uint32color(img)
        @fact buf --> gray32
        b = blue(img)
        @fact b --> gray
    end

    context("Map and indexed images") do
        img = Images.ImageCmap([1 2 3; 3 2 1], [RGB{UFixed16}(1.0,0.6,0.4), RGB{UFixed16}(0.2, 0.4, 0.6), RGB{UFixed16}(0.5,0.5,0.5)])
        mapi = MapNone(RGB{UFixed8})
        imgd = map(mapi, img)
        cmap = [RGB{UFixed8}(1.0,0.6,0.4), RGB{UFixed8}(0.2, 0.4, 0.6), RGB{UFixed8}(0.5,0.5,0.5)]
        @fact imgd --> reshape(cmap[[1,3,2,2,3,1]], (2,3))
    end
end


text:
if occursin(r"RvSpectMLEcoSystem$", pwd())   cd("EXPRES-StellarSignals")   end
using Pkg
Pkg.activate(".")

verbose = true
 make_plots = false
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectML
 if verbose   println("# Loading other packages")    end
 using CSV, DataFrames, Query, StatsBase, Statistics, Dates
 if verbose   println("# Loading unpackaged code from this repo")    end
 include("io.jl")

all_spectra = include(joinpath(pkgdir(EchelleInstruments),"examples/read_expres_data_101501.jl"))
 if typeof(get_inst(all_spectra)) <: AnyEXPRES   continuum_normalize_spectra!(all_spectra)   end
 #if typeof(get_inst(all_spectra)) <: AnyEXPRES   EXPRES.blaze_normalize_spectra!(all_spectra)   end
 #RvSpectMLBase.discard_tellurics(all_spectra)
 RvSpectMLBase.discard_blaze(all_spectra)
 RvSpectMLBase.discard_continuum(all_spectra)
 GC.gc()
 ref_obs_idx = RvSpectMLBase.choose_obs_idx_for_init_guess(df_files_use,get_inst(all_spectra))
 tidx_epoch4 = df_files_use.expres_epoch.==4
 tidx_epoch5 = df_files_use.expres_epoch.==5

# Something buggy if include order 1.  Order 86 essentially ignored due to tellurics. First and last orders have NaN issues
lfc_orders = 43:72
  #max_orders = 12:83
  order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=lfc_orders, recalc=true )
  #order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=max_orders, recalc=true )
  order_list_timeseries_ref_obs_only = extract_chunklist_timeseries_with_subset_obs(order_list_timeseries, ref_obs_idx)

function calc_epoch_offset_vector(rvs_ccf, σ_rvs, tidx_epoch4, tidx_epoch5)
    vcat(rvs_ccf[tidx_epoch4].-mean(rvs_ccf[tidx_epoch4], weights(1.0 ./σ_rvs[tidx_epoch4].^2)), (rvs_ccf[tidx_epoch5].-mean(rvs_ccf[tidx_epoch5],weights(1.0 ./σ_rvs[tidx_epoch5].^2))))
end

linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks","G9.espresso.mas")
 line_list_espresso = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc, recalc=true, verbose=false)
 (ccfs, v_grid) = ccf_total(order_list_timeseries_ref_obs_only, line_list_espresso, pipeline_plan,  mask_scale_factor=2.0, ccf_mid_velocity=ccf_mid_velocity, v_step=100, recalc=true)
 line_width_50 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.5)
 line_sigma = line_width_50/sqrt(2*log(2))
 line_width_05 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.05)
 σ_rvs_formal_espresso = map(obsid->calc_formal_rv_precission(all_spectra[obsid],RvSpectMLBase.make_chunk_list_from_loc_df(all_spectra[obsid],get_inst(all_spectra), line_list_espresso),smooth_factor=2.0),1:length(all_spectra) )
 mean(σ_rvs_formal_espresso), median(σ_rvs_formal_espresso), extrema(σ_rvs_formal_espresso)

((ccfs_espresso, ccf_vars_espresso), v_grid) = ccf_total(order_list_timeseries, line_list_espresso, pipeline_plan,  mask_scale_factor=8.0, range_no_mask_change=line_width_05, ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=true, recalc=true)
 alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
 rvs_ccf_espresso = RvSpectML.calc_rvs_from_ccf_total(ccfs_espresso, ccf_vars_espresso, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
 #σ_rvs_espresso = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))
 chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_telluric_free(all_spectra, line_list_espresso, min_pixels_in_chunk=5, verbose=false)
 dont_need_to!(pipeline_plan,:extract_orders)
 rv_offset_epoch4_minus_epoch5 = calc_epoch_offset_vector(rvs_ccf_espresso, σ_rvs_formal_espresso, tidx_epoch4, tidx_epoch5 )
 set_rv_est!(chunk_list_timeseries, rvs_ccf_espresso .+ tidx_epoch5.*rv_offset_epoch4_minus_epoch5 )
 chunk_grids = RvSpectML.make_grids_for_chunklist_timeseries(chunk_list_timeseries)
 idx_use_pca = trues(length(all_spectra)); idx_use_pca[2] = false
 dcpca_result_espresso = RvSpectML.DCPCA.clean_rvs_dcpca(chunk_list_timeseries, rvs_ccf_espresso, idx_use_pca=idx_use_pca, chunk_grids=chunk_grids, num_samples_σ_scores=0)

 σ_rvs_clean_espresso = σ_rvs_formal_espresso # sqrt.(σ_rvs_formal_espresso.^2 .+ dcpca_result_espresso.σ_rvs_dirty.^2 )
 df_espresso = make_results_df(order_list_timeseries.times,dcpca_result_espresso.rvs_clean,σ_rvs_clean_espresso,dcpca_result_espresso.rvs_dirty,NaN.+dcpca_result_espresso.σ_rvs_dirty,dcpca_result_espresso.scores,NaN.+dcpca_result_espresso.σ_scores)
 CSV.write(joinpath("pennstate","101501","101501_pennstate_dcpca0_results.csv"), df_espresso)

linelist_for_ccf_fn_w_path = joinpath("data","101501","HD101501_f75s75v3_mask.csv")
 line_list_clean1 = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc, recalc=true, verbose=false)
 σ_rvs_formal_clean1 = map(obsid->calc_formal_rv_precission(all_spectra[obsid],RvSpectMLBase.make_chunk_list_from_loc_df(all_spectra[obsid],get_inst(all_spectra), line_list_clean1),smooth_factor=2.0),1:length(all_spectra) )
 mean(σ_rvs_formal_clean1), median(σ_rvs_formal_clean1), extrema(σ_rvs_formal_clean1)
 t_start_1 = now()
 ((ccfs_clean1, ccf_vars_clean1), v_grid) = ccf_total(order_list_timeseries, line_list_clean1, pipeline_plan,  mask_scale_factor=8.0, range_no_mask_change=line_width_05, ccf_mid_velocity=ccf_mid_velocity, v_step=100, calc_ccf_var=true, recalc=true)
 alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
 rvs_ccf_clean1 = calc_rvs_from_ccf_total(ccfs_clean1, ccf_vars_clean1, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
 #σ_rvs_clean1 = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))
 chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_telluric_free(all_spectra, line_list_clean1, min_pixels_in_chunk=5, verbose=false)
 dont_need_to!(pipeline_plan,:extract_orders)
 rv_offset_epoch4_minus_epoch5 = calc_epoch_offset_vector(rvs_ccf_clean1, σ_rvs_formal_clean1, tidx_epoch4, tidx_epoch5 )
 set_rv_est!(chunk_list_timeseries, rvs_ccf_clean1 .+ tidx_epoch5.*rv_offset_epoch4_minus_epoch5 )
 chunk_grids = RvSpectML.make_grids_for_chunklist_timeseries(chunk_list_timeseries)
 idx_use_pca = trues(length(all_spectra)); idx_use_pca[2] = false
 dcpca_result_clean1 = RvSpectML.DCPCA.clean_rvs_dcpca(chunk_list_timeseries, rvs_ccf_clean1, idx_use_pca=idx_use_pca, chunk_grids=chunk_grids, num_samples_σ_scores=0)
 σ_rvs_clean_clean1 = σ_rvs_formal_clean1 # sqrt.(σ_rvs_formal_clean1.^2 .+ dcpca_result_clean1.σ_rvs_dirty.^2 )
 df_clean1 = make_results_df(order_list_timeseries.times,dcpca_result_clean1.rvs_clean,σ_rvs_clean_clean1,dcpca_result_clean1.rvs_dirty,NaN.+dcpca_result_clean1.σ_rvs_dirty,dcpca_result_clean1.scores,NaN.+dcpca_result_clean1.σ_scores)
 t_stop_1 = now()
 CSV.write(joinpath("pennstate","101501","101501_pennstate_dcpca1_results.csv"), df_clean1)
 t_stop_1-t_start_1

linelist_for_ccf_fn_w_path = joinpath("data","101501","HD101501_f80s85v2_mask.csv")
  line_list_clean2 = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = RvSpectMLBase.max_bc, recalc=true, verbose=false)
  σ_rvs_formal_clean2 = map(obsid->calc_formal_rv_precission(all_spectra[obsid],RvSpectMLBase.make_chunk_list_from_loc_df(all_spectra[obsid],get_inst(all_spectra), line_list_clean2),smooth_factor=2.0),1:length(all_spectra) )
  mean(σ_rvs_formal_clean2), median(σ_rvs_formal_clean2), extrema(σ_rvs_formal_clean2)
  t_start_2 = now()
  ((ccfs_clean2, ccf_vars_clean2), v_grid) = ccf_total(order_list_timeseries, line_list_clean2, pipeline_plan,  mask_scale_factor=8.0, range_no_mask_change=line_width_05, ccf_mid_velocity=ccf_mid_velocity, v_step=100, calc_ccf_var=true, recalc=true)
  alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2)
  rvs_ccf_clean2 = calc_rvs_from_ccf_total(ccfs_clean2, ccf_vars_clean2, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
  #σ_rvs_clean2 = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))
  chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_telluric_free(all_spectra, line_list_clean2, min_pixels_in_chunk=5, verbose=false)
  dont_need_to!(pipeline_plan,:extract_orders)
  rv_offset_epoch4_minus_epoch5 = calc_epoch_offset_vector(rvs_ccf_clean2, σ_rvs_formal_clean2, tidx_epoch4, tidx_epoch5 )
  set_rv_est!(chunk_list_timeseries, rvs_ccf_clean2 .+ tidx_epoch5.*rv_offset_epoch4_minus_epoch5 )
  chunk_grids = RvSpectML.make_grids_for_chunklist_timeseries(chunk_list_timeseries)
  idx_use_pca = trues(length(all_spectra)); idx_use_pca[2] = false
  dcpca_result_clean2 = RvSpectML.DCPCA.clean_rvs_dcpca(chunk_list_timeseries, rvs_ccf_clean2, idx_use_pca=idx_use_pca, chunk_grids=chunk_grids, num_samples_σ_scores=0)
  σ_rvs_clean_clean2 = σ_rvs_formal_clean2 # sqrt.(σ_rvs_formal_clean2.^2 .+ dcpca_result_clean2.σ_rvs_dirty.^2 )
  df_clean2 = make_results_df(order_list_timeseries.times,dcpca_result_clean2.rvs_clean,σ_rvs_clean_clean2,dcpca_result_clean2.rvs_dirty,NaN.+dcpca_result_clean2.σ_rvs_dirty,dcpca_result_clean2.scores,NaN.+dcpca_result_clean2.σ_scores)
  t_stop_2 = now()
  CSV.write(joinpath("pennstate","101501","101501_pennstate_dcpca2_results.csv"), df_clean2)
  t_stop_2-t_start_2

std(dcpca_result_espresso.rvs_clean), std(dcpca_result_clean1.rvs_clean), std(dcpca_result_clean2.rvs_clean)

using Scalpels
function compute_scalpels_output(mnb::Integer)
  (rv_scalpels, scores, ) = calc_clean_rvs_scores_basis_scalpels(rvs_ccf, ccfs, num_basis=mnb, σ_rvs=σ_rvs)
  rv_scalpels_offsets = (rv_scalpels .- vcat(fill(mean(rv_scalpels[tidx_epoch4], weights(1.0 ./σ_rvs[tidx_epoch4].^2)), sum(tidx_epoch4)), fill(mean(rv_scalpels[tidx_epoch5],weights(1.0 ./σ_rvs[tidx_epoch5].^2)), sum(tidx_epoch5)) ) )

  println("max_num_basis = ", mnb, "  std RV = ", std(rv_scalpels_offsets) )

  output_df = DataFrame("Time [MJD]"=>order_list_timeseries.times, "RV_C [m/s]"=>rv_scalpels_offsets, "e_RV_C [m/s]"=>σ_rvs,
                      "RV_D [m/s]"=>rvs_ccf.-rv_scalpels_offsets, "e_RV_D [m/s]"=>σ_rvs )
  for i in 1:mnb
     output_df[!,"Score_" * string(i) * " [nondim]"] = vec(scores[:,i])
     output_df[!,"e_Score_" * string(i) * " [nondim]"] = fill(NaN,length(rv_scalpels_offsets))
  end
  return output_df
end

num_basis_scalpels = 2
 ccfs = ccfs_espresso; ccf_vars = ccf_vars_espresso; rvs_ccf = rvs_ccf_espresso; σ_rvs = σ_rvs_formal_espresso; println("# Setting ccfs to ESPRESSO mask ", num_basis_scalpels, " basis vectors");
 @time output_df = compute_scalpels_output(num_basis_scalpels)
 CSV.write(joinpath("pennstate","101501","101501_pennstate_scalpels0_results.csv"), output_df)

 ccfs = ccfs_clean1; ccf_vars = ccf_vars_clean1; rvs_ccf = rvs_ccf_clean1; σ_rvs = σ_rvs_formal_clean1;  println("# Setting ccfs to clean1 mask");
 output_df = compute_scalpels_output(num_basis_scalpels)
 CSV.write(joinpath("pennstate","101501","101501_pennstate_scalpels1_results.csv"), output_df)

 ccfs = ccfs_clean2; ccf_vars = ccf_vars_clean2; rvs_ccf = rvs_ccf_clean2; σ_rvs = σ_rvs_formal_clean2;  println("# Setting ccfs to clean2 mask");
 output_df = compute_scalpels_output(num_basis_scalpels)
 CSV.write(joinpath("pennstate","101501","101501_pennstate_scalpels2_results.csv"), output_df)

println("# Done")


text:
facts("Testing simple vote with single voter") do
    srand(1234)
    agents = [Agent(100)]
    @fact TcrCore.simpleVote(0, 30, agents) --> [true, agents, []]
    @fact TcrCore.simpleVote(50, 30, agents) --> [false, [], agents]
end

facts("Testing simple vote with diversified voters") do
    srand(1234)
    highAccAgent = Agent(100)
    lowAccAgent = Agent(20)
    agents = [lowAccAgent, highAccAgent]

    @fact TcrCore.simpleVote(50, 50, agents) --> [true, agents, []]
    @fact TcrCore.simpleVote(50, 40, agents) --> [false, [lowAccAgent], [highAccAgent]]
    @fact TcrCore.simpleVote(50, 0, agents) --> [false, [], agents]
end

facts("Testing binary token vote with diversified voters") do
    srand(1234)
    highAccAgent = Agent(100)
    lowAccAgent = Agent(20)
    agents = [lowAccAgent, highAccAgent]

    #Both agents vote
    @fact TcrCore.binaryTokenVote(50, 40, agents) --> [false, [lowAccAgent], [highAccAgent]]

    #Only high accuracy votes
    lowAccAgent.balance = 0
    @fact TcrCore.binaryTokenVote(50, 40, agents) --> [false, [], [highAccAgent]]

    #Only high accuracy votes
    lowAccAgent.balance = 50
    highAccAgent.balance = 0
    @fact TcrCore.binaryTokenVote(50, 40, agents) --> [true, [lowAccAgent], []]
end

facts("Testing pro-rata token vote with diversified voters") do
    srand(1234)
    highAccAgent = Agent(100)
    lowAccAgent = Agent(20)
    agents = [lowAccAgent, highAccAgent]

    #Two agents with the same amount of tokens
    @fact TcrCore.proRataTokenVote(50, 40, agents) --> [false, [lowAccAgent], [highAccAgent]]

    #Agents differs in tokens
    lowAccAgent.balance = 51
    @fact TcrCore.proRataTokenVote(50, 40, agents) --> [true, [lowAccAgent], [highAccAgent]]

    #One agent has got more tokens than two others in total
    srand(1234)
    highAccAgent.balance = 50
    lowAccAgent.balance = 101
    @fact TcrCore.proRataTokenVote(50, 40, [lowAccAgent, highAccAgent, highAccAgent]) --> [true, [lowAccAgent], [highAccAgent, highAccAgent]]
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related with highlights.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Highlight searches in strings" begin
    str      = "Test high\e[1mlight\e[0m in a string with no underlines."
    expected = "Test \e[7mhighlight\e[0m\e[0m in a string with no underlines."
    hstr     = highlight_search(str, r"highlight")
    @test expected == hstr

    str      = "Test high\e[1mlight\e[0m in a string with no underlines."
    expected = "Test \e[30m\e[43mhighlight\e[0m\e[0m in a string with no underlines."
    hstr     = highlight_search(str, r"highlight"; active_match = 1)
    @test expected == hstr

    str      = "Test high\e[4mlight in a string with underlines."
    expected = "Test \e[7mhighlight\e[0m\e[4m in a string with underlines."
    hstr     = highlight_search(str, r"highlight")
    @test expected == hstr

    # No matches in the string.
    hstr = highlight_search(str, r"Nothing to match")
    @test hstr == str
end

@testset "Highlight searches in texts with multiple lines" begin
    str = """
        Test high\e[1mlight\e[0m in a string with no underlines.
        Test high\e[4mlight in a string with underlines\e[0m.
        Test another high\e[33mlight with colors.
        This is the last line."""
    lines = split(str, '\n')

    expected = """
        Test \e[7mhighlight\e[0m\e[0m in a string with no underlines.
        Test \e[7mhighlight\e[0m\e[4m in a string with underlines\e[0m.
        Test another \e[7mhighlight\e[0m\e[33m with colors.
        This is the last line."""
    hstr = highlight_search(lines, r"highlight")
    @test hstr == expected

    expected = """
        Test \e[7mhighlight\e[0m\e[0m in a string with no underlines.
        Test \e[30m\e[43mhighlight\e[0m\e[4m in a string with underlines\e[0m.
        Test another \e[7mhighlight\e[0m\e[33m with colors.
        This is the last line."""
    hstr = highlight_search(lines, r"highlight"; active_match = 2)
    @test hstr == expected

    expected = """
        Test \e[7mhighlight\e[0m\e[4m in a string with underlines\e[0m.
        Test another \e[30m\e[43mhighlight\e[0m\e[33m with colors."""
    hstr = highlight_search(
        lines,
        r"highlight";
        active_match = 3,
        start_line = 2,
        end_line = 3
    )
    @test hstr == expected

    # No matches in the string.
    expected = """
        Test high\e[1mlight\e[0m in a string with no underlines.
        Test high\e[4mlight in a string with underlines\e[0m.
        Test another high\e[33mlight with colors."""
    hstr = highlight_search(
        lines,
        r"nothing to match";
        active_match = 3,
        start_line = 1,
        end_line = 3
    )
    @test hstr == expected
end


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ 43f2b466-8c3b-40a3-8e33-5579401aaf1f
begin
	using Pkg
	Pkg.activate(temp=true)
	Pkg.add(url="https://github.com/queryverse/VegaLite.jl.git", rev="master")
	Pkg.add("RDatasets")

	using VegaLite, RDatasets
end

# ╔═╡ 94662366-dfff-499b-9aed-db167544d243
iris = dataset("datasets", "iris")

# ╔═╡ d9d15c0e-b00d-4558-a1da-183833f89215
iris |> @vlplot(:circle,
	x="PetalWidth", y="PetalLength",
	size="SepalLength")

# ╔═╡ e45a73a4-0c0d-4e52-a623-85528bf583e5
min_sl, max_sl = extrema(iris.SepalLength)

# ╔═╡ e9df01c9-4a82-44a5-a2a0-49836a582a4f
iris |> @vlplot(:circle,
	x="PetalWidth", y="PetalLength",
	size={"SepalLength",
		scale={domain=[min_sl, max_sl]}})

# ╔═╡ 7351e879-fb50-40b3-81c8-4a8ae22f839f
iris |> @vlplot(:circle,
	x="PetalWidth", y="PetalLength",
	tooltip="SepalLength")

# ╔═╡ 5e613450-dd3b-48d6-87c0-c140173c5b3e
iris |> @vlplot(:circle,
	x="PetalWidth", y="PetalLength",
	selection={brush={type="interval"}},
	color={condition={selection="brush",
		value="red"}, value="lightgrey"})

# ╔═╡ 30cd2d23-3a5a-4baf-aadf-de64566b4907
iris |> @vlplot(repeat={
	column=["PetalLength", "PetalWidth"],
	row=["SepalLength", "SepalWidth"]}) +
	@vlplot(:circle,
		x={field={repeat="column"},
			type="quantitative"},
		y={field={repeat="row"},
			type="quantitative"},
		selection={brush={type="interval"}},
    	color={condition={selection="brush",
			value="red"}, value="lightgrey"},
		width=250, height=200)

# ╔═╡ Cell order:
# ╠═43f2b466-8c3b-40a3-8e33-5579401aaf1f
# ╠═94662366-dfff-499b-9aed-db167544d243
# ╠═d9d15c0e-b00d-4558-a1da-183833f89215
# ╠═e45a73a4-0c0d-4e52-a623-85528bf583e5
# ╠═e9df01c9-4a82-44a5-a2a0-49836a582a4f
# ╠═7351e879-fb50-40b3-81c8-4a8ae22f839f
# ╠═5e613450-dd3b-48d6-87c0-c140173c5b3e
# ╠═30cd2d23-3a5a-4baf-aadf-de64566b4907


text:
export Polytope, Edge, LineSegment, QuadraticSegment, Face, Polygon, QuadraticPolygon,
       Triangle, Quadrilateral, QuadraticTriangle, QuadraticQuadrilateral, Cell,
       Polyhedron, QuadraticPolyhedron, Tetrahedron, Hexahedron, QuadraticTetrahedron,
       QuadraticHexahedron
export vertices, facets, ridges, peaks, alias_string, vertex_type, paramdim,
       isstraight

"""

    Polytope{K,P,N,T}

A `K`-polytope of order `P` with `N` vertices of type `T`.

## Aliases

```julia
# 1-polytope
const Edge                      = Polytope{1}
const LineSegment               = Edge{1,2}
const QuadraticSegment          = Edge{2,3}
# 2-polytope
const Face                      = Polytope{2}
const Polygon                   = Face{1}
const QuadraticPolygon          = Face{2}
const Triangle                  = Polygon{3}
const Quadrilateral             = Polygon{4}
const QuadraticTriangle         = QuadraticPolygon{6}
const QuadraticQuadrilateral    = QuadraticPolygon{8}
# 3-polytope
const Cell                      = Polytope{3}
const Polyhedron                = Cell{1}
const QuadraticPolyhedron       = Cell{2}
const Tetrahedron               = Polyhedron{4}
const Hexahedron                = Polyhedron{8}
const QuadraticTetrahedron      = QuadraticPolyhedron{10}
const QuadraticHexahedron       = QuadraticPolyhedron{20}
```
### Notes
- These are Lagrangian finite elements.
- This struct only supports the shapes found in "The Visualization Toolkit:
  An Object-Oriented Approach to 3D Graphics, 4th Edition, Chapter 8, Advanced
  Data Representation".
- See the VTK book for specific vertex ordering info, but generally vertices are
  ordered in a counterclockwise fashion, with vertices of the linear shape given
  first.
- See https://en.wikipedia.org/wiki/Polytope for help with terminology.
"""
struct Polytope{K,P,N,T}
    vertices::Vec{N,T}
    Polytope{K,P,N,T}(vertices::Vec{N,T}) where {K,P,N,T} = new{K,P,N,T}(vertices)
end

# type aliases
# 1-polytope
const Edge                      = Polytope{1}
const LineSegment               = Edge{1,2}
const QuadraticSegment          = Edge{2,3}
# 2-polytope
const Face                      = Polytope{2}
const Polygon                   = Face{1}
const QuadraticPolygon          = Face{2}
const Triangle                  = Polygon{3}
const Quadrilateral             = Polygon{4}
const QuadraticTriangle         = QuadraticPolygon{6}
const QuadraticQuadrilateral    = QuadraticPolygon{8}
# 3-polytope
const Cell                      = Polytope{3}
const Polyhedron                = Cell{1}
const QuadraticPolyhedron       = Cell{2}
const Tetrahedron               = Polyhedron{4}
const Hexahedron                = Polyhedron{8}
const QuadraticTetrahedron      = QuadraticPolyhedron{10}
const QuadraticHexahedron       = QuadraticPolyhedron{20}

# constructors
Polytope{K,P,N,T}(vertices...) where {K,P,N,T} = Polytope{K,P,N,T}(Vec{N,T}(vertices))
Polytope{K,P,N}(vertices::Vec{N,T}) where {K,P,N,T} = Polytope{K,P,N,T}(vertices)
Polytope{K,P,N}(vertices...) where {K,P,N} = Polytope{K,P,N}(Vec(vertices))
function Polytope{K,P,N}(v::Vector{T}) where {K,P,N,T}
    return Polytope{K,P,N}(Vec{length(v),T}(v))
end
function Polytope{K,P,N,T}(v::Vector) where {K,P,N,T}
    return Polytope{K,P,N,T}(Vec{length(v),T}(v))
end
function Polytope{K,P,N,T}(v::SVector{N}) where {K,P,N,T}
    return Polytope{K,P,N,T}(Vec{N,T}(v...))
end

# Convert SVector vals
function Base.convert(::Type{Polytope{K,P,N,T}}, v::SVector{N}) where {K,P,N,T}
    return Polytope{K,P,N,T}(v...)
end
# Convert Polytope vals
function Base.convert(::Type{Polytope{K,P,N,T}},
                          p::Polytope{K,P,N}) where {K,P,N,T}
    return Polytope{K,P,N,T}(p.vertices...)
end
# Convert SVector to general polytope
# Primarily for non-homogenous meshes
function Base.convert(::Type{Polytope{K,P,N,T} where {N}},
                      v::SVector{V}) where {K,P,T,V}
    return Polytope{K,P,V,T}(v...)
end

Base.getindex(poly::Polytope, i::Int) = Base.getindex(poly.vertices, i)

paramdim(::Type{<:Polytope{K}}) where {K} = K
vertex_type(::Type{Polytope{K,P,N,T}}) where {K,P,N,T} = T
vertex_type(::Polytope{K,P,N,T}) where {K,P,N,T} = T

vertices(p::Polytope) = p.vertices

peaks(p::Polytope{3}) = vertices(p)

ridges(p::Polytope{2}) = vertices(p)
ridges(p::Polytope{3}) = edges(p)

facets(p::Polytope{1}) = vertices(p)
facets(p::Polytope{2}) = edges(p)
facets(p::Polytope{3}) = faces(p)

function alias_string(::Type{P}) where {P<:Polytope}
    P <: LineSegment            && return "LineSegment"
    P <: QuadraticSegment       && return "QuadraticSegment"
    P <: Triangle               && return "Triangle"
    P <: Quadrilateral          && return "Quadrilateral"
    P <: QuadraticTriangle      && return "QuadraticTriangle"
    P <: QuadraticQuadrilateral && return "QuadraticQuadrilateral"
    P <: Tetrahedron            && return "Tetrahedron"
    P <: Hexahedron             && return "Hexahedron"
    P <: QuadraticTetrahedron   && return "QuadraticTetrahedron"
    P <: QuadraticHexahedron    && return "QuadraticHexahedron"
    # fallback on default
    return "$(P)"
end

# If we think of the polytopes as sets, p₁ ∩ p₂ = p₁ and p₁ ∩ p₂ = p₂ implies p₁ = p₂
function Base.:(==)(l₁::LineSegment{T}, l₂::LineSegment{T}) where {T}
    return (l₁[1] === l₂[1] && l₁[2] === l₂[2]) ||
           (l₁[1] === l₂[2] && l₁[2] === l₂[1])
end
Base.:(==)(t₁::Triangle, t₂::Triangle) = return all(v->v ∈ t₂.vertices, t₁.vertices)
Base.:(==)(t₁::Tetrahedron, t₂::Tetrahedron) = return all(v->v ∈ t₂.vertices, t₁.vertices)
function Base.:(==)(q₁::QuadraticSegment{T}, q₂::QuadraticSegment{T}) where {T}
    return q₁[3] === q₂[3] &&
          (q₁[1] === q₂[1] && q₁[2] === q₂[2])  ||
          (q₁[1] === q₂[2] && q₁[2] === q₂[1])
end

isstraight(::LineSegment) = true

"""
    isstraight(q::QuadraticSegment)

Return if the quadratic segment is effectively straight.
(If P₃ is at most ϵ_Point distance from LineSegment(P₁,P₂))
"""
function isstraight(q::QuadraticSegment{T}) where {T<:Point}
    # Project P₃ onto the line from P₁ to P₂, call it P₄
    𝘃₁₃ = q[3] - q[1]
    𝘃₁₂ = q[2] - q[1]
    v₁₂ = norm²(𝘃₁₂)
    𝘃₁₄ = (𝘃₁₃ ⋅ 𝘃₁₂)*inv(v₁₂)*𝘃₁₂
    # Determine the distance from P₃ to P₄ (P₄ - P₃ = P₁ + 𝘃₁₄ - P₃ = 𝘃₁₄ - 𝘃₁₃)
    d² = norm²(𝘃₁₄ - 𝘃₁₃)
    return d² < ϵ_Point^2
end

# Show aliases when printing
function Base.show(io::IO, poly::Polytope)
    print(io, alias_string(typeof(poly)),"(",vertices(poly), ")")
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
module transD_GP
include("GP.jl")
include("TransD_GP_MCMC.jl")
include("AbstractOperator.jl")
include("MCMC_Driver.jl")
include("GradientInversion.jl")
include("GeophysOperator.jl")
end


text:
using HypothesisTests, Base.Test
using Distributions

# This is always the null in our tests.
null = Normal(0.0, 1.0)

## ONE SAMPLE T-TEST

# One sample
x = -5:5
@test pvalue(OneSampleZTest(x)) == 1

x = -5:10
m, s, n = mean(x), std(x), length(x)
se = s / sqrt(n)
z = (m - 0) / se

tst = OneSampleZTest(x)
@test_approx_eq(pvalue(tst), 2 * min(cdf(null, z), ccdf(null, z)))
@test_approx_eq(pvalue(tst; tail=:left), cdf(null, z))
@test_approx_eq(pvalue(tst; tail=:right), ccdf(null, z))
show(IOBuffer(), tst)

tst = OneSampleZTest(m, s, n)
@test_approx_eq(pvalue(tst), 2 * min(cdf(null, z), ccdf(null, z)))
@test_approx_eq(ci(tst)[1], m + quantile(null, 0.05 / 2) * se)
@test_approx_eq(ci(tst)[2], m + cquantile(null, 0.05 / 2) * se)
@test_approx_eq(ci(tst, 0.10)[1], m + quantile(null, 0.10 / 2) * se)
@test_approx_eq(ci(tst, 0.10)[2], m + cquantile(null, 0.10 / 2) * se)
@test_approx_eq(ci(tst; tail=:left)[1], -Inf)
@test_approx_eq(ci(tst; tail=:left)[2], m + cquantile(null, 0.05) * se)
@test_approx_eq(ci(tst; tail=:right)[1], m + quantile(null, 0.05) * se)
@test_approx_eq(ci(tst; tail=:right)[2], Inf)
show(IOBuffer(), tst)

x = -10:5
m, s, n = mean(x), std(x), length(x)
se = s / sqrt(n)
z = (m - 0) / se

tst = OneSampleZTest(x)
@test_approx_eq(pvalue(tst), 2 * min(cdf(null, z), ccdf(null, z)))
@test_approx_eq(pvalue(tst; tail=:left), cdf(null, z))
@test_approx_eq(pvalue(tst; tail=:right), ccdf(null, z))
show(IOBuffer(), tst)

tst = OneSampleZTest(m, s, n)
@test_approx_eq(pvalue(tst), 2 * min(cdf(null, z), ccdf(null, z)))
@test_approx_eq(ci(tst)[1], m + quantile(null, 0.05 / 2) * se)
@test_approx_eq(ci(tst)[2], m + cquantile(null, 0.05 / 2) * se)
@test_approx_eq(ci(tst, 0.10)[1], m + quantile(null, 0.10 / 2) * se)
@test_approx_eq(ci(tst, 0.10)[2], m + cquantile(null, 0.10 / 2) * se)
@test_approx_eq(ci(tst; tail=:left)[1], -Inf)
@test_approx_eq(ci(tst; tail=:left)[2], m + cquantile(null, 0.05) * se)
@test_approx_eq(ci(tst; tail=:right)[1], m + quantile(null, 0.05) * se)
@test_approx_eq(ci(tst; tail=:right)[2], Inf)
show(IOBuffer(), tst)

# Paired samples
x, y = [1, 1, 2, 1, 0], [0, 1, 1, 1, 0]
m, s, n = mean(x - y), std(x - y), length(x - y)
se = s / sqrt(n)
z = (m - 0) / se
tst = OneSampleZTest(x, y)
@test_approx_eq(pvalue(tst), 2 * min(cdf(null, z), ccdf(null, z)))

## TWO SAMPLE Z-TESTS

a1 = [30.02, 29.99, 30.11, 29.97, 30.01, 29.99]
a2 = [29.89, 29.93, 29.72, 29.98, 30.02, 29.98]

tst = EqualVarianceZTest(a1, a2)
m1, s1sq, n1 = mean(a1), var(a1), length(a1)
m2, s2sq, n2 = mean(a2), var(a2), length(a2)
xbar = (m1 - m2)
avg_var = (n1 - 1) / (n1 + n2 - 2) * s1sq + (n2 - 1) / (n1 + n2 - 2) * s2sq
se = sqrt(avg_var / n1 + avg_var / n2)
z = xbar / se

@test_approx_eq(tst.z, z)
@test_approx_eq(pvalue(tst), 2 * min(cdf(null, z), ccdf(null, z)))
@test_approx_eq(pvalue(tst; tail=:left), cdf(null, z))
@test_approx_eq(pvalue(tst; tail=:right), ccdf(null, z))
@test_approx_eq(ci(tst)[1], xbar + quantile(null, 0.05 / 2) * se)
@test_approx_eq(ci(tst)[2], xbar + cquantile(null, 0.05 / 2) * se)
show(IOBuffer(), tst)

tst = UnequalVarianceZTest(a1, a2)
se = sqrt(s1sq / n1 + s2sq / n2)
z = xbar / se
@test_approx_eq(tst.z, z)
@test_approx_eq(pvalue(tst), 2 * min(cdf(null, z), ccdf(null, z)))
@test_approx_eq(pvalue(tst; tail=:left), cdf(null, z))
@test_approx_eq(pvalue(tst; tail=:right), ccdf(null, z))
@test_approx_eq(ci(tst)[1], xbar + quantile(null, 0.05 / 2) * se)
@test_approx_eq(ci(tst)[2], xbar + cquantile(null, 0.05 / 2) * se)
show(IOBuffer(), tst)


text:
using KernelAbstractions
using Test

import KernelAbstractions.NDIteration: NDRange, StaticSize, DynamicSize

@kernel function index(A)
    I  = @index(Global, NTuple)
    @show A[I...]
end
kernel = index(CPU(), DynamicSize(), DynamicSize())
iterspace = NDRange{1, StaticSize{(128,)}, StaticSize{(8,)}}();
ctx = KernelAbstractions.mkcontext(kernel, 1, nothing, iterspace, Val(KernelAbstractions.NoDynamicCheck()))

@test KernelAbstractions.Cassette.overdub(ctx, KernelAbstractions.__index_Global_NTuple, CartesianIndex(1)) == (1,)


text:
module GPS2XY

export gps_dist
export gps_bearing
export gps_offset
#export gps_sim
#export parse_log

# deprecated in Julia v0.6
#typealias LocTuple NTuple{2, Float64}
LocTuple  = NTuple{2, Float64}

"""
`gps_dist(latlon1::LocTuple, latlon2::LocTuple)`

Returns a scalar distance between these two points (I think it won't be -ve)
"""
function gps_dist(latlon1::LocTuple, latlon2::LocTuple)
	lat1 = latlon1[1]
	lon1 = latlon1[2]
	lat2 = latlon2[1]
	lon2 = latlon2[2]

	a = sind( (lat2-lat1)/2 )^2 + cosd(lat1)*cosd(lat2) * sind((lon2-lon1)/2)^2
	c = 2atan2(sqrt(a), sqrt(1-a))
	R = 6371000   # earth distance in meters
	d = R*c
	return d
end

"""
`gps_bearing(latlon1::LocTuple, latlon2::LocTuple)`

where `LocTuple` is a tuple of two floats.
Returns bearing from latlon1 to latlon2, east of north.
"""
function gps_bearing(latlon1::LocTuple, latlon2::LocTuple)
	lat1 = latlon1[1]
	lon1 = latlon1[2]
	lat2 = latlon2[1]
	lon2 = latlon2[2]

	y = sind(lon2-lon1)*cosd(lat2)
	x = cosd(lat1)*sind(lat2) - sind(lat1)*cosd(lat2)*cosd(lon2-lon1)
	return mod(rad2deg(atan2(y,x)), 360.0)
end

"""
`gps_offset(latlon1::LocTuple, latlon2::LocTuple)`

Returns the dx, dy offset from `latlon1` to `latlon2`.
"""
function gps_offset(latlon1::LocTuple, latlon2::LocTuple)
	d = gps_dist(latlon1, latlon2)
	b = gps_bearing(latlon1, latlon2)
	dx = d * sind(b)
	dy = d * cosd(b)
	return dx,dy
end

end # module


text:
# From StatsMakie
function _linear(x::AbstractVector{T}, y::AbstractVector;
                 npoints = 100, wts = similar(x, 0),
                 interval = length(wts) > 0 ? nothing : :confidence) where T
    # Note: confidence intervals are currently not supported for WLS in GLM.jl
    try
        y = collect(y)
        x = collect(x)
        wts = collect(wts)
        lin_model = GLM.lm([ones(T, length(x)) x], y, wts=wts)
        x_min, x_max = extrema(x)
        x_new = range(x_min, x_max, length = npoints)
        pred = GLM.predict(lin_model,
                           [ones(T, npoints) x_new],
                           interval=interval)
        if !isnothing(interval)
            y_new, lower, upper = pred
        else
            y_new = pred
        end
        # the GLM predictions always return matrices
        x, y = x_new, vec(y_new)
        lines = style(x, y)
        if !isnothing(interval)
            band = style(x, vec(lower), vec(upper))
            return spec(:Lines) * lines + spec(:Band, alpha = 0.2) * band
        else
            return spec(:Lines) * lines
        end
    catch e
        @warn "Linear fit not possible for the given data"
        return AlgebraicList()
    end
end

const linear = Analysis(_linear)

function _smooth(x, y; npoints = 100, kwargs...)
    min, max = extrema(x)
    min < max || return AlgebraicList()
    model = Loess.loess(Float64.(x), Float64.(y); kwargs...)
    us = collect(range(min, stop = max, length = npoints))
    vs = Loess.predict(model, us)
    return spec(:Lines) * style(us, vs)
end

const smooth = Analysis(_smooth)


text:
# uniform interface for model estimation

export Estimator, MLEstimator
export nsamples, estimate

@compat abstract type Estimator{D<:Distribution} end

nsamples{D<:UnivariateDistribution}(e::Estimator{D}, x::Array) = length(x)
nsamples{D<:MultivariateDistribution}(e::Estimator{D}, x::Matrix) = size(x, 2)

type MLEstimator{D<:Distribution} <: Estimator{D} end
MLEstimator{D<:Distribution}(::Type{D}) = MLEstimator{D}()

estimate{D<:Distribution}(e::MLEstimator{D}, x) = fit_mle(D, x)
estimate{D<:Distribution}(e::MLEstimator{D}, x, w) = fit_mle(D, x, w)


text:
test2{T, T2}(a::T, b::T2) = a + b
#

#function ast_macro(f, types)
function test{T}(a::T, b)
    c = a+b::T
    if c == 10
        for i=1:7
            c += rand(Bool) ? 1 : 2
        end
    end
    f = (x,y) -> x^y
    c, f
end

function get_ast(li::LambdaInfo)
    ast = li.code
    if isa(ast, Vector{UInt8})
        return Base.uncompressed_ast(li)
    end
    ast
end

f = test
types = (Int, Int)
typed_lam = get_lambda(code_typed, f, types)
ast = get_ast(typed_lam)

ast, str = macro_form(f, types)
params = get_static_parameters(f, types)

# typed_ast = get_ast(typed_lam)
# f = test
# types = (Int, Int)
# lambda = get_lambda(code_lowered, f, types)

# ast = get_ast(lambda)

# function map_if_assignement(x::Expr)
#     if x.head == :(=)
#     end
# end