text:
###############################################################################
#
#   LaurentSeries.jl : Laurent series over rings and fields,
#                      capped relative precision
#
###############################################################################

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

doc"""
    O{T <: RingElement}(a::LaurentSeriesElem{T})
> Returns $0 + O(x^\mbox{deg}(a))$. Usually this function is called with $x^n$
> as parameter. Then the function returns the power series $0 + O(x^n)$, which
> can be used to set the precision of a power series when constructing it.
"""
function O(a::LaurentSeriesElem{T}) where T <: RingElement
   val = pol_length(a) + valuation(a) - 1
   return parent(a)(Array{T}(0), 0, val, val)
end

parent_type(::Type{T}) where {S <: RingElement, T <: LaurentSeriesRingElem{S}} = LaurentSeriesRing{S}

parent_type(::Type{T}) where {S <: FieldElement, T <: LaurentSeriesFieldElem{S}} = LaurentSeriesField{S}

doc"""
    parent(a::LaurentSeriesElem)
> Return the parent of the given power series.
"""
parent(a::LaurentSeriesElem) = a.parent

elem_type(::Type{T}) where {S <: RingElement, T <: LaurentSeriesRing{S}} = LaurentSeriesRingElem{S}

elem_type(::Type{T}) where {S <: FieldElement, T <: LaurentSeriesField{S}} = LaurentSeriesFieldElem{S}

doc"""
    base_ring(R::LaurentSeriesRing)
> Return the base ring of the given power series ring.
"""
base_ring(R::LaurentSeriesRing{T}) where T <: RingElement = R.base_ring::parent_type(T)

doc"""
    base_ring(R::LaurentSeriesField)
> Return the base ring of the given power series ring.
"""
base_ring(R::LaurentSeriesField{T}) where T <: FieldElement = R.base_ring::parent_type(T)

doc"""
    base_ring(a::LaurentSeriesElem)
> Return the base ring of the power series ring of the given power series.
"""
base_ring(a::LaurentSeriesElem) = base_ring(parent(a))

function isdomain_type(::Type{T}) where {S <: RingElement, T <: LaurentSeriesElem{S}}
   return isdomain_type(S)
end

isexact_type(a::Type{T}) where T <: LaurentSeriesElem = false

doc"""
    var(a::LaurentSeriesRing)
> Return the internal name of the generator of the power series ring. Note that
> this is returned as a `Symbol` not a `String`.
"""
var(a::LaurentSeriesRing) = a.S

doc"""
    var(a::LaurentSeriesField)
> Return the internal name of the generator of the power series ring. Note that
> this is returned as a `Symbol` not a `String`.
"""
var(a::LaurentSeriesField) = a.S

function check_parent(a::LaurentSeriesElem, b::LaurentSeriesElem)
   parent(a) != parent(b) &&
             error("Incompatible power series rings in power series operation")
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function Base.hash(a::LaurentSeriesElem, h::UInt)
   b = 0xb163af5694734274%UInt
   for i in 0:pol_length(a) - 1
      b = xor(b, hash(polcoeff(a, i), h), h)
      b = (b << 1) | (b >> (sizeof(Int)*8 - 1))
   end
   return b
end

doc"""
    pol_length(a::LaurentSeriesElem)
> Return the length of the polynomial underlying the given power series. This
> will be zero if the power series has no nonzero terms.
"""
pol_length(a::LaurentSeriesElem) = a.length

doc"""
    precision(a::LaurentSeriesElem)
> Return the precision of the given power series in absolute terms. This will
> be the sum of the valuation and the length of the underlying polynomial.
"""
precision(a::LaurentSeriesElem) = a.prec

doc"""
    valuation(a::LaurentSeriesElem)
> Return the valuation of the given power series, i.e. the degree of the first
> nonzero term (or the precision if it is arithmetically zero).
"""
valuation(a::LaurentSeriesElem) = a.val

doc"""
    max_precision(R::LaurentSeriesRing)
> Return the maximum relative precision of power series in the given power
> series ring.
"""
max_precision(R::LaurentSeriesRing) = R.prec_max

doc"""
    max_precision(R::LaurentSeriesField)
> Return the maximum relative precision of power series in the given power
> series ring.
"""
max_precision(R::LaurentSeriesField) = R.prec_max

function normalise(a::LaurentSeriesElem, len::Int)
   while len > 0 && iszero(a.coeffs[len])
      len -= 1
   end
   return len
end

function set_length!(a::LaurentSeriesElem, len::Int)
   a.length = len
end

function set_prec!(a::LaurentSeriesElem, prec::Int)
   a.prec = prec
end

function set_val!(a::LaurentSeriesElem, val::Int)
   a.val = val
end

function polcoeff(a::LaurentSeriesElem, n::Int)
   n < 0  && throw(DomainError())
   return n >= pol_length(a) ? zero(base_ring(a)) : a.coeffs[n + 1]
end

function coeff(a::LaurentSeriesElem, n::Int)
   if n < valuation(a)
      return base_ring(a)()
   else
      return polcoeff(a, n - valuation(a))
   end
end

doc"""
    zero(R::LaurentSeriesRing)
> Return $0 + O(x^n)$ where $n$ is the maximum precision of the power series
> ring $R$.
"""
zero(R::LaurentSeriesRing) = R(0)

doc"""
    zero(R::LaurentSeriesField)
> Return $0 + O(x^n)$ where $n$ is the maximum precision of the power series
> ring $R$.
"""
zero(R::LaurentSeriesField) = R(0)

doc"""
    one(R::LaurentSeriesRing)
> Return $1 + O(x^n)$ where $n$ is the maximum precision of the power series
> ring $R$.
"""
one(R::LaurentSeriesField) = R(1)

doc"""
    one(R::LaurentSeriesField)
> Return $1 + O(x^n)$ where $n$ is the maximum precision of the power series
> ring $R$.
"""
one(R::LaurentSeriesRing) = R(1)

doc"""
    gen(R::LaurentSeriesRing)
> Return the generator of the power series ring, i.e. $x + O(x^{n + 1})$ where
> $n$ is the maximum precision of the power series ring $R$.
"""
function gen(R::LaurentSeriesRing)
   S = base_ring(R)
   return R([S(1)], 1, max_precision(R) + 1, 1)
end

doc"""
    gen(R::LaurentSeriesField)
> Return the generator of the power series ring, i.e. $x + O(x^{n + 1})$ where
> $n$ is the maximum precision of the power series ring $R$.
"""
function gen(R::LaurentSeriesField)
   S = base_ring(R)
   return R([S(1)], 1, max_precision(R) + 1, 1)
end

doc"""
    iszero(a::LaurentSeriesElem)
> Return `true` if the given power series is arithmetically equal to zero to
> its current precision, otherwise return `false`.
"""
iszero(a::LaurentSeriesElem) = pol_length(a) == 0

doc"""
    isone(a::LaurentSeriesElem)
> Return `true` if the given power series is arithmetically equal to one to
> its current precision, otherwise return `false`.
"""
function isone(a::LaurentSeriesElem)
   return valuation(a) == 0 && pol_length(a) == 1 && isone(polcoeff(a, 0))
end

doc"""
    isgen(a::LaurentSeriesElem)
> Return `true` if the given power series is arithmetically equal to the
> generator of its power series ring to its current precision, otherwise return
> `false`.
"""
function isgen(a::LaurentSeriesElem)
   return valuation(a) == 1 && pol_length(a) == 1 && isone(polcoeff(a, 0))
end

doc"""
    isunit(a::LaurentSeriesElem)
> Return `true` if the given power series is arithmetically equal to a unit,
> i.e. is invertible, otherwise return `false`.
"""
isunit(a::LaurentSeriesElem) = valuation(a) == 0 && isunit(polcoeff(a, 0))

doc"""
    modulus{T <: ResElem}(a::LaurentSeriesElem{T})
> Return the modulus of the coefficients of the given power series.
"""
modulus(a::LaurentSeriesElem{T}) where {T <: ResElem} = modulus(base_ring(a))

function deepcopy_internal(a::LaurentSeriesElem{T}, dict::ObjectIdDict) where {T <: RingElement}
   coeffs = Array{T}(pol_length(a))
   for i = 1:pol_length(a)
      coeffs[i] = deepcopy(polcoeff(a, i - 1))
   end
   return parent(a)(coeffs, pol_length(a), precision(a), valuation(a))
end

function renormalize!(z::LaurentSeriesElem)
   i = 0
   zlen = pol_length(z)
   zval = valuation(z)
   zprec = precision(z)
   while i < zlen && iszero(polcoeff(z, i))
      i += 1
   end
   set_prec!(z, zprec)
   if i == zlen
      set_length!(z, 0)
      set_val!(z, zprec)
   else
      set_val!(z, zval + i)
      for j = 1:zlen - i
         z = setcoeff!(z, j - 1, polcoeff(z, j + i - 1))
      end
      set_length!(z, zlen - i)
   end
   return nothing
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, x::LaurentSeriesElem)
   len = pol_length(x)
   if len == 0
      print(io, zero(base_ring(x)))
   else
      coeff_printed = false
      for i = 0:len - 1
         c = polcoeff(x, i)
         bracket = needs_parentheses(c)
         if !iszero(c)
            if coeff_printed && !isnegative(c)
               print(io, "+")
            end
            if i + valuation(x) != 0
               if !isone(c) && (c != -1 || show_minus_one(elem_type(base_ring(x))))
                  if bracket
                     print(io, "(")
                  end
                  print(io, c)
                  if bracket
                     print(io, ")")
                  end
                  if i + valuation(x) != 0
                     print(io, "*")
                  end
               end
               if c == -1 && !show_minus_one(elem_type(base_ring(x)))
                  print(io, "-")
               end
               print(io, string(var(parent(x))))
               if i + valuation(x) != 1
                  print(io, "^")
                  print(io, valuation(x) + i)
               end
            else
               print(io, c)
            end
            coeff_printed = true
         end
      end
   end
   print(io, "+O(", string(var(parent(x))), "^", precision(x), ")")
end

function show(io::IO, a::LaurentSeriesRing)
   print(io, "Laurent series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

function show(io::IO, a::LaurentSeriesField)
   print(io, "Laurent series field in ", var(a), " over ")
   show(io, base_ring(a))
end

needs_parentheses(x::LaurentSeriesElem) = pol_length(x) > 1

isnegative(x::LaurentSeriesElem) = pol_length(x) <= 1 && isnegative(polcoeff(x, 0))

show_minus_one(::Type{LaurentSeriesElem{T}}) where {T <: RingElement} = show_minus_one(T)

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

doc"""
    -(a::LaurentSeriesElem)
> Return $-a$.
"""
function -(a::LaurentSeriesElem)
   len = pol_length(a)
   z = parent(a)()
   set_prec!(z, precision(a))
   set_val!(z, valuation(a))
   fit!(z, len)
   for i = 1:len
      z = setcoeff!(z, i - 1, -polcoeff(a, i - 1))
   end
   return z
end

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

doc"""
    +{T <: RingElement}(a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T})
> Return $a + b$.
"""
function +(a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T}) where {T <: RingElement}
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   vala = valuation(a)
   valb = valuation(b)
   valz = min(vala, valb)
   prec = min(precision(a), precision(b))
   mina = min(vala + lena, prec)
   minb = min(valb + lenb, prec)
   lenz = max(mina, minb) - valz
   R = base_ring(a)
   z = parent(a)()
   fit!(z, lenz)
   set_prec!(z, prec)
   set_val!(z, valz)
   if vala >= valb
      for i = 1:min(lenb, vala - valb)
         z = setcoeff!(z, i - 1, polcoeff(b, i - 1))
      end
      for i = lenb + 1:min(vala - valb, lenz)
         z = setcoeff!(z, i - 1, R())
      end
      for i = vala - valb + 1:lenb
         z = setcoeff!(z, i - 1, polcoeff(a, i - vala + valb - 1) + polcoeff(b, i - 1))
      end
      for i = max(lenb, vala - valb) + 1:lena + vala - valb
         z = setcoeff!(z, i - 1, polcoeff(a, i - vala + valb - 1))
      end
      for i = lena + vala - valb + 1:lenb
         z = setcoeff!(z, i - 1, polcoeff(b, i - 1))
      end
   else
      for i = 1:min(lena, valb - vala)
         z = setcoeff!(z, i - 1, polcoeff(a, i - 1))
      end
      for i = lena + 1:min(valb - vala, lenz)
         z = setcoeff!(z, i - 1, R())
      end
      for i = valb - vala + 1:lena
         z = setcoeff!(z, i - 1, polcoeff(a, i - 1) + polcoeff(b, i - valb + vala - 1))
      end
      for i = max(lena, valb - vala) + 1:lenb + valb - vala
         z = setcoeff!(z, i - 1, polcoeff(b, i - valb + vala - 1))
      end
      for i = lenb + valb - vala + 1:lena
         z = setcoeff!(z, i - 1, polcoeff(a, i - 1))
      end
   end
   set_length!(z, normalise(z, lenz))
   renormalize!(z)
   return z
end

doc"""
    -{T <: RingElement}(a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T})
> Return $a - b$.
"""
function -(a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T}) where {T <: RingElement}
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   vala = valuation(a)
   valb = valuation(b)
   valz = min(vala, valb)
   prec = min(precision(a), precision(b))
   mina = min(vala + lena, prec)
   minb = min(valb + lenb, prec)
   lenz = max(mina, minb) - valz
   R = base_ring(a)
   z = parent(a)()
   fit!(z, lenz)
   set_prec!(z, prec)
   set_val!(z, valz)
   if vala >= valb
      for i = 1:min(lenb, vala - valb)
         z = setcoeff!(z, i - 1, -polcoeff(b, i - 1))
      end
      for i = lenb + 1:min(vala - valb, lenz)
         z = setcoeff!(z, i - 1, R())
      end
      for i = vala - valb + 1:lenb
         z = setcoeff!(z, i - 1, polcoeff(a, i - vala + valb - 1) - polcoeff(b, i - 1))
      end
      for i = max(lenb, vala - valb) + 1:lena + vala - valb
         z = setcoeff!(z, i - 1, polcoeff(a, i - vala + valb - 1))
      end
      for i = lena + vala - valb + 1:lenb
         z = setcoeff!(z, i - 1, -polcoeff(b, i - 1))
      end
   else
      for i = 1:min(lena, valb - vala)
         z = setcoeff!(z, i - 1, polcoeff(a, i - 1))
      end
      for i = lena + 1:min(valb - vala, lenz)
         z = setcoeff!(z, i - 1, R())
      end
      for i = valb - vala + 1:lena
         z = setcoeff!(z, i - 1, polcoeff(a, i - 1) - polcoeff(b, i - valb + vala - 1))
      end
      for i = max(lena, valb - vala) + 1:lenb + valb - vala
         z = setcoeff!(z, i - 1, -polcoeff(b, i - valb + vala - 1))
      end
      for i = lenb + valb - vala + 1:lena
         z = setcoeff!(z, i - 1, polcoeff(a, i - 1))
      end
   end
   set_length!(z, normalise(z, lenz))
   renormalize!(z)
   return z
end

doc"""
    *{T <: RingElement}(a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T})
> Return $a\times b$.
"""
function *(a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T}) where {T <: RingElement}
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   zval = aval + bval
   prec = min(precision(a) - aval, precision(b) - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   if lena == 0 || lenb == 0
      return parent(a)(Array{T}(0), 0, prec + zval, zval)
   end
   t = base_ring(a)()
   lenz = min(lena + lenb - 1, prec)
   d = Array{T}(lenz)
   for i = 1:min(lena, lenz)
      d[i] = polcoeff(a, i - 1)*polcoeff(b, 0)
   end
   if lenz > lena
      for j = 2:min(lenb, lenz - lena + 1)
          d[lena + j - 1] = polcoeff(a, lena - 1)*polcoeff(b, j - 1)
      end
   end
   for i = 1:lena - 1
      if lenz > i
         for j = 2:min(lenb, lenz - i + 1)
            t = mul!(t, polcoeff(a, i - 1), polcoeff(b, j - 1))
            d[i + j - 1] = addeq!(d[i + j - 1], t)
         end
      end
   end
   z = parent(a)(d, lenz, prec + zval, zval)
   set_length!(z, normalise(z, lenz))
   renormalize!(z)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

doc"""
    *{T <: RingElem}(a::T, b::LaurentSeriesElem{T})
> Return $a\times b$.
"""
function *(a::T, b::LaurentSeriesElem{T}) where {T <: RingElem}
   len = pol_length(b)
   z = parent(b)()
   fit!(z, len)
   set_prec!(z, precision(b))
   set_val!(z, valuation(b))
   for i = 1:len
      z = setcoeff!(z, i - 1, a*polcoeff(b, i - 1))
   end
   set_length!(z, normalise(z, len))
   renormalize!(z)
   return z
end

doc"""
    *(a::Union{Integer, Rational, AbstractFloat}, b::LaurentSeriesElem)
> Return $a\times b$.
"""
function *(a::Union{Integer, Rational, AbstractFloat}, b::LaurentSeriesElem)
   len = pol_length(b)
   z = parent(b)()
   fit!(z, len)
   set_prec!(z, precision(b))
   set_val!(z, valuation(b))
   for i = 1:len
      z = setcoeff!(z, i - 1, a*polcoeff(b, i - 1))
   end
   set_length!(z, normalise(z, len))
   renormalize!(z)
   return z
end

doc"""
    *{T <: RingElem}(a::LaurentSeriesElem{T}, b::T)
> Return $a\times b$.
"""
*(a::LaurentSeriesElem{T}, b::T) where {T <: RingElem} = b*a

doc"""
    *(a::LaurentSeriesElem, b::Union{Integer, Rational, AbstractFloat})
> Return $a\times b$.
"""
*(a::LaurentSeriesElem, b::Union{Integer, Rational, AbstractFloat}) = b*a

###############################################################################
#
#   Shifting
#
###############################################################################

doc"""
    shift_left(x::LaurentSeriesElem, n::Int)
> Return the power series $f$ shifted left by $n$ terms, i.e. multiplied by
> $x^n$.
"""
function shift_left(x::LaurentSeriesElem{T}, len::Int) where {T <: RingElement}
   z = deepcopy(x)
   set_prec!(z, precision(x) + len)
   set_val!(z, valuation(x) + len)
   return z
end

doc"""
    shift_right(f::LaurentSeriesElem, n::Int)
> Return the power series $f$ shifted right by $n$ terms, i.e. divided by
> $x^n$.
"""
function shift_right(x::LaurentSeriesElem{T}, len::Int) where {T <: RingElement}
   z = deepcopy(x)
   set_prec!(z, precision(x) - len)
   set_val!(z, valuation(x) - len)
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

doc"""
    truncate(a::LaurentSeriesElem, n::Int)
> Return $a$ truncated to (absolute) precision $n$.
"""
function truncate(a::LaurentSeriesElem{T}, prec::Int) where {T <: RingElement}
   alen = pol_length(a)
   aprec = precision(a)
   aval = valuation(a)
   if aprec <= prec
      return a
   end
   z = parent(a)()
   set_prec!(z, prec)
   if prec <= aval
      set_length!(z, 0)
      set_val!(z, prec)
   else
      fit!(z, prec - aval)
      set_val!(z, aval)
      for i = 1:min(prec - aval, alen)
         z = setcoeff!(z, i - 1, polcoeff(a, i - 1))
      end
      set_length!(z, normalise(z, prec - aval))
   end
   return z
end

# Intended only for internal use, does not renormalize, assumes n >= 0
# Only efficient if valuation(a) == valuation(b) == 0
function mullow(a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T}, n::Int) where {T <: RingElement}
   lena = pol_length(a)
   lenb = pol_length(b)
   if lena == 0 || lenb == 0
      return zero(parent(a))
   end
   prec = min(precision(a), precision(b))
   t = base_ring(a)()
   lenz = min(lena + lenb - 1, n)
   d = Array{T}(lenz)
   for i = 1:min(lena, lenz)
      d[i] = coeff(a, i - 1)*coeff(b, 0)
   end
   if lenz > lena
      for j = 2:min(lenb, lenz - lena + 1)
          d[lena + j - 1] = coeff(a, lena - 1)*coeff(b, j - 1)
      end
   end
   for i = 1:lena - 1
      if lenz > i
         for j = 2:min(lenb, lenz - i + 1)
            t = mul!(t, coeff(a, i - 1), coeff(b, j - 1))
            d[i + j - 1] = addeq!(d[i + j - 1], t)
         end
      end
   end
   z = parent(a)(d, lenz, prec, 0)
   set_length!(z, normalise(z, lenz))
   return z
end

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

doc"""
    ^{T <: RingElement}(a::LaurentSeriesElem{T}, b::Int)
> Return $a^b$. We require $b \geq 0$.
"""
function ^(a::LaurentSeriesElem{T}, b::Int) where {T <: RingElement}
   # special case powers of x for constructing power series efficiently
   if pol_length(a) == 0
      z = parent(a)()
      set_prec!(z, b*valuation(a))
      set_val!(z, b*valuation(a))
      return z
   elseif b == 0
      # in fact, the result would be exact 1 if we had exact series
      z = one(parent(a))
      return z
   elseif isgen(a)
      z = parent(a)()
      fit!(z, 1)
      set_prec!(z, b + precision(a) - 1)
      set_val!(z, b)
      z = setcoeff!(z, 0, polcoeff(a, 0))
      set_length!(z, 1)
      return z
   elseif pol_length(a) == 1
      z = parent(a)(polcoeff(a, 0)^b)
      set_prec!(z, (b - 1)*valuation(a) + precision(a))
      set_val!(z, b*valuation(a))
      return z
   elseif b == 1
      return deepcopy(a)
   elseif b == -1
      return inv(a)
   end

   if b < 0
      a = inv(a)
      b = -b
   end

   bit = ~((~UInt(0)) >> 1)
   while (UInt(bit) & b) == 0
      bit >>= 1
   end
   val = valuation(a)
   a = shift_right(a, val)
   prec = precision(a)
   z = a
   bit >>= 1
   while bit !=0
      z = mullow(z, z, prec)
      if (UInt(bit) & b) != 0
         z = mullow(z, a, prec)
      end
      bit >>= 1
   end
   set_val!(z, b*val)
   set_prec!(z, b*val + prec)
   renormalize!(z)
   return z
end

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

doc"""
    =={T <: RingElement}(x::LaurentSeriesElem{T}, y::LaurentSeriesElem{T})
> Return `true` if $x == y$ arithmetically, otherwise return `false`. Recall
> that power series to different precisions may still be arithmetically
> equal to the minimum of the two precisions.
"""
function ==(x::LaurentSeriesElem{T}, y::LaurentSeriesElem{T}) where {T <: RingElement}
   check_parent(x, y)
   xval = valuation(x)
   xprec = precision(x)
   yval = valuation(y)
   yprec = precision(y)
   prec = min(xprec, yprec)
   if prec <= xval && prec <= yval
      return true
   end
   if xval != yval
      return false
   end
   xlen = normalise(x, min(pol_length(x), prec - xval))
   ylen = normalise(y, min(pol_length(y), prec - yval))
   if xlen != ylen
      return false
   end
   for i = 1:xlen
      if polcoeff(x, i - 1) != polcoeff(y, i - 1)
         return false
      end
   end
   return true
end

doc"""
    isequal{T <: RingElement}(x::LaurentSeriesElem{T}, y::LaurentSeriesElem{T})
> Return `true` if $x == y$ exactly, otherwise return `false`. Only if the
> power series are precisely the same, to the same precision, are they declared
> equal by this function.
"""
function isequal(x::LaurentSeriesElem{T}, y::LaurentSeriesElem{T}) where {T <: RingElement}
   if parent(x) != parent(y)
      return false
   end
   if precision(x) != precision(y) || pol_length(x) != pol_length(y) ||
      valuation(x) != valuation(y)
      return false
   end
   for i = 1:pol_length(x)
      if !isequal(polcoeff(x, i - 1), polcoeff(y, i - 1))
         return false
      end
   end
   return true
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

doc"""
    =={T <: RingElem}(x::LaurentSeriesElem{T}, y::T)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::LaurentSeriesElem{T}, y::T) where {T <: RingElem} = precision(x) == 0 ||
           ((pol_length(x) == 0 && iszero(y)) || (pol_length(x) == 1 &&
             valuation(x) == 0 && polcoeff(x, 0) == y))

doc"""
    =={T <: RingElem}(x::T, y::LaurentSeriesElem{T})
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::T, y::LaurentSeriesElem{T}) where {T <: RingElem} = y == x

doc"""
    ==(x::LaurentSeriesElem, y::Union{Integer, Rational, AbstractFloat})
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::LaurentSeriesElem, y::Union{Integer, Rational, AbstractFloat}) = precision(x) == 0 ||
                  ((pol_length(x) == 0 && iszero(y)) || (pol_length(x) == 1 &&
                    valuation(x) == 0 && polcoeff(x, 0) == y))

doc"""
    ==(x::Union{Integer, Rational, AbstractFloat}, y::LaurentSeriesElem)
> Return `true` if $x == y$ arithmetically, otherwise return `false`.
"""
==(x::Union{Integer, Rational, AbstractFloat}, y::LaurentSeriesElem) = y == x

###############################################################################
#
#   Approximation
#
###############################################################################

function Base.isapprox(f::LaurentSeriesElem, g::LaurentSeriesElem; atol::Real=sqrt(eps()))
   check_parent(f, g)
   nmin = min(precision(f), precision(g))
   i = 1
   while i <= nmin
      if !isapprox(coeff(f, i - 1), coeff(g, i - 1); atol=atol)
         return false
      end
      i += 1
   end
   return true
end

###############################################################################
#
#   Exact division
#
###############################################################################

doc"""
    divexact{T <: RingElement}(a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T})
> Return $a/b$. Requires $b$ to be invertible.
"""
function divexact(x::LaurentSeriesElem{T}, y::LaurentSeriesElem{T}) where {T <: RingElement}
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   v2 = valuation(y)
   if v2 != 0
      x = shift_right(x, v2)
      y = shift_right(y, v2)
   end
   y = truncate(y, precision(x) - valuation(x) + valuation(y))
   return x*inv(y)
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

doc"""
    divexact(a::LaurentSeriesElem, b::Union{Integer, Rational, AbstractFloat})
> Return $a/b$ where the quotient is expected to be exact.
"""
function divexact(x::LaurentSeriesElem, y::Union{Integer, Rational, AbstractFloat})
   y == 0 && throw(DivideError())
   lenx = pol_length(x)
   z = parent(x)()
   fit!(z, lenx)
   set_prec!(z, precision(x))
   set_val!(z, valuation(x))
   for i = 1:lenx
      z = setcoeff!(z, i - 1, divexact(polcoeff(x, i - 1), y))
   end
   return z
end

doc"""
    divexact{T <: RingElem}(a::LaurentSeriesElem{T}, b::T)
> Return $a/b$ where the quotient is expected to be exact.
"""
function divexact(x::LaurentSeriesElem{T}, y::T) where {T <: RingElem}
   iszero(y) && throw(DivideError())
   lenx = pol_length(x)
   z = parent(x)()
   fit!(z, lenx)
   set_prec!(z, precision(x))
   set_val!(z, valuation(x))
   for i = 1:lenx
      z = setcoeff!(z, i - 1, divexact(polcoeff(x, i - 1), y))
   end
   return z
end

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

doc"""
   inv(a::LaurentSeriesElem)
> Return the inverse of the power series $a$, i.e. $1/a$.
"""
function inv(a::LaurentSeriesElem)
   iszero(a) && throw(DivideError())
   a1 = polcoeff(a, 0)
   ainv = parent(a)()
   fit!(ainv, precision(a) - valuation(a))
   set_prec!(ainv, precision(a) - 2*valuation(a))
   set_val!(ainv, -valuation(a))
   !isunit(a1) && error("Unable to invert power series")
   if precision(a) - valuation(a) != 0
      ainv = setcoeff!(ainv, 0, divexact(one(base_ring(a)), a1))
   end
   a1 = -a1
   for n = 2:precision(a) - valuation(a)
      s = polcoeff(a, 1)*polcoeff(ainv, n - 2)
      for i = 2:min(n, pol_length(a)) - 1
         s += polcoeff(a, i)*polcoeff(ainv, n - i - 1)
      end
      ainv = setcoeff!(ainv, n - 1, divexact(s, a1))
   end
   set_length!(ainv, normalise(ainv, precision(a) - valuation(a)))
   return ainv
end

###############################################################################
#
#   Special functions
#
###############################################################################

doc"""
    exp(a::LaurentSeriesElem)
> Return the exponential of the power series $a$.
"""
function Base.exp(a::LaurentSeriesElem)
   if iszero(a)
      z = one(parent(a))
      set_prec!(z, precision(a))
      return z
   end
   vala = valuation(a)
   preca = precision(a)
   vala < 0 && error("Valuation must be non-negative in exp")
   z = parent(a)()
   R = base_ring(a)
   fit!(z, preca)
   set_prec!(z, preca)
   set_val!(z, 0)
   c = vala == 0 ? polcoeff(a, 0) : R()
   z = setcoeff!(z, 0, Nemo.exp(c))
   len = pol_length(a) + vala
   for k = 1 : preca - 1
      s = R()
      for j = 1 : min(k + 1, len) - 1
         c = j >= vala ? polcoeff(a, j - vala) : R()
         s += j * c * polcoeff(z, k - j)
      end
      !isunit(R(k)) && error("Unable to divide in exp")
      z = setcoeff!(z, k, divexact(s, k))
   end
   set_length!(z, normalise(z, preca))
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(a::LaurentSeriesElem)
   a.length = 0
   a.prec = parent(a).prec_max
   return a
end

function fit!(c::LaurentSeriesElem{T}, n::Int) where {T <: RingElement}
   if length(c.coeffs) < n
      t = c.coeffs
      c.coeffs = Array{T}(n)
      for i = 1:c.length
         c.coeffs[i] = t[i]
      end
      for i = pol_length(c) + 1:n
         c.coeffs[i] = zero(base_ring(c))
      end
   end
   return nothing
end

function setcoeff!(c::LaurentSeriesElem{T}, n::Int, a::T) where {T <: RingElement}
   if (a != 0 && precision(c) - valuation(c) > n) || n + 1 <= c.length
      fit!(c, n + 1)
      c.coeffs[n + 1] = a
      c.length = max(pol_length(c), n + 1)
      # don't normalise
   end
   return c
end

function mul!(c::LaurentSeriesElem{T}, a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T}) where {T <: RingElement}
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(precision(a) - aval, precision(b) - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   if lena <= 0 || lenb <= 0
      c.length = 0
   else
      t = base_ring(a)()
      lenc = min(lena + lenb - 1, prec)
      fit!(c, lenc)
      for i = 1:min(lena, lenc)
         c.coeffs[i] = mul!(c.coeffs[i], polcoeff(a, i - 1), polcoeff(b, 0))
      end
      if lenc > lena
         for i = 2:min(lenb, lenc - lena + 1)
            c.coeffs[lena + i - 1] = mul!(c.coeffs[lena + i - 1], polcoeff(a, lena - 1), polcoeff(b, i - 1))
         end
      end
      for i = 1:lena - 1
         if lenc > i
            for j = 2:min(lenb, lenc - i + 1)
               t = mul!(t, polcoeff(a, i - 1), polcoeff(b, j - 1))
               c.coeffs[i + j - 1] = addeq!(c.coeffs[i + j - 1], t)
            end
         end
      end
      c.length = normalise(c, lenc)
   end
   c.val = a.val + b.val
   c.prec = prec + c.val
   renormalize!(c)
   return c
end

function addeq!(c::LaurentSeriesElem{T}, a::LaurentSeriesElem{T}) where {T <: RingElement}
   lenc = pol_length(c)
   lena = pol_length(a)
   valc = valuation(c)
   vala = valuation(a)
   valr = min(vala, valc)
   precc = precision(c)
   preca = precision(a)
   prec = min(precc, preca)
   mina = min(vala + lena, prec)
   minc = min(valc + lenc, prec)
   lenr = max(mina, minc) - valr
   R = base_ring(c)
   fit!(c, lenr)
   if valc >= vala
      for i = lenc + valc - vala:-1:max(lena, valc - vala) + 1
         t = c.coeffs[i]
         c.coeffs[i] = c.coeffs[i - valc + vala]
         c.coeffs[i - valc + vala] = t
      end
      for i = lena:-1:valc - vala + 1
         c.coeffs[i] = add!(c.coeffs[i], c.coeffs[i - valc + vala], a.coeffs[i])
      end
      for i = 1:min(lena, valc - vala)
         c.coeffs[i] = a.coeffs[i]
      end
      for i = lena + 1:min(valc - vala, lenr)
         c.coeffs[i] = R()
      end
      for i = lenc + valc - vala + 1:lena
         c.coeffs[i] = a.coeffs[i]
      end
   else
      for i = lenc + 1:min(vala - valc, lenr)
         c.coeffs[i] = R()
      end
      for i = vala - valc + 1:lenc
         c.coeffs[i] = addeq!(c.coeffs[i], a.coeffs[i - vala + valc])
      end
      for i = max(lenc, vala - valc) + 1:lena + vala - valc
         c.coeffs[i] = a.coeffs[i - vala + valc]
      end
   end
   c.length = normalise(c, lenr)
   c.prec = prec
   c.val = valr
   renormalize!(c)
   return c
end

function add!(c::LaurentSeriesElem{T}, a::LaurentSeriesElem{T}, b::LaurentSeriesElem{T}) where {T <: RingElement}
   if c === a
      return addeq!(c, b)
   elseif c === b
      return addeq!(c, a)
   end
   lena = pol_length(a)
   lenb = pol_length(b)
   valb = valuation(b)
   vala = valuation(a)
   valr = min(vala, valb)
   precb = precision(b)
   preca = precision(a)
   prec = min(precb, preca)
   mina = min(vala + lena, prec)
   minb = min(valb + lenb, prec)
   lenr = max(mina, minb) - valr
   R = base_ring(c)
   fit!(c, lenr)
   c.prec = prec
   c.val = valr
   if vala > valb
      for i = 1:min(lenb, vala - valb)
         c.coeffs[i] = b.coeffs[i]
      end
      for i = lenb + 1:vala - valb
         c.coeffs[i] = R()
      end
      for i = vala - valb + 1:lenb
         c.coeffs[i] = add!(c.coeffs[i], a.coeffs[i - vala + valb], b.coeffs[i])
      end
      for i = max(lenb, vala - valb) + 1:lena + vala - valb
         c.coeffs[i] = a.coeffs[i - vala + valb]
      end
      for i = lena + vala - valb + 1:lenb
         c.coeffs[i] = b.coeffs[i]
      end
   else
      for i = 1:min(lena, valb - vala)
         c.coeffs[i] = a.coeffs[i]
      end
      for i = lena + 1:valb - vala
         c.coeffs[i] = R()
      end
      for i = valb - vala + 1:lena
         c.coeffs[i] = add!(c.coeffs[i], a.coeffs[i], b.coeffs[i - valb + vala])
      end
      for i = max(lena, valb - vala) + 1:lenb + valb - vala
         c.coeffs[i] = b.coeffs[i - valb + vala]
      end
      for i = lenb + valb - vala + 1:lena
         c.coeffs[i] = a.coeffs[i]
      end
   end
   set_length!(c, normalise(c, lenr))
   renormalize!(c)
   return c
end

###############################################################################
#
#   Random elements
#
###############################################################################

function rand(S::LaurentSeriesRing, val_range::UnitRange{Int}, v...)
   R = base_ring(S)
   f = S()
   x = gen(S)
   for i = 0:S.prec_max - 1
      f += rand(R, v...)*x^i
   end
   return shift_left(f, rand(val_range))
end

function rand(S::LaurentSeriesField, val_range::UnitRange{Int}, v...)
   R = base_ring(S)
   f = S()
   x = gen(S)
   for i = 0:S.prec_max - 1
      f += rand(R, v...)*x^i
   end
   return shift_left(f, rand(val_range))
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{LaurentSeriesElem{T}}, ::Type{LaurentSeriesElem{T}}) where T <: RingElement = RelSeries{T}

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

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (R::LaurentSeriesRing{T})(b::RingElement) where {T <: RingElement}
   return R(base_ring(R)(b))
end

function (R::LaurentSeriesField{T})(b::RingElement) where {T <: FieldElement}
   return R(base_ring(R)(b))
end

function (R::LaurentSeriesRing{T})() where {T <: RingElement}
   z = LaurentSeriesRingElem{T}(Array{T}(0), 0, R.prec_max, R.prec_max)
   z.parent = R
   return z
end

function (R::LaurentSeriesField{T})() where {T <: FieldElement}
   z = LaurentSeriesFieldElem{T}(Array{T}(0), 0, R.prec_max, R.prec_max)
   z.parent = R
   return z
end

function (R::LaurentSeriesRing{T})(b::Union{Integer, Rational, AbstractFloat}) where {T <: RingElement}
   if b == 0
      z = LaurentSeriesRingElem{T}(Array{T}(0), 0, R.prec_max, R.prec_max)
   else
      z = LaurentSeriesRingElem{T}([base_ring(R)(b)], 1, R.prec_max, 0)
   end
   z.parent = R
   return z
end

function (R::LaurentSeriesField{T})(b::Union{Rational, AbstractFloat}) where {T <: FieldElement}
   if b == 0
      z = LaurentSeriesFieldElem{T}(Array{T}(0), 0, R.prec_max, R.prec_max)
   else
      z = LaurentSeriesFieldElem{T}([base_ring(R)(b)], 1, R.prec_max, 0)
   end
   z.parent = R
   return z
end

function (R::LaurentSeriesRing{T})(b::T) where {T <: RingElem}
   parent(b) != base_ring(R) && error("Unable to coerce to power series")
   if iszero(b)
      z = LaurentSeriesRingElem{T}(Array{T}(0), 0, R.prec_max, R.prec_max)
   else
      z = LaurentSeriesRingElem{T}([b], 1, R.prec_max, 0)
   end
   z.parent = R
   return z
end

function (R::LaurentSeriesField{T})(b::T) where {T <: FieldElem}
   parent(b) != base_ring(R) && error("Unable to coerce to power series")
   if iszero(b)
      z = LaurentSeriesFieldElem{T}(Array{T}(0), 0, R.prec_max, R.prec_max)
   else
      z = LaurentSeriesFieldElem{T}([b], 1, R.prec_max, 0)
   end
   z.parent = R
   return z
end

function (R::LaurentSeriesRing{T})(b::LaurentSeriesElem{T}) where {T <: RingElement}
   parent(b) != R && error("Unable to coerce power series")
   return b
end

function (R::LaurentSeriesField{T})(b::LaurentSeriesElem{T}) where {T <: FieldElement}
   parent(b) != R && error("Unable to coerce power series")
   return b
end

function (R::LaurentSeriesRing{T})(b::Array{T, 1}, len::Int, prec::Int, val::Int) where {T <: RingElement}
   if length(b) > 0
      parent(b[1]) != base_ring(R) && error("Unable to coerce to power series")
   end
   z = LaurentSeriesRingElem{T}(b, len, prec, val)
   z.parent = R
   return z
end

function (R::LaurentSeriesField{T})(b::Array{T, 1}, len::Int, prec::Int, val::Int) where {T <: RingElement}
   if length(b) > 0
      parent(b[1]) != base_ring(R) && error("Unable to coerce to power series")
   end
   z = LaurentSeriesFieldElem{T}(b, len, prec, val)
   z.parent = R
   return z
end

###############################################################################
#
#   PowerSeriesRing constructor
#
###############################################################################

doc"""
   PowerSeriesRing(R::Nemo.Ring, prec::Int, s::AbstractString; cached=true, model=:capped_relative)
> Return a tuple $(S, x)$ consisting of the parent object `S` of a power series
> ring over the given base ring and a generator `x` for the power series ring.
> The maximum precision of power series in the ring is set to `prec`. If the
> model is set to `:capped_relative` this is taken as a maximum relative
> precision, and if it is set to `:capped_absolute` this is take to be a
> maximum absolute precision. The supplied string `s` specifies the way the
> generator of the power series ring will be printed. By default, the parent
> object `S` will be cached so that supplying the same base ring, string and
> precision in future will return the same parent object and generator. If
> caching of the parent object is not required, `cached` can be set to `false`.
"""
function LaurentSeriesRing(R::Nemo.Ring, prec::Int, s::AbstractString; cached=true)
   S = Symbol(s)
   T = elem_type(R)

   parent_obj = LaurentSeriesRing{T}(R, prec, S, cached)

   return parent_obj, gen(parent_obj)
end

function LaurentSeriesRing(R::Nemo.Field, prec::Int, s::AbstractString; cached=true)
   S = Symbol(s)
   T = elem_type(R)

   parent_obj = LaurentSeriesField{T}(R, prec, S, cached)

   return parent_obj, gen(parent_obj)
end

function LaurentSeriesField(R::Nemo.Field, prec::Int, s::AbstractString; cached=true)
   S = Symbol(s)
   T = elem_type(R)

   parent_obj = LaurentSeriesField{T}(R, prec, S, cached)

   return parent_obj, gen(parent_obj)
end


text:
function setup_etc(backend::GPUBackend, layer::DropoutLayer, inputs::Vector{Blob})
  cuda_rand_states = CuPtr
  kernel = backend.mocha.dropout_init
  rnd_state_size_blob = make_blob(backend, Float64, 1, 1, 1, 1)
  CUDA.launch(backend.mocha.dropout_alloc_size, 1, 1, (rnd_state_size_blob.ptr.p, ))
  rnd_state_size = Float64[0]
  copy!(rnd_state_size, rnd_state_size_blob)
  destroy(rnd_state_size_blob)
  rnd_state_size = round(Int, rnd_state_size[1])

  len = length(inputs[1])
  cuda_rand_states = CUDA.cualloc(UInt8, rnd_state_size*len)
  x_block = round(Int, ceil(convert(Float64, len)/CUDA.THREADS_PER_BLOCK_X))
  CUDA.launch(kernel, x_block, CUDA.THREADS_PER_BLOCK_X, (cuda_rand_states, len))

  # hold copy of input blob, we will restore the inputs after backward computing
  # this is because when used with the cuDNN pooling layer, the in-place modification
  # of the pooling layer output will cause error in computed gradient due to unknown
  # implementation of the pooling backward. So as a workaround, we restore the input
  # after finishing dropout backward computing. See #96 on github for more details.
  input_copy = make_blob(backend, eltype(inputs[1]), size(inputs[1]))

  return (cuda_rand_states, input_copy)
end

function destroy_etc(backend::GPUBackend, state::DropoutLayerState)
  cuda_rand_states, input_copy = state.etc
  CUDA.free(cuda_rand_states)
  destroy(input_copy)
end

function forward(backend::GPUBackend, state::DropoutLayerState, inputs::Vector{Blob})
  # make copy of input blob
  copy!(state.etc[2], inputs[1])

  len = length(inputs[1])
  x_block = round(Int, ceil(convert(Float64, len)/CUDA.THREADS_PER_BLOCK_X))
  data_type = eltype(inputs[1])
  if data_type == Float32
    kernel = backend.mocha.dropout_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.dropout_forward_double
  end

  CUDA.launch(kernel, x_block, CUDA.THREADS_PER_BLOCK_X,
      (state.etc[1], length(inputs[1]), inputs[1].ptr.p,
      state.rand_vals.ptr.p, state.ratio, state.scale))
end

function backward(backend::GPUBackend, state::DropoutLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  if !isa(diffs[1], NullBlob)
    len = length(inputs[1])
    x_block = round(Int, ceil(convert(Float64, len)/CUDA.THREADS_PER_BLOCK_X))
    data_type = eltype(inputs[1])
    if data_type == Float32
      kernel = backend.mocha.dropout_backward_float
    elseif data_type == Float64
      kernel = backend.mocha.dropout_backward_double
    end

    CUDA.launch(kernel, x_block, CUDA.THREADS_PER_BLOCK_X,
        (state.etc[1], length(inputs[1]), diffs[1].ptr.p,
        state.rand_vals.ptr.p, state.ratio, state.scale))

    # restore the input blob
    copy!(inputs[1], state.etc[2])
  end
end

ext: jl


text:
using JImages
using Images
using Test

function test_image_fom(x::AbstractArray, y::AbstractArray)
    dff = x .- y
    adff = broadcast(abs,dff)
    sdff = sum(adff) /sum(x)
    return sdff
end

function test_limit_boundaries_low(y::AbstractArray)
    yc = copy(y)
    yc2 = copy(y)
    for i in 1:5
        yc[i,i] = 0
        yc2[i,i] = 0.5
    end
    JImages.limit_boundaries(yc, 0.5)
    return yc, yc2
end

function test_zero_protection(y::AbstractArray)
    yc = copy(y)
    yc2 = copy(y)
    for i in 1:5
        yc[i,i] = 0
        yc2[i,i] = 0.5
    end
    println("in test zero protection")
    JImages.zero_protection(yc, 0.5)
    return yc, yc2
end

function test_limit_boundaries_high(y::AbstractArray)
    yc = copy(y)
    yc2 = copy(y)
    for i in 1:5
        yc[i,i] = 200.
        yc2[i,i] = 100.
    end
    JImages.limit_boundaries(yc, eps(), 100.)
    return yc,yc2
end
@testset "JImages.jl" begin
    x = reshape([i+(j-1)* 1. * i for j in 1:5 for i in 1:5],(5,5))
    y = reshape([i+(j-1)*2. * i for j in 1:5 for i in 1:5],(5,5))
    xy = reshape([(i,j) for j in -2:2 for i in -2:2], (5,5))
    xc = reshape([[4.,6.,8.] [6.,9.,12.] [8.,12.,16.]],(3,3))
    xf = reshape([[25,
                20, 15, 10, 5] [20,
                16, 12, 8, 4] [15,
                12, 9, 6, 3] [10,
                8, 6, 4, 2] [5,
                4, 3, 2, 1]],
               (5,5)) .*1.0
    kg1 = Kernel.gaussian(1)
    ktest = JImages.gaussianxy(xy)
    @test JImages.dxy(10, -3) == 13
    @test JImages.image_fom(x, y) == test_image_fom(x, y)
    @test JImages.center_view(x, 1) == xc
    @test JImages.flip2d(x) == xf
    @test isapprox(ktest, collect(kg1))

    yl, yl2 = test_limit_boundaries_low(y)
    yh, yh2 = test_limit_boundaries_high(y)
    yz, yz2 = test_zero_protection(y)
    @test isapprox(yl,yl2)
    @test isapprox(yh,yh2)
    @test isapprox(yz,yz2)
end

max_stars_repo_name: jjgomezcadenas/JImages

max_issues_repo_name: jjgomezcadenas/JImages

max_forks_repo_name: jjgomezcadenas/JImages


text:
using ChanVese
using Test
using TestSetExtensions


text:
include("connectnodes.jl")
include("helixnodes.jl")
include("viagroup.jl")
include("coilcraft1010vsgroup.jl")


text:
"""
    collect_dependent_system_nodes!(dag::ExprDAG, node::ExprNode, system_nodes)
Collect dependent nodes in the system headers.
"""
function collect_dependent_system_nodes! end

collect_dependent_system_nodes!(dag::ExprDAG, node::ExprNode, system_nodes) = dag

function collect_dependent_system_nodes!(dag::ExprDAG, node::ExprNode{FunctionProto}, system_nodes)
    cursor = node.cursor
    types = CLType[getArgType(getCursorType(cursor), i - 1) for i in 1:getNumArguments(cursor)]
    push!(types, getCursorResultType(cursor))
    for ty in types
        # for arrays and pointers, meaningful types need to be recursively retrieved
        jlty = tojulia(ty)
        leaf_ty = get_jl_leaf_type(jlty)

        # skip julia basic types
        is_jl_basic(leaf_ty) && continue

        hasref = has_elaborated_reference(ty)
        if (hasref && haskey(dag.tags, leaf_ty.sym)) ||
            (!hasref && haskey(dag.ids, leaf_ty.sym)) ||
            haskey(dag.ids_extra, leaf_ty.sym) ||
            (hasref && is_jl_pointer(jlty))
            # pass
        else
            # add all sys nodes with the same id
            for (i, n) in enumerate(dag.sys)
                if n.id == leaf_ty.sym
                    system_nodes[n] = i
                end
            end
        end
    end
    return system_nodes
end

function collect_dependent_system_nodes!(dag::ExprDAG, node::ExprNode{FunctionNoProto}, system_nodes)
    cursor = node.cursor
    @assert getNumArguments(cursor) == 0
    ty = getCursorResultType(cursor)
    jlty = tojulia(ty)
    leaf_ty = get_jl_leaf_type(jlty)

    is_jl_basic(leaf_ty) && return system_nodes

    hasref = has_elaborated_reference(ty)
    if (hasref && haskey(dag.tags, leaf_ty.sym)) ||
        (!hasref && haskey(dag.ids, leaf_ty.sym)) ||
        haskey(dag.ids_extra, leaf_ty.sym)
        # pass
    else
        # add all sys nodes with the same id
        for (i, n) in enumerate(dag.sys)
            if n.id == leaf_ty.sym
                system_nodes[n] = i
            end
        end
    end

    return system_nodes
end

function collect_dependent_system_nodes!(dag::ExprDAG, node::ExprNode{TypedefElaborated}, system_nodes)
    cursor = node.cursor
    ty = getTypedefDeclUnderlyingType(cursor)
    jlty = tojulia(ty)
    leaf_ty = get_jl_leaf_type(jlty)

    is_jl_basic(leaf_ty) && return system_nodes

    if !(haskey(dag.tags, leaf_ty.sym) || haskey(dag.ids_extra, leaf_ty.sym))
        # add all sys nodes with the same id
        for (i, n) in enumerate(dag.sys)
            if n.id == leaf_ty.sym
                system_nodes[n] = i
            end
        end
    end

    return system_nodes
end

function collect_dependent_system_nodes!(dag::ExprDAG, node::ExprNode{TypedefDefault}, system_nodes)
    cursor = node.cursor
    ty = getTypedefDeclUnderlyingType(cursor)
    jlty = tojulia(ty)
    leaf_ty = get_jl_leaf_type(jlty)

    is_jl_basic(leaf_ty) && return system_nodes

    # do nothing for unknowns since we just skip them in the downstream passes
    is_jl_unknown(leaf_ty) && return system_nodes

    if !(haskey(dag.ids, leaf_ty.sym) || haskey(dag.ids_extra, leaf_ty.sym))
        # add all sys nodes with the same id
        for (i, n) in enumerate(dag.sys)
            if n.id == leaf_ty.sym
                system_nodes[n] = i
            end
        end
    end

    return system_nodes
end

function collect_dependent_system_nodes!(dag::ExprDAG, node::ExprNode{TypedefToAnonymous}, system_nodes)
    # add all sys nodes with the same id
    for (i, n) in enumerate(dag.sys)
        if n.id == node.type.sym
            system_nodes[n] = i
        end
    end
    return system_nodes
end

function collect_dependent_system_nodes!(dag::ExprDAG, node::ExprNode{<:Union{AbstractStructNodeType, <:AbstractUnionNodeType}}, system_nodes)
    collect_dependent_system_nodes!(dag, getCursorType(node.cursor), system_nodes)
end

function collect_dependent_system_nodes!(dag::ExprDAG, type::CLType, system_nodes)
    for c in fields(type)
        ty = getCursorType(c)
        jlty = tojulia(ty)
        leaf_ty = get_jl_leaf_type(jlty)

        is_jl_basic(leaf_ty) && continue

        hasref = has_elaborated_reference(ty)
        if (hasref && haskey(dag.tags, leaf_ty.sym)) ||
            (!hasref && haskey(dag.ids, leaf_ty.sym)) ||
            haskey(dag.ids_extra, leaf_ty.sym) ||
            occursin("anonymous", spelling(ty))
            # nested tags may also be from system headers
            collect_dependent_system_nodes!(dag, ty, system_nodes)
        else
            # add all sys nodes with the same id
            for (i, n) in enumerate(dag.sys)
                if n.id == leaf_ty.sym
                    system_nodes[n] = i
                end
            end
            if Clang.is_elaborated(ty)
                collect_dependent_system_nodes!(dag, ty, system_nodes)
            end
        end
    end
    return dag
end


text:
# bracketing intervals for RSSA-based aggregators
# see: "On the rejection-based algorithm for simulation and analysis of
#       large-scale reaction networks", Thanh et al, J. Chem. Phys., 2015
# note, expects the type of the bracketing variables [ulow,uhigh] to be the
# same as the fluct_rate and ushift.
struct BracketData{T1,T2}
    fluctrate::T1         # interval should be [1-fluctrate,1+fluctrate] * u
    threshold::T2         # for u below threshold interval is:
    Δu::T2                #   [max(u-Δu,0),u+Δu]
end

# # suggested defaults
# BracketData{T1}() = BracketData(.1, 25, 4)
BracketData{T1,T2}() where {T1,T2} = BracketData(T1(.1),T2(25),T2(4))

# support either vectors of data for each field, or scalars
@inline getfr(bd::BracketData{AbstractVector{T1},T2}, i) where {T1,T2} = bd.fluctrate[i]
@inline getfr(bd::BracketData{T1,T2}, i) where {T1 <: Number,T2} = bd.fluctrate

@inline gettv(bd::BracketData{T1,AbstractVector{T2}}, i) where {T1,T2} = bd.threshold[i]
@inline gettv(bd::BracketData{T1,T2}, i) where {T1,T2 <: Number} = bd.threshold

@inline getΔu(bd::BracketData{T1,AbstractVector{T2}}, i) where {T1,T2} = bd.Δu[i]
@inline getΔu(bd::BracketData{T1,T2}, i) where {T1,T2 <: Number} = bd.Δu

@inline function get_spec_brackets(bd, i, u)
    if u == zero(u)
        return zero(u), zero(u)
    elseif u < gettv(bd, i)
        Δu = getΔu(bd,i)
        return max(zero(Δu), u - Δu), u + Δu
    else
        δ = getfr(bd, i)
        return trunc(typeof(u), (one(δ) - δ) * u), trunc(typeof(u), (one(δ) + δ) * u)
    end
end

@inline get_majump_brackets(ulow, uhigh, k, majumps) = evalrxrate(ulow, k, majumps), evalrxrate(uhigh, k, majumps)

# for constant rate jumps we must check the ordering of the bracket values
@inline function get_cjump_brackets(ulow, uhigh, rate, params, t)
    rlow  = rate(ulow, params, t)
    rhigh = rate(uhigh, params, t)
    if rlow <= rhigh
        return rlow,rhigh
    else
        return rhigh,rlow
    end
end

"get brackets for reaction rx by first checking if the reaction is a massaction reaction"
@inline function get_jump_brackets(rx, p :: AbstractSSAJumpAggregator, params, t)
    ma_jumps = p.ma_jumps
    num_majumps = get_num_majumps(ma_jumps)
    if rx <= num_majumps
        return get_majump_brackets(p.ulow, p.uhigh, rx, ma_jumps)
    else
        @inbounds return get_cjump_brackets(p.ulow, p.uhigh, p.rates[rx - num_majumps], params, t)
    end
end

# set up bracketing
function set_bracketing!(p :: AbstractSSAJumpAggregator, u, params, t)
    # species bracketing interval
    ubnds = p.cur_u_bnds
    @inbounds for (i,uval) in enumerate(u)
        ubnds[1,i], ubnds[2,i] = get_spec_brackets(p.bracket_data, i, uval)
    end

    # reaction rate bracketing interval
    # mass action jumps
    sum_rate = zero(p.sum_rate)
    majumps  = p.ma_jumps
    crlow    = p.cur_rate_low
    crhigh   = p.cur_rate_high
    @inbounds for k = 1:get_num_majumps(majumps)
        crlow[k], crhigh[k] = get_majump_brackets(p.ulow, p.uhigh, k, majumps)
        sum_rate += crhigh[k]
    end

    # constant rate jumps
    k = get_num_majumps(majumps) + 1
    @inbounds for rate in p.rates
        crlow[k], crhigh[k] = get_cjump_brackets(p.ulow, p.uhigh, rate, params, t)
        sum_rate += crhigh[k]
        k += 1
    end
    p.sum_rate = sum_rate

    nothing
end


text:
using CLIMA.PlanetParameters
export BoundaryCondition, PeriodicBC, NoFluxBC, InitStateBC

function atmos_boundary_flux_diffusive!(nf::NumericalFluxDiffusive,
                                        bc,
                                        atmos::AtmosModel,
                                        F,
                                        state⁺, diff⁺, aux⁺, n⁻,
                                        state⁻, diff⁻, aux⁻,
                                        bctype, t, state1⁻, diff1⁻, aux1⁻)
  atmos_boundary_state!(nf, bc, atmos,
                        state⁺, diff⁺, aux⁺, n⁻,
                        state⁻, diff⁻, aux⁻,
                        bctype, t,
                        state1⁻, diff1⁻, aux1⁻)
  flux_diffusive!(atmos, F, state⁺, diff⁺, aux⁺, t)
end

#TODO: figure out a better interface for this.
# at the moment we can just pass a function, but we should do something better
# need to figure out how subcomponents will interact.
function atmos_boundary_state!(::Union{NumericalFluxNonDiffusive, NumericalFluxGradient},
                               f::Function, m::AtmosModel, state⁺::Vars,
                               aux⁺::Vars, n⁻, state⁻::Vars, aux⁻::Vars, bctype,
                               t, _...)
  f(state⁺, aux⁺, n⁻, state⁻, aux⁻, bctype, t)
end

function atmos_boundary_state!(::NumericalFluxDiffusive, f::Function,
                               m::AtmosModel, state⁺::Vars, diff⁺::Vars,
                               aux⁺::Vars, n⁻, state⁻::Vars, diff⁻::Vars,
                               aux⁻::Vars, bctype, t, _...)
  f(state⁺, diff⁺, aux⁺, n⁻, state⁻, diff⁻, aux⁻, bctype, t)
end

# lookup boundary condition by face
function atmos_boundary_state!(nf::Union{NumericalFluxNonDiffusive, NumericalFluxGradient},
                               bctup::Tuple, m::AtmosModel, state⁺::Vars,
                               aux⁺::Vars, n⁻, state⁻::Vars, aux⁻::Vars, bctype,
                               t, _...)
  atmos_boundary_state!(nf, bctup[bctype], m, state⁺, aux⁺, n⁻, state⁻, aux⁻,
                        bctype, t)
end

function atmos_boundary_state!(nf::NumericalFluxDiffusive,
                               bctup::Tuple, m::AtmosModel, state⁺::Vars,
                               diff⁺::Vars, aux⁺::Vars, n⁻, state⁻::Vars,
                               diff⁻::Vars, aux⁻::Vars, bctype, t, _...)
  atmos_boundary_state!(nf, bctup[bctype], m, state⁺, diff⁺, aux⁺, n⁻, state⁻,
                        diff⁻, aux⁻, bctype, t)
end

abstract type BoundaryCondition
end

"""
    PeriodicBC <: BoundaryCondition

Assume that the topology is periodic and hence nothing special needs to be done at the boundaries.
"""
struct PeriodicBC <: BoundaryCondition end

# TODO: assert somewhere that the topology is actually periodic when using those
atmos_boundary_state!(_, ::PeriodicBC, _...) = nothing

"""
    NoFluxBC <: BoundaryCondition

Set the momentum at the boundary to be zero.

# TODO: This should be fixed later once BCs are figured out (likely want
# different things here?)

"""
struct NoFluxBC <: BoundaryCondition
end

function atmos_boundary_state!(nf::Union{NumericalFluxNonDiffusive, NumericalFluxGradient},
                               bc::NoFluxBC, m::AtmosModel, state⁺::Vars,
                               aux⁺::Vars, n⁻, state⁻::Vars, aux⁻::Vars, bctype,
                               t, _...)
  FT = eltype(state⁻)
  state⁺.ρ = state⁻.ρ
  if typeof(nf) <: NumericalFluxNonDiffusive
    state⁺.ρu -= 2 * dot(state⁻.ρu, n⁻) * SVector(n⁻)
  else
    state⁺.ρu -=  dot(state⁻.ρu, n⁻) * SVector(n⁻)
  end
end

function atmos_boundary_state!(::NumericalFluxDiffusive, bc::NoFluxBC,
                               m::AtmosModel, state⁺::Vars, diff⁺::Vars,
                               aux⁺::Vars, n⁻, state⁻::Vars, diff⁻::Vars,
                               aux⁻::Vars, bctype, t, _...)
  FT = eltype(state⁻)
  state⁺.ρ = state⁻.ρ
  state⁺.ρu -= dot(state⁻.ρu, n⁻) * SVector(n⁻)

  fill!(getfield(diff⁺, :array), FT(0))
end

"""
    InitStateBC <: BoundaryCondition

Set the value at the boundary to match the `init_state!` function. This is
mainly useful for cases where the problem has an explicit solution.

# TODO: This should be fixed later once BCs are figured out (likely want
# different things here?)
"""
struct InitStateBC <: BoundaryCondition
end
function atmos_boundary_state!(::Union{NumericalFluxNonDiffusive, NumericalFluxGradient},
                               bc::InitStateBC, m::AtmosModel, state⁺::Vars,
                               aux⁺::Vars, n⁻, state⁻::Vars, aux⁻::Vars, bctype,
                               t, _...)
  init_state!(m, state⁺, aux⁺, aux⁺.coord, t)
end
function atmos_boundary_state!(::NumericalFluxDiffusive, bc::InitStateBC,
                               m::AtmosModel, state⁺::Vars, diff⁺::Vars,
                               aux⁺::Vars, n⁻, state⁻::Vars, diff⁻::Vars,
                               aux⁻::Vars, bctype, t, _...)
  init_state!(m, state⁺, aux⁺, aux⁺.coord, t)
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using ONNX, Flux, ProtoBuf
include("ops_tests.jl")

# Maxpool 1D default:
ip = read_input("$ONNX_TEST_PATH/test_maxpool_1d_default")
main_test("$ONNX_TEST_PATH/test_maxpool_1d_default",
          read_output("$ONNX_TEST_PATH/test_maxpool_1d_default"),
          read_input("$ONNX_TEST_PATH/test_maxpool_1d_default")[1])

# Maxpool 2D default
ip = read_input("$ONNX_TEST_PATH/test_maxpool_2d_default")
main_test("$ONNX_TEST_PATH/test_maxpool_2d_default",
          read_output("$ONNX_TEST_PATH/test_maxpool_2d_default"),
          read_input("$ONNX_TEST_PATH/test_maxpool_2d_default")[1])

# Maxpool 2D precomputed pads:
ip = read_input("$ONNX_TEST_PATH/test_maxpool_2d_precomputed_pads")
main_test("$ONNX_TEST_PATH/test_maxpool_2d_precomputed_pads",
          read_output("$ONNX_TEST_PATH/test_maxpool_2d_precomputed_pads"),
          read_input("$ONNX_TEST_PATH/test_maxpool_2d_precomputed_pads")[1])

# Maxpool 2D precomputed strides:
ip = read_input("$ONNX_TEST_PATH/test_maxpool_2d_precomputed_strides")
main_test("$ONNX_TEST_PATH/test_maxpool_2d_precomputed_strides",
          read_output("$ONNX_TEST_PATH/test_maxpool_2d_precomputed_strides"),
          read_input("$ONNX_TEST_PATH/test_maxpool_2d_precomputed_strides")[1])

# Maxpool 2D strides:
ip = read_input("$ONNX_TEST_PATH/test_maxpool_2d_strides")
main_test("$ONNX_TEST_PATH/test_maxpool_2d_strides",
          read_output("$ONNX_TEST_PATH/test_maxpool_2d_strides"),
          read_input("$ONNX_TEST_PATH/test_maxpool_2d_strides")[1])

# Averagepool 1D
ip = read_input("$ONNX_TEST_PATH/test_averagepool_1d_default")
    main_test("$ONNX_TEST_PATH/test_averagepool_1d_default",
        read_output("$ONNX_TEST_PATH/test_averagepool_1d_default"),
            read_input("$ONNX_TEST_PATH/test_averagepool_1d_default"))

# AveragePool 2D Default
ip = read_input("$ONNX_TEST_PATH/test_averagepool_2d_default")
main_test("$ONNX_TEST_PATH/test_averagepool_2d_default",
            read_output("$ONNX_TEST_PATH/test_averagepool_2d_default"),
                     read_input("$ONNX_TEST_PATH/test_averagepool_2d_default")[1])

# Averagepool 2d pads count include pad
ip = read_input("$ONNX_TEST_PATH/test_averagepool_2d_pads_count_include_pad")
main_test("$ONNX_TEST_PATH/test_averagepool_2d_pads_count_include_pad",
    read_output("$ONNX_TEST_PATH/test_averagepool_2d_pads_count_include_pad"),
        read_input("$ONNX_TEST_PATH/test_averagepool_2d_pads_count_include_pad")[1]);

# AveragePool 2d precomputed pads count include pad
ip = read_input("$ONNX_TEST_PATH/test_averagepool_2d_precomputed_pads_count_include_pad")
main_test("$ONNX_TEST_PATH/test_averagepool_2d_precomputed_pads_count_include_pad",
    read_output("$ONNX_TEST_PATH/test_averagepool_2d_precomputed_pads_count_include_pad"),
        read_input("$ONNX_TEST_PATH/test_averagepool_2d_precomputed_pads_count_include_pad")[1])

# Averagepool precomputed strides
ip = read_input("$ONNX_TEST_PATH/test_averagepool_2d_precomputed_strides")
    main_test("$ONNX_TEST_PATH/test_averagepool_2d_precomputed_strides",
        read_output("$ONNX_TEST_PATH/test_averagepool_2d_precomputed_strides"),
            read_input("$ONNX_TEST_PATH/test_averagepool_2d_precomputed_strides"))

# AveragePool 2D Strides
ip = read_input("$ONNX_TEST_PATH/test_averagepool_2d_strides")
main_test("$ONNX_TEST_PATH/test_averagepool_2d_strides",
            read_output("$ONNX_TEST_PATH/test_averagepool_2d_strides"),
                                     read_input("$ONNX_TEST_PATH/test_averagepool_2d_strides")[1])

# Test globalaveragepool
ip = read_input("$ONNX_TEST_PATH/test_globalaveragepool")
main_test("$ONNX_TEST_PATH/test_globalaveragepool",
    read_output("$ONNX_TEST_PATH/test_globalaveragepool"),
        read_input("$ONNX_TEST_PATH/test_globalaveragepool")[1])

# Test globalaveragepool precomputed
ip = read_input("$ONNX_TEST_PATH/test_globalaveragepool_precomputed")
main_test("$ONNX_TEST_PATH/test_globalaveragepool_precomputed",
    read_output("$ONNX_TEST_PATH/test_globalaveragepool_precomputed"),
        read_input("$ONNX_TEST_PATH/test_globalaveragepool_precomputed")[1])

# Test globalmaxpool
ip = read_input("$ONNX_TEST_PATH/test_globalmaxpool")
main_test("$ONNX_TEST_PATH/test_globalmaxpool",
    read_output("$ONNX_TEST_PATH/test_globalmaxpool"),
        read_input("$ONNX_TEST_PATH/test_globalmaxpool")[1])

# Test globalmaxpool precomputed
ip = read_input("$ONNX_TEST_PATH/test_globalmaxpool_precomputed")
main_test("$ONNX_TEST_PATH/test_globalmaxpool_precomputed",
    read_output("$ONNX_TEST_PATH/test_globalmaxpool_precomputed"),
        read_input("$ONNX_TEST_PATH/test_globalmaxpool_precomputed")[1])


text:
import Tables

Tables.isrowtable(::Type{<:StructArray}) = true

Tables.columnaccess(::Type{<:StructArray}) = true
Tables.columns(s::StructArray) = components(s)
Tables.schema(s::StructArray) = _schema(staticschema(eltype(s)))

_schema(::Type{NT}) where {NT<:NamedTuple} = Tables.Schema(NT)
# make schema for unnamed case
function _schema(::Type{T}) where {T<:NTuple{N, Any}} where N
    return Tables.Schema{ntuple(identity, N), T}
end

function try_compatible_columns(rows::R, s::StructArray) where {R}
    Tables.isrowtable(rows) && Tables.columnaccess(rows) || return nothing
    T = eltype(rows)
    hasfields(T) || return nothing
    NT = staticschema(T)
    _schema(NT) == Tables.schema(rows) || return nothing
    return Tables.columntable(rows)
end

for (f, g) in zip((:append!, :prepend!), (:push!, :pushfirst!))
    @eval function Base.$f(s::StructVector, rows)
        table = try_compatible_columns(rows, s)
        if table !== nothing
            # Input `rows` is a container of rows _and_ satisfies column
            # table interface.  Thus, we can add the input column-by-column.
            foreachfield($f, s, table)
            return s
        else
            # Otherwise, fallback to a generic implementation expecting
            # that `rows` is an iterator:
            return foldl($g, rows; init = s)
        end
    end
end


text:
using JSON
using JuMP, Ipopt
using Base

function read_json_file(fname)
    #TODO: write description
    dict2 = Dict()
    open(fname, "r") do f
        dicttxt = read(f, String)
        dict2 = JSON.parse(dicttxt)
    end
    return dict2
end

function read_nx_file(fname, flow)
    #TODO: write description
    G_dict = read_json_file(fname)
    G = Dict()
    nodes = []
    for edge in G_dict["links"]
        s = string(edge["source"])
        t = string(edge["target"])
        G[(s,t)] = Dict()
        push!(nodes, s)
        push!(nodes, t)
        nodes = unique(nodes)

        G[(s,t)]["source"] = s
        G[(s,t)]["target"] = t

        if flow == 0
            for att in ["capacity", "length", "t_0"]
                G[(s,t)][att] = edge[att]

            end
        else
            for att in ["capacity", "length", "t_0", "flow"]
                G[(s,t)][att] = edge[att]
            end
        end
    end
    return G, sort(nodes)

end

function read_fcoeffs(fname)
    #TODO: write description
    return read_json_file(fname)
end

function read_g(fname)
    #TODO: write description
    g_dict = Dict()
    dict2 = read_json_file(fname)
    for od in keys(dict2)
        a = split(split(split(od, "(")[2],")")[1],",")
        g_dict[(a[1], a[2][2:end])] = dict2[od]
    end
    return g_dict
end

function dict2json(dict, fname)
    json_string = JSON.json(dict)
    open(fname,"w") do f
        write(f, json_string)
    end
end

function get_incoming_edges(G, i)
    a = []
    a = unique([if G[k]["target"]==i G[k]["source"] else 0 end for k in keys(G)])
    filter!(x->x≠0,a)
    return a
end

function get_outgoing_edges(G, i)
    a = []
    a = unique([if G[k]["source"]==i G[k]["target"] else 0 end for k in keys(G)])
    filter!(x->x≠0,a)
    return a
end

function add_demand_cnstr(model, G, g, xc, nodes, bush)
    links = keys(G)
    od_pairs = keys(g)
    if bush==0
        for j in nodes
            for od in od_pairs
                if od[1] == j
                    @constraint(model, sum(xc[od][(i,j)] for i in get_incoming_edges(G, j)) + g[od] == sum(xc[od][(j,i)] for i in get_outgoing_edges(G, j)) )
                elseif od[2] == j
                    @constraint(model, sum(xc[od][(i,j)] for i in get_incoming_edges(G, j)) == g[od] + sum(xc[od][(j,i)] for i in get_outgoing_edges(G, j)) )
                else
                    @constraint(model, sum(xc[od][(i,j)] for i in get_incoming_edges(G, j)) == sum(xc[od][(j,i)] for i in get_outgoing_edges(G, j)) )
                end
            end
        end
    elseif bush==1
        O = unique([od[1] for od in od_pairs if g[od]>0])
        D = unique([od[2] for od in od_pairs if g[od]>0])
        od_nodes = unique([O;D])
        p = Dict()
        for j in od_nodes
            p[j] = sum(g[od] for od in od_pairs if od[2]==j) - sum(g[od] for od in od_pairs if od[1]==j)
        end
        # Global
        for j in od_nodes
            @constraint(model, sum(xc[s][(i,j)] for s in O for i in get_incoming_edges(G,j) ) - sum( xc[s][(j,k)] for s in O for k in get_outgoing_edges(G,j) ) == p[j])
        end
        # Local
        for s in O
            dsum = sum([g[od] for od in od_pairs if od[1]==s ])
            for j in nodes
                D = [od[2] for od in od_pairs if od[1]==s ]
                if j in D
                    d1 = g[(s,j)]
                    @constraint(model, sum([xc[s][(i,j)] for i in get_incoming_edges(G,j)]) - d1 == sum([xc[s][(j,k)] for k in get_outgoing_edges(G,j)]))
                elseif j == s
                    @constraint(model, sum([xc[s][(i,j)] for i in get_incoming_edges(G,j)]) == sum([xc[s][(j,k)] for k in get_outgoing_edges(G,j)])-dsum)
                else
                    @constraint(model, sum([xc[s][(i,j)] for i in get_incoming_edges(G,j)]) == sum([xc[s][(j,k)] for k in get_outgoing_edges(G,j)]))
                end
            end
        end
   end
end

#       for j in tnet.O:
#            p[j] = sum([tnet.g[(s,t)] for s,t in tnet.g.keys() if t==j]) - sum([tnet.g[(s,t)] for s,t in tnet.g.keys() if s==j])
        # Global
 #       for j in tnet.G_supergraph.nodes():
 #           m.addConstr(quicksum([x[(i,j)] for i,l in tnet.G_supergraph.in_edges(nbunch=j)]) - quicksum([x[(j,k)] for l,k in tnet.G_supergraph.out_edges(nbunch=j)]) == p[j] )
        # Local
  #      for s in tnet.O:
  #          dsum = sum([v for k,v in tnet.g.items() if k[0]==s])
  #          for j in tnet.G_supergraph.nodes():
  #              D = [d for o,d in tnet.g.keys() if o==s]
  #              if j in D:
  #                  d1 = tnet.g[(s,j)]
  #                  m.addConstr(quicksum(m.getVarByName('x^' + str(s) + '_' + str(i) + '_' + str(j)) for i, l in tnet.G_supergraph.in_edges(nbunch=j)) - d1 \
  #                      == quicksum(m.getVarByName('x^' + str(s) + '_' + str(j) + '_' + str(k)) for l, k in tnet.G_supergraph.out_edges(nbunch=j)))
  #              elif j == s:
  #                  m.addConstr(quicksum(m.getVarByName('x^' + str(s) + '_' + str(i) + '_' + str(j)) for i, l in tnet.G_supergraph.in_edges(nbunch=j)) \
  #                      == quicksum(m.getVarByName('x^' + str(s) + '_' + str(j) + '_' + str(k)) for l, k in tnet.G_supergraph.out_edges(nbunch=j)) - dsum)
  #              else:
  #                  m.addConstr(quicksum(m.getVarByName('x^' + str(s) + '_' + str(i) + '_' + str(j)) for i, l in tnet.G_supergraph.in_edges(nbunch=j)) \
  #                              == quicksum(m.getVarByName('x^' + str(s) + '_' + str(j) + '_' + str(k)) for l, k in tnet.G_supergraph.out_edges(nbunch=j)) )

function add_rebalancing_cnstr(m,G,xc,xr,nodes)
    od_pairs = keys(g)
    @constraint(m, [j in nodes], sum(xr[(i,j)] + sum(xc[od][(i,j)] for od in od_pairs) for i in get_incoming_edges(G, j)) == sum(xr[(j,k)] + sum(xc[od][(j,k)] for od in od_pairs) for k in get_outgoing_edges(G, j)))
end

function define_xc_vars(od_pairs, links, g, bush)
    xc = Dict()
    if bush==0
        for od in od_pairs
            xc[od] = Dict()
            for link in links
                xc[od][link] = @variable(m, lower_bound=0)
                set_name(xc[od][link], string("xc-(",od[1],",",od[2],")-(",link[1],",",link[2],")"))
            end
        end
    else
        O = unique([od[1] for od in od_pairs if g[od]>0])
        for o in O
            xc[o] = Dict()
            for link in links
                xc[o][link] = @variable(m, lower_bound=0)
                set_name(xc[o][link], string("xc-(",o,")-(",link[1],",",link[2],")"))
            end
        end
    end
    return xc
end

function define_xr_vars(links)
    xr = Dict()
    for link in links
        xr[link] = @variable(m, lower_bound=0)
        set_name(xr[link], string("xr-(",link[1],",",link[2],")"))
    end
    return xr
end

function sol2dict(m, G, xc, od_pairs, links, bush)
    dict = Dict()
    for link in links
        flow = 0
        if bush == 0
            for od in od_pairs
                var = string("xc-(",od[1],",",od[2],")-(",link[1],",",link[2],")")
                flow +=  value(variable_by_name(m, var))
            end
        else
            for o in O
                var = string("xc-(",o,")-(",link[1],",",link[2],")")
                flow +=  value(variable_by_name(m, var))
            end
        end
        #reb_var = value(variable_by_name(m, string("xr-(",link[1],",",link[2],")")))
        dict[link] = Dict()
        dict[link]["flow"] = flow
       # dict[link]["flowNoRebalancing"] = flow - reb_var
   end
   return dict
end

function export_results(m, G, xc, od_pairs, links, fname, bush)
    dict = sol2dict(m, G, xc, od_pairs, links, bush)
    dict2json(dict, fname)
end

function print_sol(m, xc, od_pairs, links, flag=1)
    if flag == 1
        for od in od_pairs
            xc[od] = Dict()
            for link in links
                var = string("xc-(",od[1],",",od[2],")-(",link[1],",",link[2],")")
                println(var, "= ", value(variable_by_name(m, var)))
            end
        end

        for link in links
            var =
            println(value(variable_by_name(m, string("xr-(",link[1],",",link[2],")"))))
        end
    else
        return 0
    end
end

#G, nodes = read_nx_file("tmp/G.json", 0)
G, nodes = read_nx_file("tmp/G_supergraph.json", 0)
G_exogenous, nodes_ex = read_nx_file("tmp/exogenous_G.json", 1)
#G_exogenous = "False"
fcoeffs = read_fcoeffs("tmp/fcoeffs.json")
g = read_g("tmp/g.json")
out_file = "tmp/out.json"

od_pairs = keys(g)
links = keys(G)
N = length(fcoeffs)

bush = 1

m = Model(with_optimizer(Ipopt.Optimizer))

# Define variables
xc = define_xc_vars(od_pairs, links, g, bush)
#xr = define_xr_vars(links)

# Define sum od flow
O = unique([od[1] for od in od_pairs if g[od]>0])
x = Dict()
for link in links
    if bush==0
        x[link] = @NLexpression(m, sum(xc[od][link] for od in od_pairs))
    else
        x[link] = @NLexpression(m, sum(xc[o][link] for o in O))
    end
end

## Define Objective
@NLobjective(m, Min, sum( x[link]* G[link]["t_0"]* sum( fcoeffs[n]*((x[link] + G_exogenous[link]["flow"])/G[link]["capacity"])^(n-1) for n=1:N ) for link in links))

## Define Constraints
add_demand_cnstr(m, G, g, xc, nodes, bush)

optimize!(m)

export_results(m, G, xc, od_pairs, links, out_file, bush)

# write runnning time to file
solvetime = MOI.get(m, MOI.SolveTime())
ofile= open("tmp/solvetime.txt", "w")
println(ofile, solvetime)
close(ofile)


text:
"""
    spgr! = SPGRBlochSim(TR, TE, rf, [spoiling], [nTR], [save_transients])
    spgr!(spin, [workspace])

Simulate a spoiled gradient-recalled echo (SPGR) scan on `spin`, overwriting the
spin's magnetization vector. The resultant magnetization is stored in `spin.M`.
If `nTR > 0` and `save_transients === true`, then `spgr!(...)` returns a
`Vector` with the magnetization vectors at the echo time for each of the `nTR`
simulated TRs.

# Arguments
- `TR::Real`: Repetition time (ms)
- `TE::Real`: Echo time (ms)
- `rf`:
  - `::Real`: Excitation flip angle (rad) (assumes an instantaneous RF pulse)
  - `::AbstractRF`: Excitation RF pulse
- `spoiling::AbstractSpoiling = IdealSpoiling()`: Type of spoiling to apply
- `nTR::Val = Val(0)`: Number of TRs to simulate; `Val(0)` indicates to simulate
  a steady-state scan
- `save_transients::Val = Val(false)`: Whether or not to return the
  magnetization vectors at the TE for each of the `nTR` simulated TRs; does
  nothing if `nTR == 0`

`workspace isa SPGRBlochSimWorkspace`.
"""
struct SPGRBlochSim{T1<:AbstractRF,T2<:AbstractSpoiling,nTR,save_transients}
    TR::Float64
    TE::Float64
    rf::T1
    spoiling::T2

    function SPGRBlochSim(
        TR,
        TE,
        rf::T1,
        spoiling::T2,
        nTR::Val{T3},
        save_transients::Val{T4}
    ) where {T1<:AbstractRF,T2<:AbstractSpoiling,T3,T4}

        Tg = spoiler_gradient_duration(spoiling)
        TR >= TE + Tg + duration(rf) / 2 ||
            error("TR must be greater than or equal to TE + Tg + duration(rf) / 2")
        TE >= duration(rf) / 2 || error("TE must not be during the excitation pulse")
        (T3 isa Int && T3 >= 0) || error("nTR must be a nonnegative Int")
        T2 <: Union{<:RFSpoiling,<:RFandGradientSpoiling} && (T3 > 0 ||
            error("nTR must be positive when simulating RF spoiling"))
        T4 isa Bool || error("save_transients must be a Bool")
        T3 == 0 && T4 &&
            @warn("save_transients is true, but nTR = 0; no transients will be saved")
        new{T1,T2,T3,T4}(TR, TE, rf, spoiling)

    end
end

SPGRBlochSim(TR, TE, rf, spoiling::AbstractSpoiling, nTR::Val = Val(0)) = SPGRBlochSim(TR, TE, rf, spoiling, nTR, Val(false))
SPGRBlochSim(TR, TE, rf, nTR::Val, save_transients::Val = Val(false)) = SPGRBlochSim(TR, TE, rf, IdealSpoiling(), nTR, save_transients)
SPGRBlochSim(TR, TE, rf) = SPGRBlochSim(TR, TE, rf, IdealSpoiling(), Val(0), Val(false))
SPGRBlochSim(TR, TE, α::Real, spoiling, nTR, save_transients) = SPGRBlochSim(TR, TE, InstantaneousRF(α), spoiling, nTR, save_transients)

Base.show(io::IO, spgr::SPGRBlochSim{<:AbstractRF,<:AbstractSpoiling,nTR,save}) where {nTR,save} =
    print(io, "SPGRBlochSim(", spgr.TR, ", ", spgr.TE, ", ", spgr.rf, ", ", spgr.spoiling, ", Val(", nTR, "), Val(", save, "))")

function Base.show(io::IO, ::MIME"text/plain", spgr::SPGRBlochSim{<:AbstractRF,<:AbstractSpoiling,nTR,save}) where {nTR,save}

    print(io, "Spoiled Gradient-Recalled Echo (SPGR) Bloch Simulation:")
    print(io, "\n TR = ", spgr.TR, " ms")
    print(io, "\n TE = ", spgr.TE, " ms")
    print(io, "\n rf (excitation pulse) = ")
    show(io, "text/plain", spgr.rf)
    print(io, "\n spoiling = ")
    show(io, "text/plain", spgr.spoiling)
    if nTR == 0
        print(io, "\n steady-state")
    else
        print(io, "\n nTR = ", nTR)
        print(io, "\n save transients = ", save)
    end

end

struct SPGRBlochSimWorkspace{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}
    Aex::T1
    Bex::T2
    Atr::T3
    Btr::T4
    Atg::T5
    Btg::T6
    tmpA1::T7
    tmpB1::T4
    tmpA2::T7
    tmpB2::T4
    mat::T8
    vec::T9
    bm_workspace::T10
    ex_workspace::T11
end

function SPGRBlochSimWorkspace(
    spin::AbstractSpin,
    scan::SPGRBlochSim,
    bm_workspace = spin isa Spin ? nothing : BlochMcConnellWorkspace(spin)
)

    SPGRBlochSimWorkspace(typeof(spin), typeof(scan), bm_workspace)

end

function SPGRBlochSimWorkspace(
    spin::Union{Type{Spin{T}},Type{SpinMC{T,N}}},
    scan::Type{<:SPGRBlochSim{T1,T2}},
    bm_workspace = spin isa Spin ? nothing : BlochMcConnellWorkspace(spin)
) where {T,N,T1,T2}

    if T1 <: InstantaneousRF
        Aex = ExcitationMatrix{T}()
        Bex = nothing
        ex_workspace = nothing
    else
        if spin <: Spin
            Aex = BlochMatrix{T}()
            Bex = Magnetization{T}()
        else
            Aex = BlochMcConnellMatrix{T}(N)
            Bex = MagnetizationMC{T}(N)
        end
        ex_workspace = ExcitationWorkspace(spin, bm_workspace)
    end
    if T2 <: IdealSpoiling
        Atg = idealspoiling
        Btg = nothing
    elseif T2 <: RFSpoiling
        Atg = nothing
        Btg = nothing
    elseif spin <: Spin
        Atg = FreePrecessionMatrix{T}()
        Btg = Magnetization{T}()
    else
        Atg = BlochMcConnellMatrix{T}(N)
        Btg = MagnetizationMC{T}(N)
    end
    if spin <: Spin
        Atr = FreePrecessionMatrix{T}()
        Btr = Magnetization{T}()
        tmpA1 = BlochMatrix{T}()
        tmpB1 = Magnetization{T}()
        tmpA2 = BlochMatrix{T}()
        tmpB2 = Magnetization{T}()
        mat = Matrix{T}(undef, 3, 3)
        vec = Vector{T}(undef, 3)
    else
        Atr = BlochMcConnellMatrix{T}(N)
        Btr = MagnetizationMC{T}(N)
        tmpA1 = BlochMcConnellMatrix{T}(N)
        tmpB1 = MagnetizationMC{T}(N)
        tmpA2 = BlochMcConnellMatrix{T}(N)
        tmpB2 = MagnetizationMC{T}(N)
        mat = Matrix{T}(undef, 3N, 3N)
        vec = Vector{T}(undef, 3N)
    end
    SPGRBlochSimWorkspace(Aex, Bex, Atr, Btr, Atg, Btg, tmpA1, tmpB1, tmpA2,
                          tmpB2, mat, vec, bm_workspace, ex_workspace)

end

# Case when nTR = 0
function (scan::SPGRBlochSim{<:AbstractRF,<:AbstractSpoiling,0})(
    spin::AbstractSpin,
    workspace::SPGRBlochSimWorkspace = SPGRBlochSimWorkspace(spin, scan)
)

    Tg = spoiler_gradient_duration(scan.spoiling)

    excite!(workspace.Aex, workspace.Bex, spin, scan.rf, workspace.ex_workspace)
    freeprecess!(workspace.Atr, workspace.Btr, spin, scan.TR - Tg - duration(scan.rf), workspace.bm_workspace)
    spoil!(workspace.Atg, workspace.Btg, spin, scan.spoiling, workspace.bm_workspace)

    combine!(workspace.tmpA1, workspace.tmpB1, workspace.Atr, workspace.Btr, workspace.Atg, workspace.Btg)
    combine!(workspace.tmpA2, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1, workspace.Aex, workspace.Bex)
    subtract!(workspace.mat, I, workspace.tmpA2)
    copyto!(workspace.vec, workspace.tmpB2)
    F = lu!(workspace.mat)
    ldiv!(F, workspace.vec)
    copyto!(spin.M, workspace.vec)

    freeprecess!(workspace.Atr, workspace.Btr, spin, scan.TE - duration(scan.rf) / 2, workspace.bm_workspace)
    applydynamics!(spin, workspace.tmpB1, workspace.Atr, workspace.Btr)

end

# Case when nTR > 0
function (scan::SPGRBlochSim{<:AbstractRF,T,nTR,save})(
    spin::AbstractSpin,
    workspace::SPGRBlochSimWorkspace = SPGRBlochSimWorkspace(spin, scan)
) where {T,nTR,save}

    rf = scan.rf
    rfspoiling = T <: Union{<:RFSpoiling,<:RFandGradientSpoiling}
    Tg = spoiler_gradient_duration(scan.spoiling)

    if rfspoiling
        rf isa RF && (rf.Δθ[] = rf.Δθ_initial)
        Δθinc = rfspoiling_increment(scan.spoiling)
        θ = zero(Δθinc) # For knowing how much phase to remove when recording signal
        Δθ = Δθinc
    else
        excite!(workspace.Aex, workspace.Bex, spin, rf, workspace.ex_workspace)
    end

    if save
        M = Vector{typeof(spin.M)}(undef, nTR)
        freeprecess!(workspace.Atr, workspace.Btr, spin, scan.TR - scan.TE - Tg - duration(scan.rf) / 2, workspace.bm_workspace)
    else
        freeprecess!(workspace.Atr, workspace.Btr, spin, scan.TR - Tg - duration(scan.rf), workspace.bm_workspace)
    end
    spoil!(workspace.Atg, workspace.Btg, spin, scan.spoiling, workspace.bm_workspace)
    combine!(workspace.tmpA1, workspace.tmpB1, workspace.Atr, workspace.Btr, workspace.Atg, workspace.Btg)
    freeprecess!(workspace.Atr, workspace.Btr, spin, scan.TE - duration(scan.rf) / 2, workspace.bm_workspace)

    for rep = 1:nTR-1

        rfspoiling && excite!(workspace.Aex, workspace.Bex, spin, rf, workspace.ex_workspace)
        applydynamics!(spin, workspace.tmpB2, workspace.Aex, workspace.Bex)
        if save
            applydynamics!(spin, workspace.tmpB2, workspace.Atr, workspace.Btr)
            M[rep] = copy(spin.M)
            if rfspoiling
                modulation = exp(im * θ)
                if spin isa Spin
                    tmp = signal(spin.M) * modulation
                    M[rep].x = real(tmp)
                    M[rep].y = imag(tmp)
                else
                    for i = 1:spin.N
                        tmp = signal(spin.M[i]) * modulation
                        M[rep][i].x = real(tmp)
                        M[rep][i].y = imag(tmp)
                    end
                end
            end
        end
        applydynamics!(spin, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1)

        if rfspoiling
            if rf isa InstantaneousRF
                rf = InstantaneousRF(rf.α, rf.θ + Δθ)
            else
                rf.Δθ[] += Δθ
            end
            θ += Δθ
            Δθ += Δθinc
        end

    end

    rfspoiling && excite!(workspace.Aex, workspace.Bex, spin, rf, workspace.ex_workspace)
    applydynamics!(spin, workspace.tmpB2, workspace.Aex, workspace.Bex)
    applydynamics!(spin, workspace.tmpB2, workspace.Atr, workspace.Btr)
    if rfspoiling
        modulation = exp(im * θ)
        if spin isa Spin
            tmp = signal(spin.M) * modulation
            spin.M.x = real(tmp)
            spin.M.y = imag(tmp)
        else
            for i = 1:spin.N
                tmp = signal(spin.M[i]) * modulation
                spin.M[i].x = real(tmp)
                spin.M[i].y = imag(tmp)
            end
        end
    end
    if save
        M[nTR] = copy(spin.M)
        return M
    end

end


text:
using CthulhuVision.Image
using CthulhuVision.Math
using CthulhuVision.Light
using CthulhuVision.Materials
using CthulhuVision.Rendering
using CthulhuVision.Camera
using CthulhuVision.Worlds
using CthulhuVision.Scenes
using CthulhuVision.Random

################
# Define image #
################

width = 512
height = 512
image = PPM(Dimension(width, height))

####################
# Define materials #
####################

red = lambertian(RGB(0.9f0, 0.1f0, 0.1f0))
green = lambertian(RGB(0.1f0, 0.9f0, 0.1f0))
blue = lambertian(RGB(0.1f0, 0.1f0, 0.9f0))
white = lambertian(RGB(1.0f0, 1.0f0, 1.0f0))
light = lambertian(RGB(1f0, 1f0, 1f0); emission = RGB(10f0, 10f0, 10f0))
glass = dielectric(1.5f0)

#################
# SceneSettings #
#################

ambientemission = RGB(0.0f0, 0.0f0, 0.0f0)
settings = SceneSettings(ambientemission)

##########
# Camera #
##########

aspect = Float32(image.dimension.width / image.dimension.height)
vfov = 40.0f0

lookfrom = Vector3(0f0, 0f0, 800.0f0)
lookat = Vector3(0f0, 0f0, 0.0f0)
vup = Vector3(0.0f0, 1.0f0, 0.0f0)
focusdist = 10.0f0
aperture = 0.0f0

camera = FovCamera(lookfrom, lookat, vup, vfov, aspect, aperture, focusdist)

###################
# Construct scene #
###################

function constructscene() :: Scene
    scene = Scene(settings)

    wallmesh = Mesh(
        # Vertexes
        [
            Vector3(-277.5f0, -277.5f0, 0f0),
            Vector3( 277.5f0, -277.5f0, 0f0),
            Vector3( 277.5f0,  277.5f0, 0f0),
            Vector3(-277.5f0,  277.5f0, 0f0),
        ],

        # Triangles
        [
            MeshTriangle(1, 2, 4),
            MeshTriangle(2, 3, 4),
        ]
    )
    wallmeshindex = addmesh!(scene, wallmesh)

    shortmesh = Mesh(
        # Vertexes
        [
            # Front vertexes
            Vector3(-82.5f0, -82.5f0, 82.5f0),
            Vector3( 82.5f0, -82.5f0, 82.5f0),
            Vector3( 82.5f0,  82.5f0, 82.5f0),
            Vector3(-82.5f0,  82.5f0, 82.5f0),

            # Back vertexes
            Vector3(-82.5f0, -82.5f0, -82.5f0),
            Vector3( 82.5f0, -82.5f0, -82.5f0),
            Vector3( 82.5f0,  82.5f0, -82.5f0),
            Vector3(-82.5f0,  82.5f0, -82.5f0),
        ],

        # Triangles
        [
            # Front face
            MeshTriangle(1, 2, 4),
            MeshTriangle(2, 3, 4),

            # Back face
            MeshTriangle(6, 5, 7),
            MeshTriangle(5, 8, 7),

            # Left face
            MeshTriangle(5, 1, 8),
            MeshTriangle(1, 4, 8),

            # Right face
            MeshTriangle(2, 6, 3),
            MeshTriangle(6, 7, 3),

            # Top face
            MeshTriangle(4, 3, 8),
            MeshTriangle(3, 7, 8),

            # Bottom face
            MeshTriangle(5, 6, 1),
            MeshTriangle(6, 2, 1),
        ]
    )
    shortmeshindex = addmesh!(scene, shortmesh)

    lightmesh = Mesh(
        # Vertexes
        [
            Vector3(-65f0, -52.5f0, 0f0),
            Vector3( 65f0, -52.5f0, 0f0),
            Vector3( 65f0,  52.5f0, 0f0),
            Vector3(-65f0,  52.5f0, 0f0),
        ],

        # Triangles
        [
            MeshTriangle(1, 2, 4),
            MeshTriangle(2, 3, 4),
        ]
    )
    lightmeshindex = addmesh!(scene, lightmesh)

    # Back wall
    addinstance!(
        scene,
        MeshInstance(
            wallmeshindex,
            white,
            translation(Vector3(0f0, 0f0, -555f0)),
        )
    )

    # Bottom
    addinstance!(
        scene,
        MeshInstance(
            wallmeshindex,
            white,
            translation(Vector3(0f0, -277.5f0, -277.5f0)) * rotation(-Float32(π) / 2f0, Vector3(1f0, 0f0, 0f0)),
        )
    )

    # Top
    addinstance!(
        scene,
        MeshInstance(
            wallmeshindex,
            white,
            translation(Vector3(0f0,  277.5f0, -277.5f0)) * rotation(Float32(π) / 2f0, Vector3(1f0, 0f0, 0f0)),
        )
    )

    # Left wall
    addinstance!(
        scene,
        MeshInstance(
            wallmeshindex,
            red,
            translation(Vector3(-277.5f0, 0f0, -277.5f0)) * rotation(Float32(π) / 2f0, Vector3(0f0, 1f0, 0f0))
        )
    )

    # Right wall
    addinstance!(
        scene,
        MeshInstance(
            wallmeshindex,
            green,
            translation(Vector3( 277.5f0, 0f0, -277.5f0)) * rotation(-Float32(π) / 2f0, Vector3(0f0, 1f0, 0f0))
        )
    )

    # Short block
    addinstance!(
        scene,
        MeshInstance(
            shortmeshindex,
            white,
            # identitytransform(),
            translation(Vector3(87.5f0, -195f0, -147.5f0)) * rotation(-Float32(2.0 * π) * 0.05f0, Vector3(0f0, 1f0, 0f0)),
        )
    )

    # Glass block
    addinstance!(
        scene,
        MeshInstance(
            shortmeshindex,
            glass,
            # identitytransform(),
            translation(Vector3(-97.5f0, -195f0, -147.5f0)) * rotation(-Float32(2.0 * π) * 0.05f0, Vector3(0f0, 1f0, 0f0)) * scale(Vector3(0.7f0, 0.7f0, 0.7f0)),
        )
    )

    # Tall block
    addinstance!(
        scene,
        MeshInstance(
            shortmeshindex,
            white,
            # identitytransform(),
            translation(Vector3(-97.5f0, -112.5f0, -347.5f0)) * rotation(Float32(2.0 * π) * 0.041667f0, Vector3(0f0, 1f0, 0f0)) * scale(Vector3(1f0, 2f0, 1f0)),
        )
    )

    # Light
    addinstance!(
        scene,
        MeshInstance(
            lightmeshindex,
            light,
            translation(Vector3(0f0,  277.4f0, -277.5f0)) * rotation(Float32(π) / 2f0, Vector3(1f0, 0f0, 0f0)),
        )
    )

    scene
end

##################
# Perform render #
##################

scene = constructscene()
rendersettings = RenderSettings(30000)

render(image, camera, scene, rendersettings)

saveimage(image, examplefilename("cornellbox"))