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

export NfOrdQuoRing, NfOrdQuoRingElem, quo, *, -, ==, deepcopy, divrem,
       gcd, inv, parent, show, divexact, isone, iszero, howell_form,
       strong_echelon_form, triangularize, det, euclid, xxgcd, annihilator

################################################################################
#
#  Assert
#
################################################################################

add_assert_scope(:NfOrdQuoRing)

set_assert_level(:NfOrdQuoRing, 0)

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

function elem_type(::Type{AbsOrdQuoRing{S, T}}) where {S, T}
  U = elem_type(S)
  return AbsOrdQuoRingElem{S, T, U}
end

function elem_type(::AbsOrdQuoRing{S, T}) where {S, T}
  U = elem_type(S)
  return AbsOrdQuoRingElem{S, T, U}
end

base_ring(Q::AbsOrdQuoRing) = Q.base_ring

ideal(Q::AbsOrdQuoRing) = Q.ideal

basis_matrix(Q::AbsOrdQuoRing) = Q.basis_matrix

parent(x::AbsOrdQuoRingElem) = x.parent

parent_type(::Type{AbsOrdQuoRingElem{S, T, U}}) where {S, T, U} = AbsOrdQuoRing{S, T}

(R::AbsOrdQuoRing)() = zero(R)

################################################################################
#
#  Hashing
#
################################################################################

hash(x::AbsOrdQuoRingElem, h::UInt) = hash(mod(x.elem, parent(x)), h)

################################################################################
#
#  Functions to allow polynomial and polynomial ring constructions
#
################################################################################

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

Nemo.promote_rule(::Type{NfOrdQuoRingElem}, ::Type{fmpz}) = NfOrdQuoRingElem

################################################################################
#
#  Copying
#
################################################################################

Base.deepcopy_internal(x::AbsOrdQuoRingElem, dict::IdDict) =
        AbsOrdQuoRingElem(parent(x), Base.deepcopy_internal(x.elem, dict))

#copy(x::NfOrdQuoRingElem) = deepcopy(x)

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

function show(io::IO, Q::AbsOrdQuoRing)
  print(io, "Quotient of $(Q.base_ring)")
end

function AbstractAlgebra.expressify(x::AbsOrdQuoRingElem; context = nothing)
  return AbstractAlgebra.expressify(x.elem, context = context)
end

function show(io::IO, x::AbsOrdQuoRingElem)
  print(io, AbstractAlgebra.obj_to_string(x, context = io))
end

################################################################################
#
#  Easy reduction of elements
#
################################################################################

#TODO: Inplace versions of mod
function _easy_mod(x::NfOrdQuoRingElem)
  Q = parent(x)
  I = Q.ideal
  O = parent(x.elem)
  if isdefining_polynomial_nice(nf(O)) && contains_equation_order(O)
    x.elem = O(mod(x.elem.elem_in_nf, minimum(I, copy = false)), false)
  else
    x.elem = mod(x.elem, I)
  end
  return x
end

function _easy_mod(x::AbsOrdQuoRingElem)
  x.elem = mod(x.elem, parent(x).ideal)
  return x
end

################################################################################
#
#  Parent object overloading
#
################################################################################

function (Q::AbsOrdQuoRing{S, T})(x::U) where {S, T, U}
  parent(x) !== base_ring(Q) && error("Cannot coerce element into the quotient ring")
  res = AbsOrdQuoRingElem{S, T, U}(Q, x)
  return _easy_mod(res)
end

function (Q::AbsOrdQuoRing{S, T})(x::Integer) where {S, T}
  res = elem_type(Q)(Q, base_ring(Q)(x))
  return res
end

function (Q::AbsOrdQuoRing{S, T})(x::fmpz) where {S, T}
  res = elem_type(Q)(Q, base_ring(Q)(x))
  return res
end

################################################################################
#
#  Quotient function
#
# (and standard helpers)
#
################################################################################

@doc Markdown.doc"""
    quo(O::NfOrd, I::NfOrdIdl) -> NfOrdQuoRing, Map
    quo(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl) -> AbsOrdQuoRing, Map

The quotient ring $O/I$ as a ring together with the section $M: O/I \to O$.
The pointwise inverse of $M$ is the canonical projection $O\to O/I$.
"""
function quo(O::Union{NfAbsOrd, AlgAssAbsOrd}, I::Union{NfAbsOrdIdl, AlgAssAbsOrdIdl})
  @assert order(I) === O
  # We should check that I is not zero
  Q = AbsOrdQuoRing(O, I)
  f = AbsOrdQuoMap(O, Q)
  return Q, f
end

function quo(O::Union{NfAbsOrd, AlgAssAbsOrd}, I::Union{NfAbsOrdIdl, AlgAssAbsOrdIdl}, ::Type{GrpAbFinGen})
  f = GrpAbFinGenToNfOrdQuoNfOrd(O, I)
  return domain(f), f
end

@doc Markdown.doc"""
    ResidueRing(O::NfOrd, I::NfOrdIdl) -> NfOrdQuoRing
    ResidueRing(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl) -> AbsOrdQuoRing

The quotient ring $O$ modulo $I$ as a new ring.
"""
Nemo.ResidueRing(O::Union{NfAbsOrd, AlgAssAbsOrd}, I::Union{NfAbsOrdIdl, AlgAssAbsOrdIdl}) = AbsOrdQuoRing(O, I)

@doc Markdown.doc"""
    lift(O::NfOrd, a::NfOrdQuoRingElem) -> NfOrdElem

Returns a lift of $a$ back to $O$.
"""
function lift(O::NfOrd, a::NfOrdQuoRingElem)
  f = NfOrdQuoMap(O, parent(a))
  return preimage(f, a)
end

@doc Markdown.doc"""
    lift(a::NfOrdQuoRingElem) -> NfOrdElem

Given an element of the quotient ring $\mathcal O/I$, return a lift in
$\mathcal O$.
"""
function lift(a::NfOrdQuoRingElem)
  mod!(a.elem, parent(a))
  return a.elem
end

################################################################################
#
#  Parent check
#
################################################################################

function check_parent(x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  if parent(x) !== parent(y)
    error("Elements must have same parents")
  end
  return true
end

################################################################################
#
#  Arithmetic
#
################################################################################

function +(x::AbsOrdQuoRingElem{S, T, U}, y::AbsOrdQuoRingElem{S, T, U}) where {S, T, U}
  check_parent(x, y)
  Q = parent(x)
  return Q(x.elem + y.elem)
end

function -(x::AbsOrdQuoRingElem{S, T, U}, y::AbsOrdQuoRingElem{S, T, U}) where {S, T, U}
  check_parent(x, y)
  Q = parent(x)
  return Q(x.elem - y.elem)
end

function -(x::AbsOrdQuoRingElem{S, T, U}) where {S, T, U}
  Q = parent(x)
  return AbsOrdQuoRingElem{S, T, U}(Q, -x.elem)
end

function *(x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  check_parent(x, y)
  Q = parent(x)
  return Q(x.elem * y.elem)
end

function mul!(z::AbsOrdQuoRingElem, x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  z.elem = mul!(z.elem, x.elem, y.elem)
  return _easy_mod(z)
end

function add!(z::AbsOrdQuoRingElem, x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  z.elem = add!(z.elem, x.elem, y.elem)
  return _easy_mod(z)
end

addeq!(x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem) = add!(x, x, y)

function sub!(z::AbsOrdQuoRingElem, x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  z.elem = sub!(z.elem, x.elem, y.elem)
  return _easy_mod(z)
end

function *(x::T, y::AbsOrdQuoRingElem) where T <: Union{Integer, fmpz}
  Q = parent(y)
  return Q(x*y.elem)
end

*(x::AbsOrdQuoRingElem, y::T) where T <: Union{Integer, fmpz} = y*x

function ^(a::AbsOrdQuoRingElem, f::fmpz)
  if fits(Int, f)
    return a^Int(f)
  end
  f == 0 && return one(parent(a))
  f == 1 && return deepcopy(a)
  if f < 0
    f = -f
    a = inv(a)
  end
  b = a^(div(f, 2))
  b *= b
  if isodd(f)
    b *= a
  end
  return b
end

#^(a::NfOrdQuoRingElem, f::Integer) = a^fmpz(f)
function ^(a::NfOrdQuoRingElem, b::Int)
  if b < 0
    a = inv(a)
    b = -b
  end
  Q = parent(a)
  O = base_ring(Q)
  if !isdefining_polynomial_nice(nf(O)) || !contains_equation_order(O)
    return pow1(a, b)
  end
  m = minimum(Q.ideal, copy = false)
  x = _powermod(a.elem.elem_in_nf, b, m)
  return Q(O(x))
end

function ^(a::AbsOrdQuoRingElem, b::Int)
  return pow1(a, b)
end

function pow1(a::AbsOrdQuoRingElem, b::Int)
  if b == 0
    return one(parent(a))
  elseif b == 1
    return deepcopy(a)
  else
    if b < 0
      a = inv(a)
      b = -b
    end
    bit = ~((~UInt(0)) >> 1)
    while (UInt(bit) & b) == 0
      bit >>= 1
    end
    z = deepcopy(a)
    bit >>= 1
    while bit != 0
      z = mul!(z, z, z)
      if (UInt(bit) & b) != 0
        z = mul!(z, z, a)
      end
      bit >>= 1
    end
    return z
  end
end

################################################################################
#
#  Special elements
#
################################################################################

function iszero(x::AbsOrdQuoRingElem)
  if iszero(x.elem)
    return true
  end
  mod!(x.elem, parent(x))
  return iszero(x.elem)
end

function isone(x::AbsOrdQuoRingElem)
  if isone(x.elem)
    return true
  end
  mod!(x.elem, parent(x))
  return isone(x.elem)
end

function one(Q::AbsOrdQuoRing)
  return elem_type(Q)(Q, base_ring(Q)(1))
end

function zero(Q::AbsOrdQuoRing)
  return elem_type(Q)(Q, base_ring(Q)(0))
end

function zero!(x::AbsOrdQuoRingElem)
  zero!(x.elem)
  return x
end

################################################################################
#
#  Equality
#
################################################################################

function ==(x::AbsOrdQuoRing, y::AbsOrdQuoRing)
  return base_ring(x) === base_ring(y) && ideal(x) == ideal(y)
end

function ==(x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  mod!(x.elem, parent(x))
  mod!(y.elem, parent(x))
  return x.elem == y.elem
end

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

function divexact(x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  b, z = isdivisible(x, y)
  @assert b
  return z
end

function isdivisible2(x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  check_parent(x, y)

  iszero(y) && error("Dividing by zero")

  if iszero(x)
    return true, zero(parent(x))
  end

  R = parent(x)
  O = base_ring(R)
  d = degree(O)

  A = representation_matrix_mod(y.elem, minimum(R.ideal))
  B = parent(x).basis_matrix
  V = hcat(A', B')

  a = coordinates(x.elem, copy = false)
  rhs = matrix(FlintZZ, d, 1, a)

  fl, sol = cansolve(V, rhs)
  if !fl
    return fl, zero(R)
  end
  z = R(O(fmpz[sol[i, 1] for i = 1:degree(O)]))
  @hassert :NfOrdQuoRing 1 z*y == x
  return true, z
end

function isdivisible(x::AbsOrdQuoRingElem, y::AbsOrdQuoRingElem)
  check_parent(x, y)

  iszero(y) && error("Dividing by zero")

  if iszero(x)
    return true, zero(parent(x))
  end

  R = parent(x)
  d = degree(base_ring(R))
  # We cannot solve with non-square matrices.
  # Thus build the matrix
  # ( 1   x    0  )
  # ( 0  M_y   I  )
  # ( 0  M_I   0  ).
  # Compute the UPPER RIGHT HNF ->
  # ( 1   0   u )
  # ( *   *   * )
  # u will be the coefficient vector of the quotient

  V = R.tmp_div
  A = representation_matrix_mod(y.elem, minimum(R.ideal))
  B = parent(x).basis_matrix

  V[1, 1] = 1

  a = coordinates(x.elem, copy = false)

  for i in 1:d
    V[1, 1 + i] = a[i]
  end

  _copy_matrix_into_matrix(V, 2, 2, A)   # this really is a copy
  _copy_matrix_into_matrix(V, 2 + d, 2, B) # this really is a copy

  for i in 1:d
    V[1 + i, d + 1 + i] = 1
  end

  if typeof(base_ring(parent(x))) <: NfOrd
    hnf_modular_eldiv!(V, minimum(R.ideal))
  else
    hnf!(V)
  end

  for i in 2:(d + 1)
    if !iszero(V[1, i])
  #if !iszero(sub(V, 1:1, 2:(d + 1)))
      ccall((:fmpz_mat_zero, libflint), Nothing, (Ref{fmpz_mat}, ), V)
      return false, zero(parent(x))
    end
  end

  z = R(-base_ring(R)(fmpz[ V[1, i] for i in (d + 2):(2*d + 1)])) # V[1, i] is always a copy

  ccall((:fmpz_mat_zero, libflint), Nothing, (Ref{fmpz_mat}, ), V)

  @hassert :NfOrdQuoRing 1 z*y == x
  return true, z
end

################################################################################
#
#  Strong exact division
#
################################################################################

function _divexact_strong(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  n = euclid(x)
  m = euclid(y)
  @hassert :NfOrdQuoRing 1 mod(n, m) == 0
  target = divexact(n, m)

  #println("target valuation: $target")
  #println("doing first an ordinary divexact with\n $x \n and \n $y")
  q0 = divexact(x, y)
  #println("valuation of first quotient: $(euclid(q0))")

  if euclid(q0) == target
    return q0
  else
    i = annihilator(y)
    #println("generator of annihilator: $i")

    q = q0 + rand(parent(x))*i

    k = 0
    while euclid(q) != target
      k += 1
      q = q0 + rand(parent(x))*i
      #println("current valuation $(euclid(q))")
      if k > 500
        error("Could not find proper quotion for strong division")
      end
    end
  end

  @hassert :NfOrdQuoRing 1 q*y == x
  @hassert :NfOrdQuoRing 1 euclid(q) *euclid(y) == euclid(x)

  return q
end

################################################################################
#
#  Inverse element
#
################################################################################

function isinvertible(x::AbsOrdQuoRingElem)
  if iszero(x)
    return false, x
  end
  return isdivisible(one(parent(x)), x)
end

function inv(x::AbsOrdQuoRingElem)
  t, y = isinvertible(x)
  @assert t "Element is not invertible"
  return y
end

################################################################################
#
#  Unit
#
################################################################################

isunit(x::AbsOrdQuoRingElem) = isinvertible(x)[1]

################################################################################
#
#  Euclidean function
#
################################################################################

function euclid(x::NfOrdQuoRingElem)
  if iszero(x)
    return fmpz(-1)
  end

  U = parent(x).tmp_euc

  d = degree(base_ring(parent(x)))

  _copy_matrix_into_matrix(U, 1, 1, representation_matrix(x.elem))
  _copy_matrix_into_matrix(U, d + 1, 1, parent(x).basis_matrix)

  hnf_modular_eldiv!(U, minimum(parent(x).ideal, copy = false))

  z = fmpz(1)

  for i in 1:degree(base_ring(parent(x)))
    mul!(z, z, U[i, i])
  end

  @hassert :NfOrdQuoRing 1 z == norm(ideal(parent(x.elem), x.elem) + parent(x).ideal)

  return z
end

################################################################################
#
#  Division with remainder
#
################################################################################

function Base.divrem(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  b, q = isdivisible(x, y)
  if b
    return q, zero(parent(x))
  end

  e = euclid(y)

  cnt = 0
  while true
    cnt += 1
    q = rand(parent(x))
    r = x - q*y
    if euclid(r) < e
      @hassert :NfOrdQuoRing 1 euclid(r) < e
      return q, r
    end
    if cnt > 1000
      error("Something odd in divrem for $x $y $(parent(x))")
    end
  end
end

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

Random.gentype(::Type{NfOrdQuoRing}) = elem_type(NfOrdQuoRing)

function rand(rng::AbstractRNG, Qsp::Random.SamplerTrivial{NfOrdQuoRing})
  Q = Qsp[]
  A = basis_matrix(Q)
  B = basis(base_ring(Q))
  z = rand(rng, fmpz(1):A[1,1]) * B[1]

  for i in 2:nrows(A)
    z = z + rand(rng, fmpz(1):A[i, i]) * B[i]
  end

  return Q(z)
end

################################################################################
#
#  Annihilator
#
################################################################################

function annihilator(x::NfOrdQuoRingElem)
  I = parent(x).ideal
  O = base_ring(parent(x))
  d = degree(O)
  f = NfOrdQuoMap(O, parent(x))
  U = parent(x).tmp_ann
  # We first compute a Z-basis of the annihilator ideal
  # The basis is the kernel of the following matrix
  # ( M_I )
  # ( M_x )
   _copy_matrix_into_matrix(U, 1, 1, representation_matrix(x.elem))
   _copy_matrix_into_matrix(U, d + 1, 1, I.basis_matrix)

  m = left_kernel(U)[2]
  I = ideal(O, _hnf_modular_eldiv(sub(m, 1:degree(O), 1:degree(O)),
                                  minimum(I), :lowerleft))
  z = f(I)
  #@assert iszero(z*x)
  return z
end

################################################################################
#
#  GCD
#
################################################################################

function gcd(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  Q = parent(x)
  O = base_ring(Q)

  if iszero(x)
    if iszero(y)
      return zero(Q)
    else
      return y
    end
  elseif iszero(y)
    return x
  end

  if isone(x) || isone(y)
    return one(Q)
  end

  I = ideal(O, x.elem) + ideal(O, y.elem)

  f = NfOrdQuoMap(O, parent(x))

  return f(I)
end

################################################################################
#
#  Extended extended GCD
#
################################################################################

function xxgcd(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  Q = parent(x)
  O = base_ring(Q)

  d = degree(O)

  if iszero(x)
    return deepcopy(y), Q(0), Q(1), Q(-1), Q(0)
  elseif iszero(y)
    return deepcopy(x), Q(1), Q(0), Q(0), Q(1)
  end

  g = gcd(x, y)

  e = _divexact_strong(x, g)
  f = _divexact_strong(y, g)

  M_e = representation_matrix(e.elem)
  M_f = representation_matrix(f.elem)

  M_I = Q.basis_matrix

  # let us build
  # ( 1  Q(1) 0  0 )
  # ( 0  M_e  I  0 )
  # ( 0  M_f  0  I )
  # ( 0  M_I  0  0 )

  a = coordinates(one(O), copy = false)

  V = parent(x).tmp_xxgcd

  V[1, 1] = 1

  for i in 1:d
    V[1, 1 + i] = a[i]
  end

  Hecke._copy_matrix_into_matrix(V, 2, 2, M_e)
  Hecke._copy_matrix_into_matrix(V, d + 2, 2, M_f)
  Hecke._copy_matrix_into_matrix(V, 2*d + 2, 2, M_I)

  for i in 1:2*d
    V[1+i, 1 + d + i] = 1
  end

  #U = V

  hnf_modular_eldiv!(V, minimum(Q.ideal, copy = false))

  u = Q(-O([ V[1,i] for i in (d + 2):(2*d + 1)]))
  v = Q(-O([ V[1,i] for i in (2*d + 2):(3*d + 1)]))

  @hassert :NfOrdQuoRing 1 Q(O(1)) == u*e - (v*(-f))

  ccall((:fmpz_mat_zero, libflint), Nothing, (Ref{fmpz_mat}, ), V)

  return g, u, v, -f, e
end

function (M::Generic.MatSpace{NfOrdQuoRingElem})(x::Generic.MatSpaceElem{NfOrdElem})
  z = map(base_ring(M), x.entries)::Array{NfOrdQuoRingElem, 2}
  return M(z)::Generic.MatSpaceElem{NfOrdQuoRingElem}
end

################################################################################
#
#  Compute the probability
#
################################################################################

# Just for debugging purposes
# This governs the runtime of the probablistic algorithms.
function probablity(O::NfOrdQuoRing)
  p = 1.0
  I = O.ideal
  f = factor(norm(I))
  for (P, e) in f
    if valuation(I, P) > 0
      p = p * (1 - 1/Float64(norm(P)))
    end
  end
  return p
end

################################################################################
#
#  Group Structure
#
################################################################################

@doc Markdown.doc"""
    group_structure(Q::NfOrdQuoRing) -> GrpAbFinGen

Returns an abelian group with the structure of $(Q,+)$.
"""
function group_structure(Q::NfOrdQuoRing)
  i = ideal(Q)
  fac = factor(i)
  structure = Vector{fmpz}()
  for (p,vp) in fac
    pnum = minimum(p)
    e = valuation(pnum,p)
    f = factor(norm(p))[pnum]
    q, r = divrem(vp+e-1,e)
    structure_pvp = [repeat([pnum^q],inner=[Int((r+1)*f)]) ; repeat([pnum^(q-1)],inner=[Int((e-r-1)*f)])]
    append!(structure,structure_pvp)
  end
  G = abelian_group(structure)
  S, Smap = snf(G)
  return S
end


text:
#=
Macro: A macroeconomic model for LEAP

Authors:
  Eric Kemp-Benedict
  Jason Veysey
  Emily Ghosh
  Anisha Nazareth

Copyright ©2019-2022 Stockholm Environment Institute
=#

module LEAPMacro

using Logging, Dates, YAML

include("./IOlib.jl")
include("./LEAPfunctions.jl")
include("./MacroModel.jl")
using .IOlib, .LEAPfunctions, .MacroModel

function run(config_file::AbstractString = "LEAPMacro_params.yml"; dump_err_stack::Bool = false, include_energy_sectors::Bool = false)
	# Ensure needed folders exist
	if !isdir("inputs")
		throw(ErrorException("The \"inputs\" folder is needed, but does not exist"))
	end

	params = YAML.load_file(config_file)

	# Enable logging
	if include_energy_sectors
		logfilename = string("LEAPMacro_log_", params["output_folder"], "_full.txt")
	else
		logfilename = string("LEAPMacro_log_", params["output_folder"], ".txt")
	end
	logfile = open(logfilename, "w+")
	logger = ConsoleLogger(logfile)
	global_logger(logger)
	@info now()
	@info "Configuration file: '$config_file'"
	exit_status = 0
	try
		MacroModel.runleapmacromodel(config_file, logfile, include_energy_sectors)
	catch err
		exit_status = 1
		println(string("Macro exited with an error: Please check the log file '", logfilename,"'"))
		@error err
		if dump_err_stack
			@error stacktrace(catch_backtrace())
		end
	finally
		@info now()
		close(logfile)
		return(exit_status)
	end
end # function run

end #module

ext: jl

max_stars_repo_name: sei-international/LEAPMacro

max_issues_repo_name: sei-international/LEAPMacro

max_forks_repo_name: sei-international/LEAPMacro


text:
##  Copyright 2022 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

abstract type TrustScheme{ T<:Real } end

"""
  Implement for particular parametrizations to the trust weights and propensity model.
"""
bound_ignorance_weights(::TrustScheme{T}, ::PredictionParam{T},
    propensity_params::AbstractVector{T}, order::Int,
    sensitivity::T, treatment::T) where T = # returns (lower, upper)
  error("unimplemented")

"""
  Evaluate the weighing function centered at `treatment` on `point`.
  Possibly (probably) unnormalized for each `treatment`.
"""
weigh_treatment(::TrustScheme{T}, treatment::T, point::T) where T =
  error("unimplemented")

# may not actually be needed, depending on how the general method is implemented
integrate_weights(::TrustScheme{T}, treatment::T) where T =
  error("unimplemented")

"""
  Monte Carlo approximator of confounded ignorance bounds for any `TrustScheme` and propensity.
"""
function generally_bound_ignorance_weights(scheme::TrustScheme{T}, param::PredictionParam{T},
    propensity_params::AbstractVector{T}, order::Int, sensitivity::T, treatment::T;
    n_draws::Int)::NTuple{2,T} where T
  @assert sensitivity > 1
  @assert order in (1, 2)

end

using HypergeometricFunctions, SpecialFunctions

struct BetaTrust{ T } <: TrustScheme{T}
  # (alpha + beta := r) - 2, must be greater than zero
  precision :: T
end

# idea: decouple expectations (except for the first one out of three)
#   from the gamma-bound parts?

function bound_ignorance_weights(scheme::BetaTrust{T}, ::BetaParam{T},
    (prop_alpha, prop_beta)::AbstractVector{T}, order::Int,
    sensitivity::T, treatment::T)::NTuple{2,T} where T
  @assert sensitivity > 1
  @assert order in (1, 2)
  log_gamma = log(sensitivity)
  (hyper_up, hyper_lo) = _₁F₁.(
    scheme.precision * treatment + prop_alpha,
    scheme.precision + prop_alpha + prop_beta,
    (+log_gamma, -log_gamma) ) # \plusminus logGamma
  first_scale = sensitivity^treatment * log_gamma
  first_expectation = (
    ((1-treatment) * prop_alpha - treatment * prop_beta)
    / (prop_alpha + prop_beta + scheme.precision) )
  (first_up, first_lo) = (+1, -1) .* abs(first_scale * first_expectation)
  if order == 1
    (hyper_lo - first_up, hyper_up - first_lo)
  else
    # THIS expectation cannot be negative
    second_expectation = ((
        (prop_alpha^2 + prop_beta^2 + 2prop_alpha*prop_beta
          + prop_alpha + prop_beta - scheme.precision) * treatment^2
        - (2prop_alpha^2 + 2prop_alpha*prop_beta
          + 2prop_alpha - scheme.precision) * treatment
        + prop_alpha * (prop_alpha+1)
      ) / (
        (prop_alpha + prop_beta + scheme.precision)
        * (prop_alpha + prop_beta + scheme.precision + 1) ))
    @assert second_expectation >= 0
    second_scale_up = log_gamma^2 * sensitivity^treatment
    second_scale_lo = T(0)
    (second_up, second_lo) =
      (second_scale_up, second_scale_lo) .* second_expectation
    (both_up, both_lo) = (first_up+second_up, first_lo+second_lo)
    (hyper_lo - both_up, hyper_up - both_lo)
  end
end
# possibly make BetaParam an abstract type grounded by different transformations?

function weigh_treatment(scheme::BetaTrust{T}, treatment::T, point::T)::T where T
  first_part = point ^ (scheme.precision * treatment) # a-1 = (r-2)t
  second_part = (1-point) ^ (scheme.precision * (1-treatment))
  first_part * second_part
end

# do NOT use the general approximator on this. baseline below...
struct BernoulliTrust{ T } <: TrustScheme{T}
  # ways to control *how* stochastic or not?
  # thresholding versus sampling versus "Laplace rule" type debiasing...
end

# what about naively plugging in the continuous densities for the bernoulli propensity value?
function bound_ignorance_weights(scheme::BernoulliTrust{T},
    parametrization::BernoulliParam{T}, (probability,)::AbstractVector{T},
    order::Int, sensitivity::T, treatment::T)::NTuple{2,T} where T
  event = discretize(parametrization, treatment)
  propensity = event ? probability : (1-probability)
  # due to Jesson et al. (2021)
  upper = sensitivity / (propensity*(sensitivity-1) + 1)
  lower = 1 / (propensity*(1-sensitivity) + sensitivity)
  (lower, upper)
end

function bound_ignorance_weights(scheme::BernoulliTrust{T},
    parametrization::BetaParam{T}, (alpha, beta)::AbstractVector{T},
    order::Int, sensitivity::T, treatment::T)::NTuple{2,T} where T
  probability = alpha / (alpha + beta)
  bound_ignorance_weights(scheme, BernoulliParam(T(0)),
    (probability,) |> SVector{1}, order, sensitivity, treatment)
end

max_stars_repo_name: marmarelis/TreatmentCurves

max_issues_repo_name: marmarelis/TreatmentCurves

max_forks_repo_name: marmarelis/TreatmentCurves


text:
using SVM
using Base.Test
srand(1)
iris = readcsv(joinpath(dirname(@__FILE__), "iris.csv"))
X = map(Float64, iris[:, 1:4]')
p, n = size(X)
train = rand(Bool, n)
Y = [species == "setosa" ? 1.0 : -1.0 for species in iris[:, 5]]
model = svm(X[:,train], Y[train])
@test predict(model, X[:, ~train]) == Y[~train]


text:
using SumBassSolver
using Random

Random.seed!(121)
# Generate random data. For implementation input proper data
p = 0.003*rand(50)
q = 0.3*rand(50)
m = 1000 .+ 500*rand(50)
c0 = m .* rand(50) ./ 2
alpha = rand(50)
F = zeros(50)
T = fill(100.0,50)
Pi = 2500

x1,v1 = @time sumBassSolver(p,q,m,c0,alpha,F,T,Pi)
x2,v2 = @time sumBassSolver(p,q,m,c0,alpha,F,T,Pi,localopt=true)
x3,v3 = @time sumBassSolver(p,q,m,c0,alpha,F,T,Pi,alg="SB")
x4,v4 = @time sumBassSolver(p,q,m,c0,alpha,F,T,Pi,alg="SB",localopt=true)

vv1 = sum(alpha .* Bass.(x1,p,q,m,c0))
vv2 = sum(alpha .* Bass.(x2,p,q,m,c0))
vv3 = sum(alpha .* Bass.(x3,p,q,m,c0))
vv4 = sum(alpha .* Bass.(x4,p,q,m,c0))
sum(abs.(x1-x4))

lambdas = Bass_derivative.(x,p,q,m,c0).*alpha

i = findall(x -> abs(x-lambdas[1]) > 0.01,lambdas) # Sites at the boundary
x[i]

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
using Test, YaoExtensions, Yao

@testset "bag" begin
    bag = Bag(ConstGate.T)
    println(bag)
    reg = zero_state(1)
    @test isenabled(bag)
    @test copy(reg) |> bag ≈ copy(reg) |> ConstGate.T
    @test mat(bag) == mat(ConstGate.T)
    @test !isreflexive(bag)
    @test !ishermitian(bag)
    @test isunitary(bag)
    @test occupied_locs(bag) == (1,)

    # disable it
    disable_block!(bag)
    println(bag)
    @test !isenabled(bag)
    @test copy(reg) |> bag ≈ reg
    @test mat(bag) == mat(I2)
    @test isreflexive(bag)
    @test ishermitian(bag)
    @test isunitary(bag)
    @test occupied_locs(bag) == ()

    # enable and setcontent
    enable_block!(bag)
    println(bag)
    @test !ishermitian(bag)
    setcontent!(bag, Z)
    @test ishermitian(bag)
    @test copy(reg) |> X |> bag ≈ ArrayReg([0im, -1.0])
    @test mat(bag) == mat(Z)
    @test ishermitian(bag)
    @test occupied_locs(bag) == (1,)

    # parameters
    b = Bag(Rx(0.5))
    @test nparameters(b) == length(parameters(b)) == 1
    disable_block!(b)
    @test nparameters(b) == length(parameters(b)) == 0
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
using FileIO
using Images

"""
    ImageFilm

A film that produces a resolved image. The resolution of the image is `(resX, resY)`,
the `gain` will be multiplied to scale all the XYZ (and also the resulting RGB) values.

"""
struct ImageFilm{T<:Filter} <: Film
    resX::Int64
    resY::Int64
    gain::Float64
    filter::T
    pixels::Array{Float64,3}
    CIE_table::CIE_Table
end

function ImageFilm(x::Integer, y::Integer, gain::Real, f::Filter, CIE_table::CIE_Table)
    ImageFilm(x, y, gain, f, zeros(4,x,y), CIE_table)
end

"""
    add_sample!(F::ImageFilm, sample::Sample, L::Spectrum, isect::Union{Intersection,Nothing})

Add a ray sample to the `ImageFilm`. The image coordinates in `sample` tell which pixel we
are adding to. The spectral intensity `L` will be added to the pixel, and to nearby pixels
depending on the filter included in the film. `isect` is not used, but it's part of the
`Film` API.

"""
function add_sample!(F::ImageFilm, sample::Sample, L::Spectrum, isect::Union{Intersection,Nothing})
    dimgX = sample.imgX + 0.5
    dimgY = sample.imgY + 0.5
    x0 = ceil(Int64, dimgX - F.filter.xwidth)
    x1 = floor(Int64, dimgX + F.filter.xwidth)
    y0 = ceil(Int64, dimgY - F.filter.ywidth)
    y1 = floor(Int64, dimgY + F.filter.ywidth)
    x0 = max(x0, 1)
    x1 = min(x1, F.resX)
    y0 = max(y0, 1)
    y1 = min(y1, F.resY)
    if (x1-x0) < 0 || (y1-y0) < 0
        error("This shouldn't happen? Sample: $sample")
    end
    x, y, z = to_XYZ(L, F.CIE_table)
    for j = y0:y1
        for i = x0:x1
            w = evaluate(F.filter, i-dimgX, j-dimgY)
            F.pixels[1,i,j] += w*x
            F.pixels[2,i,j] += w*y
            F.pixels[3,i,j] += w*z
            F.pixels[4,i,j] += w
        end
    end
end

"""
    reset!(F::ImageFilm)

Reset all the pixel values of the `ImageFilm` to zero. This is mainly useful if you want to
use the same `ImageFilm` in multiple computations, writing out the image and resetting in
between.

"""
function reset!(F::ImageFilm)
    F.pixels .= 0.0
    nothing
end

"""
    write_image(F::ImageFilm, fname::String="test.png")

Write the given image into a file determined by the filename.
"""
function write_image(F::ImageFilm, fname::String="test.png")
    data = zeros(Float64, (3, size(F.pixels,2), size(F.pixels, 3)))
    @inbounds for j = 1:size(F.pixels,3)
        for i = 1:size(F.pixels,2)
            x,y,z,w = F.pixels[:,i,j]
            W = F.gain / w
            r, g, b = XYZtoRGB(W*x, W*y, W*z)
            data[1,i,j] = r
            data[2,i,j] = g
            data[3,i,j] = b
        end
    end
    img = colorview(RGB, map(clamp01nan, data))
    save(fname, img)
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
function nbytesToType(nbytes::Cssize_t)
    types = Nothing, Nothing;
    normalizing_factors = Nothing, Nothing;
    if nbytes == 1
        types = UInt8, Int8;
        normalizing_factors = 2^8 - 1, 2^7;
    elseif nbytes == 2
        types = UInt16, Int16;
        normalizing_factors = 2^16 - 1, 2^15;
    elseif nbytes == 4
        types = UInt32, Int32;
        normalizing_factors = 2^32 - 1, 2^31;
    elseif nbytes == 8
        types = UInt64, Int64;
        normalizing_factors = 2^64 - 1, 2^63;
    end
    return types..., normalizing_factors...;
end

# --- Pretty printing functions --- #

# To print fancy message with different colors with Tx and Rx
function customPrint(str, handler;style...)
	msglines = split(chomp(str), '\n');
    length(msglines) > 1 ? symb = "┌" : symb = "[";
	printstyled(symb, handler, ": ";style...);
	println(msglines[1]);
	for i in 2:length(msglines)
		(i == length(msglines)) ? symb = "└ " : symb = "|";
		printstyled(symb;style...);
		println(msglines[i]);
	end
end

"""
    @infoPluto TRX str

Prints custom info messages for either Rx or Tx.

# Example
- `@infoPluto :RX "msg"` to print a warning concerning the Rx part.
- `@infoPluto :TX "msg"` to print a warning concerning the Tx part.
"""

macro infoPluto(TRX, str)
    quote
        TRX = $(TRX);
        if TRX == :RX
            customPrint($(esc(str)), "Rx";bold=true,color=:light_blue)
        else TRX == :TX
            customPrint($(esc(str)), "Tx";bold=true,color=:light_blue)
        end
    end
end

"""
    @warnPluto TRX str

Prints custom warnings for either Rx or Tx.

# Example
- `@warnPluto :RX "msg"` to print a warning concerning the Rx part.
- `@warnPluto :TX "msg"` to print a warning concerning the Tx part.
"""
macro warnPluto(TRX, str)
    quote
        TRX = $(TRX);
        if TRX == :RX
            customPrint($(esc(str)), "Rx Warning";bold=true,color=:light_yellow)
        else TRX == :TX
            customPrint($(esc(str)), "Tx Warning";bold=true,color=:light_yellow)
        end
    end
end

"""
Returns the list of the available backends in Pluto, as a Vector of Strings
str = getBackends()
    str = ["usb";"ip";"xml"]
"""
function getBackends()
    # We do not have directly the backend names, but we can get the count, and then iterate to get the names
    nbBackend = C_iio_get_backends_count()
    str = String[]
    for n ∈ 1 : nbBackend
        # --- Try to get  backend name
        s = C_iio_get_backend(UInt32(n-1))
        # This should not be ""
        if !isempty(s)
            # --- Try to load backend to see if everything is fine
            flag = C_iio_has_backend(s)
            if flag == true
                # --- Store the backend
                push!(str,s)
            end
        end
    end
    return str
end

# Export udev rules
set_udev_rules = libIIO_jl.set_udev_rules


text:
#####
##### Transformation of the Chebyshev polynomials — generic code
#####

export univariate_basis, BoundedLinear, InfRational, SemiInfRational

####
#### simplified transformation API
####

"""
$(TYPEDEF)

An abstract type for univariate transformations. Transformations are not required to be
subtypes, this just documents the interface they need to support:

- [`to_domain`](@ref)

- [`from_domain`](@ref)
"""
abstract type UnivariateTransformation end

"""
`$(FUNCTIONNAME)(transformation, parent, x)`

!!! FIXME
    document, especially differentiability requirements at infinite endpoints
"""
function to_domain end

"""
`$(FUNCTIONNAME)(transformation, parent, x)`

!!! FIXME
    document, especially differentiability requirements at infinite endpoints
"""
function from_domain end

####
#### univariate basis
####

struct UnivariateBasis{P,T} <: FunctionBasis
    parent::P
    transformation::T
end

function Base.show(io::IO, univariate_basis::UnivariateBasis)
    @unpack parent, transformation = univariate_basis
    print(io, parent, "\n  domain ", transformation)
end

"""
$(SIGNATURES)

Create a univariate basis from `univariate_family`, with the specified `grid_kind` and
dimension `N`, transforming the domain with `transformation`.

`parent` is a univariate basis, `transformation` is a univariate transformation (supporting
the interface described by [`UnivariateTransformation`](@ref), but not necessarily a
subtype). Univariate bases support [`SpectralKit.gridpoint`](@ref).

# Example

The following is a basis with 10 transformed Chebyshev polynomials of the first kind on
``(3,∞)``, with equal amounts of nodes on both sides of `7 = 3 + 4` and an interior grid:

```jldoctest
julia> basis = univariate_basis(Chebyshev, InteriorGrid(), 10, SemiInfRational(3.0, 4.0))
Chebyshev polynomials (1st kind), interior grid, dimension: 10
  domain (3.0,∞) [rational transformation with scale 4.0]

julia> dimension(basis)
10

julia> domain(basis)
(3.0, Inf)
```
"""
function univariate_basis(univariate_family, grid_kind::AbstractGrid, N::Integer,
                          transformation)
    UnivariateBasis(univariate_family(grid_kind, Int(N)), transformation)
end

@inline dimension(basis::UnivariateBasis) = dimension(basis.parent)

function domain(basis::UnivariateBasis)
    @unpack parent, transformation = basis
    from_domain.(Ref(transformation), Ref(parent), domain(parent))
end

function basis_at(basis::UnivariateBasis, x::Real)
    @unpack parent, transformation = basis
    basis_at(parent, to_domain(transformation, parent, x))
end

function gridpoint(::Type{T}, basis::UnivariateBasis, i) where T
    @unpack parent, transformation = basis
    from_domain(transformation, parent, gridpoint(T, parent, i))
end

####
#### transformations
####

###
### bounded linear
###

struct BoundedLinear{T <: Real} <: UnivariateTransformation
    "Midpoint `m`."
    m::T
    "Scale `s`."
    s::T
    function BoundedLinear(a::T, b::T) where {T <: Real}
        @argcheck isfinite(a) && isfinite(b) DomainError
        s = (b - a) / 2
        m = (a + b) / 2
        @argcheck s > 0 DomainError((; a, b), "Need `a < b`.")
        m, s = promote(m, s)
        new{typeof(m)}(m, s)
    end
end

function Base.show(io::IO, transformation::BoundedLinear)
    @unpack m, s = transformation
    print(io, "(", m - s, ",", m + s, ") [linear transformation]")
end

"""
$(TYPEDEF)

Transform `x ∈ (-1,1)` to `y ∈ (a, b)`, using ``y = x ⋅ s + m``.

`m` and `s` are calculated and checked by the constructor; `a < b` is enforced.
"""
BoundedLinear(a::Real, b::Real) = BoundedLinear(promote(a, b)...)

function from_domain(T::BoundedLinear, ::Chebyshev, x::Real)
    @unpack m, s = T
    x * s + m
end

function to_domain(T::BoundedLinear, ::Chebyshev, y::Real)
    @unpack m, s = T
    (y - m) / s
end

###
### semi-infinite interval
###

struct SemiInfRational{T <: Real} <: UnivariateTransformation
    "The finite endpoint `A`."
    A::T
    "Scale factor `L ≠ 0`."
    L::T
    function SemiInfRational(A::T, L::T) where {T <: Real}
        @argcheck isfinite(A) DomainError
        @argcheck isfinite(L) && L ≠ 0 DomainError
        new{T}(A, L)
    end
end

function Base.show(io::IO, transformation::SemiInfRational)
    @unpack A, L = transformation
    if L > 0
        D = "($A,∞)"
    else
        D = "(-∞,A)"
    end
    print(io, D, " [rational transformation with scale ", L, "]")
end

"""
$(SIGNATURES)

`[-1,1]` transformed to the domain `[A, Inf)` (when `L > 0`) or `(-Inf,A]`
(when `L < 0`) using ``y = A + L ⋅ (1 + x) / (1 - x)``.

When used with Chebyshev polynomials, also known as a “rational Chebyshev” basis.

# Example mappings

- ``-1/2 ↦ A + L / 3``
- ``0 ↦ A + L``
- ``1/2 ↦ A + 3 ⋅ L``
"""
SemiInfRational(A::Real, L::Real) = SemiInfRational(promote(A, L)...)

from_domain(T::SemiInfRational, ::Chebyshev, x) = T.A + T.L * (1 + x) / (1 - x)

function to_domain(T::SemiInfRational, ::Chebyshev, y)
    @unpack A, L = T
    z = y - A
    x = (z - L) / (z + L)
    if (y == Inf && L > 0) || (y == -Inf && L < 0)
        one(x)
    else
        x
    end
end

###
### infinite interval
###

struct InfRational{T <: Real} <: UnivariateTransformation
    "The center `A`."
    A::T
    "Scale factor `L > 0`."
    L::T
    function InfRational(A::T, L::T) where {T <: Real}
        @argcheck isfinite(A) DomainError
        @argcheck isfinite(L) && L > 0 DomainError
        new{T}(A, L)
    end
end

function Base.show(io::IO, transformation::InfRational)
    @unpack A, L = transformation
    print(io, "(-∞,∞) [rational transformation with center ", A, ", scale ", L, "]")
end

"""
$(SIGNATURES)

Chebyshev polynomials transformed to the domain `(-Inf, Inf)`
using ``y = A + L ⋅ x / √(1 - x^2)``, with `L > 0`.

# Example mappings

- ``0 ↦ A``
- ``±0.5 ↦ A ± L / √3``
"""
InfRational(A::Real, L::Real) = InfRational(promote(A, L)...)

from_domain(T::InfRational, ::Chebyshev, x::Real) = T.A + T.L * x / √(1 - abs2(x))

function to_domain(T::InfRational, ::Chebyshev, y::Real)
    @unpack A, L = T
    z = y - A
    x = z / hypot(L, z)
    if isinf(y)
        y > 0 ? one(x) : -one(x)
    else
        x
    end
end


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

"""
    abstract type BAT.InitvalAlgorithm

Abstract type for BAT initial/starting value generation algorithms.

Many algorithms in BAT, like MCMC and optimization, need initial/starting
values.
"""
abstract type InitvalAlgorithm end
export InitvalAlgorithm

"""
    bat_initval(
        [rng::AbstractRNG,]
        target::BAT.AnyDensityLike,
        [algorithm::BAT.InitvalAlgorithm],
    )::V

    bat_initval(
        [rng::AbstractRNG,]
        target::BAT.AnyDensityLike,
        n::Integer,
        [algorithm::BAT.InitvalAlgorithm],
    )::AbstractVector{<:V}

Generate one or `n` random initial/starting value(s) suitable for
`target`.

Assuming the variates of `target` are of type `T`, returns a NamedTuple of
the shape

```julia
(result = X::AbstractVector{T}, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_initval`, add methods like

    ```julia
    bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, algorithm::InitvalAlgorithm; kwargs...)
    bat_initval_impl(rng::AbstractRNG, target::AnyDensityLike, n::Integer, algorithm::InitvalAlgorithm; kwargs...)
    ```

    to `bat_initval_impl` instead.
"""
function bat_initval end
export bat_initval

function bat_initval_impl end

@inline function bat_initval(rng::AbstractRNG, target::AnyDensityLike, algorithm::InitvalAlgorithm; kwargs...)
    r = bat_initval_impl(rng, target, algorithm; kwargs...)
    result_with_args(r, (rng = rng, algorithm = algorithm), kwargs)
end

@inline function bat_initval(target::AnyDensityLike; kwargs...)
    rng = bat_default_withinfo(bat_initval, Val(:rng), target)
    algorithm = bat_default_withinfo(bat_initval, Val(:algorithm), target)
    bat_initval(rng, target, algorithm; kwargs...)
end

@inline function bat_initval(target::AnyDensityLike, algorithm::InitvalAlgorithm; kwargs...)
    rng = bat_default_withinfo(bat_initval, Val(:rng), target)
    bat_initval(rng, target, algorithm; kwargs...)
end

@inline function bat_initval(rng::AbstractRNG, target::AnyDensityLike; kwargs...)
    algorithm = bat_default_withinfo(bat_initval, Val(:algorithm), target)
    bat_initval(rng, target, algorithm; kwargs...)
end

@inline function bat_initval(rng::AbstractRNG, target::AnyDensityLike, n::Integer, algorithm::InitvalAlgorithm; kwargs...)
    r = bat_initval_impl(rng, target, n, algorithm; kwargs...)
    result_with_args(r, (rng = rng, algorithm = algorithm), kwargs)
end

@inline function bat_initval(target::AnyDensityLike, n::Integer; kwargs...)
    rng = bat_default_withinfo(bat_initval, Val(:rng), target)
    algorithm = bat_default_withinfo(bat_initval, Val(:algorithm), target)
    bat_initval(rng, target, n, algorithm; kwargs...)
end

@inline function bat_initval(target::AnyDensityLike, n::Integer, algorithm::InitvalAlgorithm; kwargs...)
    rng = bat_default_withinfo(bat_initval, Val(:rng), target)
    bat_initval(rng, target, n, algorithm; kwargs...)
end

@inline function bat_initval(rng::AbstractRNG, target::AnyDensityLike, n::Integer; kwargs...)
    algorithm = bat_default_withinfo(bat_initval, Val(:algorithm), target)
    bat_initval(rng, target, n, algorithm; kwargs...)
end

function argchoice_msg(::typeof(bat_initval), ::Val{:rng}, x::AbstractRNG)
    "Initializing new RNG of type $(typeof(x))"
end

function argchoice_msg(::typeof(bat_initval), ::Val{:algorithm}, x::InitvalAlgorithm)
    "Using initval algorithm $x"
end


text:
function testephemeris(denum)
    ephem = SPK("$path/de$denum.bsp")
    if isfile("$path/testpo.$denum")
        lines = open(readlines, "$path/testpo.$denum")
    else
        error("Test file 'testpo.$denum' not found.")
    end

    start = findfirst(lines .== "EOT\n") + 1
    for l in lines[start:end]
        de, date, jd, target, center, index, value = split(l)
        target = parse(Int, target)
        center = parse(Int, center)
        index = parse(Int, index)
        value = float(value)
        jd = float(jd)

        if target in 14:15
            continue
        end

        try
            tr = teststate(ephem, jd, target)
            cr = teststate(ephem, jd, center)
            # From km/s to km/day
            tr[4:6] *= 86400
            cr[4:6] *= 86400
            # To AU and AU/day
            r = (tr - cr)/AU

            @test isapprox(r[index], value, atol=1e-13)
        # The test file for DE405 contains a wider range of dates than the SPK kernel provides.
        # Ignore those.
        catch err
            if isa(err, JPLEphemeris.OutOfRangeError)
                continue
            else
                rethrow(err)
            end
        end
    end
end

function teststate(kernel, tbd, target)
    if target == 3
        rv1 = state(kernel, 0, 3, tbd)
        rv2 = state(kernel, 3, 399, tbd)
        rv = rv1+rv2
    elseif target == 10
        rv1 = state(kernel, 0, 3, tbd)
        rv2 = state(kernel, 3, 301, tbd)
        rv = rv1+rv2
    elseif target == 11
        rv = state(kernel, 0, 10, tbd)
    elseif target == 12
        rv = zeros(6)
    elseif target == 13
        rv = state(kernel, 0, 3, tbd)
    else
        rv = state(kernel, 0, target, tbd)
    end
    return rv
end

# Run the JPL testsuite for every installed ephemeris.
@testset "Kernels" begin
    @testset for denum in (430, 405)
        testephemeris(denum)
    end
end


text:
module BossStarField

using ..Ahorn, Maple

placements = Maple.BossStarField

function Ahorn.canFgBg(effect::Maple.BossStarField)
    return true, true
end

end


text:
using Test
using LyapunovExponents
using LyapunovExponents: PhaseRelaxer,
    get_integrator, get_tangent_integrator, van_der_pol

@time @testset "Integrator time continuity" begin
    prob = van_der_pol().prob
    relaxer = PhaseRelaxer(get_integrator(prob.phase_prob),
                           prob.t_tran)
    dt = 1.0
    step!(relaxer.integrator, dt, true)
    @test relaxer.integrator.t > prob.phase_prob.tspan[1]

    integrator = get_tangent_integrator(prob, relaxer)
    @test integrator.t == relaxer.integrator.t

    step!(relaxer.integrator, dt, true)
    step!(integrator, dt, true)
    @test integrator.t == relaxer.integrator.t

    @test integrator.u[:, 1] ≈ relaxer.integrator.u rtol = 0.05
end


text:
using MLStyle
using StatsFuns

export @measure

# A fold over ASTs. Example usage in `replaceX`
function foldast(leaf, branch; kwargs...)
    function go(ast)
        MLStyle.@match ast begin
            Expr(head, args...) => branch(head, map(go, args); kwargs...)
            x                   => leaf(x; kwargs...)
        end
    end

    return go
end

# Walk through expr, replacing every occurrence of :X with newX
function replaceX(newX, expr)
    leaf(s) = (s == :X) ? newX : s
    branch(head, newargs) = Expr(head, newargs...)
    foldast(leaf, branch)(expr)
end

# from https://thautwarm.github.io/MLStyle.jl/latest/tutorials/capture.html
function capture(template, ex, action)
    let template = Expr(:quote, template)
        quote
            @match $ex begin
                $template => $action
                _         => nothing
            end
        end
    end
end

macro capture(template, ex, action)
    capture(template, ex, action) |> esc
end

function _measure(expr)
    @capture $rel($μ($(p...)), $base) expr begin
        q = quote
            struct $μ{P} <: MeasureTheory.AbstractMeasure
                par :: P
            end

            function $μ(nt::NamedTuple)
                P = typeof(nt)
                return $μ{P}(nt)
            end

            $μ(;kwargs...) = $μ((;kwargs...))

            function basemeasure(μ::$μ{P}) where {P}
                X = sampletype($μ{P})
                $base
            end

            $μ($(p...)) = $μ(;$(p...))

            # ((::$μ{P} ≪ ::typeof($base) ) where {P})  = true
        end

        # if rel == (:≃)
        #     push!(q.args, :((::typeof($base) ≪ ::$μ{P}) where {P} = true))
        # end

        return q
    end
end

"""
    @measure <declaration>

The <declaration> gives a measure and its default parameters, and specifies
its relation to its base measure. For example,

    @measure Normal(μ,σ) ≃ Lebesgue{X}

declares the `Normal` is a measure with default parameters `μ and σ`, and it is
equivalent to its base measure, which is `Lebesgue{X}`

You can see the generated code like this:

    julia> MacroTools.prettify(@macroexpand @measure Normal(μ,σ) ≃ Lebesgue{X})
    quote
        struct Normal{P, X} <: AbstractMeasure
            par::P
        end
        function Normal(nt::NamedTuple)
            P = typeof(nt)
            return Normal{P, eltype(Normal{P})}
        end
        Normal(; kwargs...) = Normal((; kwargs...))
        (basemeasure(μ::Normal{P, X}) where {P, X}) = Lebesgue{X}
        Normal(μ, σ) = Normal(; Any[:μ, :σ])
        ((:≪)(::Normal{P, X}, ::Lebesgue{X}) where {P, X}) = true
        ((:≪)(::Lebesgue{X}, ::Normal{P, X}) where {P, X}) = true
    end

Note that the `eltype` function needs to be defined separately by the user.
"""
macro measure(expr)
    esc(_measure(expr))
end

# (@macroexpand @measure Normal(μ,σ) ≃ (1/sqrt2π) * Lebesgue(X)) |> MacroTools.prettify


text:
# Note to self -- ZygoteRules doesn't have forward, only @adjoint, so it's no help.
# Later: it now has pullback, so that should work? Once registered... all copied in.

using .Zygote
using .Zygote: @adjoint, _zero, forward

#===== mapcols, maprows, MapCols =====#

@adjoint _mapcols(map::Function, f::Function, M::AbstractMatrix, args...) =
    ∇mapcols(map, map(col -> forward(x -> surevec(f(x, args...)), col), eachcol(M)), args)

@adjoint maprows(f::Function, M::AbstractMatrix, args...) =
    ∇maprows(map(row -> forward(x -> surevec(f(x, args...)), row), eachrow(M)), args)

@adjoint _MapCols(map::Function, f::Function, M::Matrix, dval, args...) =
    ∇MapCols(map, f, M, dval, args...)

#===== TensorCast =====#

@adjoint TensorCast.sliceview(A::AbstractArray, code::Tuple) =
    TensorCast.sliceview(A, code), Δ -> (TensorCast.glue(Δ, code), nothing)

@adjoint TensorCast.red_glue(A::AbstractArray, code::Tuple) =
    TensorCast.red_glue(A, code), Δ -> (TensorCast.sliceview(Δ, code), nothing)

@adjoint TensorCast.copy_glue(A::AbstractArray, code::Tuple) =
    TensorCast.copy_glue(A, code), Δ -> (TensorCast.sliceview(Δ, code), nothing)

#===== JuliennedArrays =====#

@adjoint JuliennedArrays.Slices(whole, along...) =
    Slices(whole, along...), Δ -> (Align(Δ, along...), map(_->nothing, along)...)

@adjoint JuliennedArrays.Align(whole, along...) =
    Align(whole, along...), Δ -> (Slices(Δ, along...), map(_->nothing, along)...)

#===== Misc Base =====#

@adjoint Base.reduce(::typeof(hcat), V::AbstractVector{<:AbstractVector}) =
    reduce(hcat, V), dV -> (nothing, collect(eachcol(dV)),)

#===== Misc experiments =====#

@adjoint gluecol(V::AbstractVector) =
    gluecol(V), ΔM -> (collect(eachcol(ΔM)),) # does work!

@adjoint function collecteachcol(x)
    collecteachcol(x), dy -> begin
        dx = _zero(x)
        foreach(copyto!, collecteachcol(dx), dy)
        (dx,)
    end
end


text:
export ilan_benchmark

using IterativeSolvers
using LinearAlgebra
using Random
using Statistics

"""
    ilan_benchmark()
"""
ilan_benchmark(nep::NEP;params...)=ilan_benchmark(ComplexF64,nep;params...)
function ilan_benchmark(
    ::Type{T},
    nep::NEP;
    orthmethod::Type{T_orth}=DGKS,
    maxit=30,
    linsolvercreator=DefaultLinSolverCreator(),
    tol=eps(real(T))*10000,
    neigs=6,
    errmeasure::ErrmeasureType = DefaultErrmeasure,
    σ=zero(T),
    γ=one(T),
    v=randn(real(T),size(nep,1)),
    displaylevel=0,
    check_error_every=1,
    proj_solve=false,
    inner_solver_method=DefaultInnerSolver())where{T<:Number,T_orth<:IterativeSolvers.OrthogonalizationMethod}

    # Ensure types σ and v are of type T
    σ=T(σ)
    v=Array{T,1}(v)

    n = size(nep,1);
    m = maxit;

    # initialization
    V=zeros(T,n,m+1)
    Q=zeros(T,n,m+1)
    Qp=zeros(T,n,m+1)
    Qn=zeros(T,n,m+1)
    Z=zeros(T,n,m+1)
    H=zeros(T,m+1,m)
    HH=zeros(T,m+1,m)
    ω=zeros(T,m+1)
    a=Vector{T}(γ.^(0:m)); a[1]=zero(T); # TODO
    local M0inv::LinSolver = create_linsolver(linsolvercreator,nep,σ)
    err=ones(m,m);
    λ=zeros(T,m+1);

    # precompute the symmetrizer coefficients
    G=zeros(m+1,m+1);
    for i=1:m+1 G[i,1]=1/i end
    for j=1:m
        for i=1:m+1
            G[i,j+1]=(G[i,j]*j)/(i+j);
        end
    end

    # getting matrices and functions
    fv=get_fv(nep); p=length(fv);    Av=get_Av(nep)

    # precompute derivatives and FDH matrices
    SS=diagm(0 => σ*ones(T,2m+2)) + diagm(-1 => (1:2m+1))
    fD=zeros(T,2*m+2,p)
    for t=1:p fD[:,t]=fv[t](SS)[:,1] end
    FDH=Vector{Array{T,2}}(undef,p)
    for t=1:p FDH[t]=zeros(T,m+1,m+1)
        for i=1:m+1 for j=1:m+1
            FDH[t][i,j]=fD[i+j,t];
        end end
    end

    # setting initial step
    Q[:,1]=v/norm(v)
    ω[1]=Q[:,1]⋅compute_Mlincomb(nep,0,hcat(Q[:,1],Q[:,1]),[0,1]);
    V[:,1]=Q[:,1];

    k=1; conv_eig=0;
    local pnep::NEP;
    if (proj_solve)
        pnep=create_proj_NEP(nep);
    end

    while (k <= m) && (conv_eig<neigs)
        if (displaylevel>0) && ((rem(k,check_error_every)==0) || (k==m))
            println("Iteration:",k, " conveig:",conv_eig)
        end

        broadcast!(/,view(Qn,:,2:k+1),view(Q,:,1:k),(1:k)')
        Qn[:,1] = compute_Mlincomb!(nep,σ,Qn[:,1:k+1],a[1:k+1]);
        Qn[:,1] = -lin_solve(M0inv,Qn[:,1]);

        # B-multiplication
        Z=0*Z;
        for t=1:p
            Z[:,1:k+1]=Z[:,1:k+1]+Av[t]*Qn[:,1:k+1]*(G[1:k+1,1:k+1].*FDH[t][1:k+1,1:k+1]);
        end

        # orthogonalization (three terms recurrence)
        α=sum(sum(conj(Z).*Q,dims=1))
        if k>1 β=sum(sum(conj(Z).*Qp,dims=1)) end
        η=sum(sum(conj(Z).*Qn,dims=1))

        H[k,k]=α/ω[k]
        if k>1 H[k-1,k]=β/ω[k-1] end
        Qn=Qn-H[k,k]*Q;
        if k>1 Qn=Qn-H[k-1,k]*Qp end

        H[k+1,k]=norm(Qn);
        Qn=Qn/H[k+1,k];

        ω[k+1]=η-2*α*H[k,k]+ω[k]*H[k,k]^2;
        if k>1
            ω[k+1]=ω[k+1]-2*β*H[k-1,k]+ω[k-1]*H[k-1,k]^2;
        end
        ω[k+1]=ω[k+1]/H[k+1,k]^2;
        V[:,k+1]=Qn[:,1];

        orthogonalize_and_normalize!(view(V,:,1:k),view(V,:,k+1), view(HH,1:k,k), orthmethod)

        k=k+1;
        # shift the vectors
        Qp=Q;   Q=Qn; Qn=0*Qn;

    end
    k=k-1

    return V[:,1:k+1], H[1:k,1:k-1], ω[1:k], HH[1:k,1:k]
end


text:
@testset "MurmurHash3" begin
      tests = [
      Dict{String,Any}("x86_32"=>"00000000","text"=>""),
      Dict{String,Any}("x86_32"=>"087fcd5c","seed"=>42,"text"=>""),
      Dict{String,Any}("x86_32"=>"1cb30d2d","text"=>"!eDoZj \$XVe>vtaxn<x^ qU]m([aWSI^Hy>@E#hX}v mbkz#gB"),
      Dict{String,Any}("x86_32"=>"e6aa80e1","seed"=>10,"text"=>"o.SH}f]tVS\$deSKiOgE?#V %EP&]V!oyjlq}{BO)\$cz&fRGPrGUtz@b]{C"),
      Dict{String,Any}("x86_32"=>"100949f7","text"=>"[)<{vbjRgcp"),
      Dict{String,Any}("x86_32"=>"79b365c1","text"=>"R"),
      Dict{String,Any}("x86_32"=>"e70207f0","seed"=>5489776,"text"=>"mB"),
      Dict{String,Any}("x86_32"=>"c4b773e5","seed"=>111117,"text"=>"mExonGybCN{n# iVmrVI!?ahk]q&}BbHyx"),
      Dict{String,Any}("x86_32"=>"92bbd066","text"=>"TqDa}YvgP!K&P\$vsbeGt<Q{zVuBHQ>h}tBxvVlYVO\$LAgxNIik{)WeJtX]i?yom!.KoNB*sOUJv/.pF,lncuM@EXBTgw%VuY>\$!l!TFGe?HEO>NFiK,RSS%o#)lhy&Mw{QFGN![qKcM!eSXxfQ&eYeTtq[m<h#rvwuaNFTv e@kXzD##[%.JDYlKhQzrGV{/bB)^dCT&GDum#V{cVgAn)xjahnT[E\$<)yyElAhV/imQ!P%JL\$&tB%oaj*L\$<qzsD\$*oKT^/ykA([\$(sI}G{blBT>SSo UnOlclQkAJ(@FLXP&w/EKY?Kn}ndbY\$!#IC%IlcWF]^UX*iVejGrzX%ox@b{ &KUm>q>dP,k\$(&mmE,f.Kp/R[.FnjDPJ?FhniNuAi<o%*m./BjypWIaao*&@<kGICn!qP(uTjavSxzKfaT(H]!NmRWgejf?#q%l&ulr oxLv])NsNM>Q#bEAXy[lZXod?Yeu][SYaO DZ^doDLmgrCuBBWy"),
      Dict{String,Any}("x86_32"=>"8ba3f73c","text"=>"#Xu"),
      Dict{String,Any}("x86_32"=>"fb73894d","seed"=>1,"text"=>"MMBh*"),
      Dict{String,Any}("x86_32"=>"c83626b3","text"=>"QDJ{gyEAJMbhd\$,iRU(!uv>zCpzIDjcel F#/?/Tux?nX{^n*?a"),
      Dict{String,Any}("x86_32"=>"91032392","text"=>"GeBcX^k?C,]d*nx##RL*#*lZ{ &b{.I?Z]#qMyMiMT"),
      Dict{String,Any}("x86_32"=>"1b9cb3e5","seed"=>160753,"text"=>"KA.f#)K?#%by#G(w?ag^at{jbW)Q(QJtvX\$Gs{@UMWhhzYdX# A{uPrvm#"),
      Dict{String,Any}("x86_32"=>"54c379d4","seed"=>2,"text"=>"%w.I#vFp"),
      Dict{String,Any}("x86_32"=>"950a02a4","seed"=>0,"text"=>"th!S?L*A[VUtZZ)[%Zw/&!ZO!MsJDcFjV/uGb^yuZLZ\$)mlf}#N\$<b,!k t@R"),
      Dict{String,Any}("x86_32"=>"c7d8827f","seed"=>999983,"text"=>"IqP ]t (OrK"),
      Dict{String,Any}("x86_32"=>"704a49be","text"=>"PKqG47bFcNTKKTE"),
      Dict{String,Any}("x86_32"=>"52042560","text"=>"Qs7i2oN0kr10dOf")]
      for test in tests
            if "seed" in keys(test)
                  x86_32 = lpad(string(Bitcoin.Murmur3.hash32(test["text"], test["seed"]), base=16), 8, "0")
            else
                  x86_32 = lpad(string(Bitcoin.Murmur3.hash32(test["text"]), base=16), 8, "0")
            end
            @test test["x86_32"] == x86_32
      end
end


text:
"""
    This example shows how to use Term's layout functionality.

When creating structured terminal output, you want to have control
over what goes well. Maybe you want to print two panels, side by side,
or you want to stack them vertically. Maybe you want to have some space
between content or use line to delimit sections of content that should be
separate. Term.jl has got your back.
"""

"""
Let's start with stacking.
You can use `*` to horizzontally stack content, and `/` to stack it
vertically.
"""

using Term

# with strings
s1 = @red "one string"
s2 = @green "and another"

println(s1 * " " * s2)
println(s1 / s2)

# with panels (and text boxes)
p1 = Panel("content [blue]content[/blue] "^10; width = 40)
tb1 = TextBox("this is a text boox! "^5; width = 66)

println(p1 * tb1)
println(p1 / tb1)

# you can stack multiple things too, and mix types
println(s1 / p1 / s2 / tb1)

"""
Now, stacking things is great, but want if you want some space between them?
`Spacer` lets you create an empty space with set width and height
"""
space = Spacer(20, p1.measure.h)  # create empty space with same height as p1!

# horizzontally stack everything and print!
top = p1 * space * p1
println(top)

"""
you can use `hLine` and `vLine` to create styled lines to separate pieces of content
"""

print("\n\n\n")
vline = vLine(p1.measure.h; style = "blue")  # create a vertical line as high as the panel
top = p1 * vline * space * vline * p1
line = hLine(top.measure.w, "some space"; style = "red", box = :DOUBLE)

# vertically stack everything and print nested in a panel
println(
    Panel(
        top / line / top;
        title = "layout is great!",
        title_style = "bold red",
        style = "red dim",
        title_justify = :left,
        fit=true
    ),
)


text:
# This code is unused
# Code written in Symata that is actually loaded is in ./symsrc

# Some code written in Symata itself is here.
#
# There are more examples in the test directory. eg. io_test.jl
# These are examples; many are not practical.
#
# It is tested in the  file code_in_Symata_test.jl.

# Examples of using down values

# With the present implmentation, this is slower than
# writing an apprule directly in Julia. Patterns are not optimized
unset_attribute(:Sin,:Protected)
unset_attribute(:Tan,:Protected)
unset_attribute(:Cot,:Protected)
unset_attribute(:Cos,:Protected)
unset_attribute(:Sec,:Protected)
unset_attribute(:Csc,:Protected)

# The multiplication rules won't work generally until AC matching is implemened
# Some are not evaled far enough.
# Unfix is a workaround for bugs that prevent evaluation
@sym  Sin(-1*x_) := -1 * Sin(x)
@sym  Tan(-1*x_) := -1 * Tan(x)
@sym  Tan(ArcCos(x_)) := Unfix(Sqrt(1-x^2)/x)
@sym  Tan(ArcSin(x_)) := Unfix(x/Sqrt(1-x^2))
@sym  Sin(ArcSin(x_)) := x
@sym  Sin(ArcCos(x_)) := Unfix((1-x^2)^(1/2))
@sym  Cos(-1*x_) := Cos(x)

@sym  Power(Cos(x_),-1) ^= Sec(x)
@sym  Power(Sec(x_),-1) ^= Unfix(Cos(x)) # need to fix bug that requires Unfix!
@sym  Power(Sin(x_),-1) ^= Csc(x)
@sym  Power(Csc(x_),-1) ^= Sin(x)
@sym  Power(Tan(x_),-1) ^= Cot(x)
@sym  Power(Cot(x_),-1) ^= Tan(x)

set_attribute(:Sin,:Protected)
set_attribute(:Cos,:Protected)
set_attribute(:Cot,:Protected)
set_attribute(:Tan,:Protected)
set_attribute(:Sec,:Protected)
set_attribute(:Csc,:Protected)

# the symbol x should be local, anyway.
@sym ClearAll(x)

# These are only for testing downrules.
const directed_infinitym1 = setfixed(mxpr(:DirectedInfinity,-1))
unprotect(:Log)
@sym Log(1) := 0
# This is slow because the julia expression is parsed every time
@sym Log(0) := DirectedInfinity(-1)
protect(:Log)

@sjdoc Log "
Log(x) represents the natural logarithm of x.
Log(b,x) represents the base \"b\" logarithm of x.
"

unprotect(:Zeta)
@sym Zeta(1) := ComplexInfinity
protect(:Zeta)

### These would work for symata system code, but they
### add a few seconds to startup because fo jit.

### TODO: find a way to autoload this code

unprotect(:StringQ)
@sym StringQ = MatchQ(_AbstractString)
protect(:StringQ)

unprotect(:ListQ)
@sym ListQ = MatchQ(_List)
protect(:ListQ)


text:
module Compat

include("compatmacro.jl")

# https://github.com/JuliaLang/julia/pull/29679
if VERSION < v"1.1.0-DEV.472"
    export isnothing
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end

# https://github.com/JuliaLang/julia/pull/29749
@static if VERSION < v"1.1.0-DEV.792"
    export eachrow, eachcol, eachslice
    eachrow(A::AbstractVecOrMat) = (view(A, i, :) for i in axes(A, 1))
    eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2))
    @inline function eachslice(A::AbstractArray; dims)
        length(dims) == 1 || throw(ArgumentError("only single dimensions are supported"))
        dim = first(dims)
        dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
        idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
        return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
    end
end

function rangeargcheck(;step=nothing, length=nothing, kwargs...)
    if step===nothing && length===nothing
        throw(ArgumentError("At least one of `length` or `step` must be specified"))
    end
end

if VERSION < v"1.1.0-DEV.506"
    function Base.range(start, stop; kwargs...)
        rangeargcheck(;kwargs...)
        range(start; stop=stop, kwargs...)
    end
end

# https://github.com/JuliaLang/julia/pull/30496
if VERSION < v"1.2.0-DEV.272"
    Base.@pure hasfield(::Type{T}, name::Symbol) where T =
        Base.fieldindex(T, name, false) > 0
    export hasfield
    hasproperty(x, s::Symbol) = s in propertynames(x)
    export hasproperty
end

# https://github.com/JuliaLang/julia/pull/29259
if VERSION < v"1.1.0-DEV.594"
    Base.merge(a::NamedTuple, b::NamedTuple, cs::NamedTuple...) = merge(merge(a, b), cs...)
    Base.merge(a::NamedTuple) = a
end

# https://github.com/JuliaLang/julia/pull/33129
if VERSION < v"1.4.0-DEV.142"
    export only

    Base.@propagate_inbounds function only(x)
        i = iterate(x)
        @boundscheck if i === nothing
            throw(ArgumentError("Collection is empty, must contain exactly 1 element"))
        end
        (ret, state) = i
        @boundscheck if iterate(x, state) !== nothing
            throw(ArgumentError("Collection has multiple elements, must contain exactly 1 element"))
        end
        return ret
    end

    # Collections of known size
    only(x::Ref) = x[]
    only(x::Number) = x
    only(x::Char) = x
    only(x::Tuple{Any}) = x[1]
    only(x::Tuple) = throw(
        ArgumentError("Tuple contains $(length(x)) elements, must contain exactly 1 element")
    )
    only(a::AbstractArray{<:Any, 0}) = @inbounds return a[]
    only(x::NamedTuple{<:Any, <:Tuple{Any}}) = first(x)
    only(x::NamedTuple) = throw(
        ArgumentError("NamedTuple contains $(length(x)) elements, must contain exactly 1 element")
    )
end

# https://github.com/JuliaLang/julia/pull/32628
if VERSION < v"1.3.0-alpha.8"
    Base.mod(i::Integer, r::Base.OneTo) = mod1(i, last(r))
    Base.mod(i::Integer, r::AbstractUnitRange{<:Integer}) = mod(i-first(r), length(r)) + first(r)
end

# https://github.com/JuliaLang/julia/pull/33568
if VERSION < v"1.4.0-DEV.329"
    Base.:∘(f, g, h...) = ∘(f ∘ g, h...)
end

include("deprecated.jl")

end # module Compat


text:
struct NLP <: JunctionLocationSolver
    optimizer          # for underlying NLP model
    timelimit::Float64 # seconds

    function NLP(solver; timelimit=Inf)
        new(solver, timelimit)
    end
end

"""
Direct implementation of nonconvex NLP model, using split-pipe cost.

Optionally force discrete diameter choice for every segment.
"""
function make_nlp(inst::Instance, topo::Topology, solver::NLP;
                  discdiam=false)
    nodes, nnodes = topo.nodes, length(topo.nodes)
    arcs, narcs = topo.arcs, length(topo.arcs)
    terms, nterms = inst.nodes, length(inst.nodes)
    termidx = termindex(nodes, terms)
    ndiams = length(inst.diameters)
    xmin, xmax = extrema(node.x for node in terms)
    ymin, ymax = extrema(node.y for node in terms)

    # pressure bounds for all nodes
    πlb_min = minimum([b.lb for b in inst.pressure])^2
    πub_max = maximum([b.ub for b in inst.pressure])^2
    πlb = fill(πlb_min, nnodes)
    πub = fill(πub_max, nnodes)
    πlb[termidx] = [b.lb^2 for b in inst.pressure]
    πub[termidx] = [b.ub^2 for b in inst.pressure]

    q = uniq_flow(inst, topo)

    c = [diam.cost for diam in inst.diameters]
    Dm5 = [diam.value^(-5) for diam in inst.diameters]
    C = inst.ploss_coeff .* q .* abs.(q)

    model = JuMP.Model(solver.optimizer)

    # node positions
    @variable(model, xmin ≤ x[1:nnodes] ≤ xmax)
    @variable(model, ymin ≤ y[1:nnodes] ≤ ymax)

    # length of arcs
    @variable(model, L[1:narcs] ≥ 0)

    # relative length of segments
    if discdiam
        @variable(model, l[1:narcs,1:ndiams], Bin)
    else
        @variable(model, l[1:narcs,1:ndiams] ≥ 0)
    end

    # squared pressure on nodes
    @variable(model, πlb[v] ≤ π[v=1:nnodes] ≤ πub[v])

    # objective proxy variable
    @variable(model, z ≥ 0)

    # fix terminal positions
    @constraint(model, fix[t=1:nterms], x[termidx[t]] == terms[t].x)
    @constraint(model, fiy[t=1:nterms], y[termidx[t]] == terms[t].y)

    # couple positions and length
    @constraint(model, xylength[a=1:narcs],
                L[a]^2 == (x[arcs[a].head] - x[arcs[a].tail])^2 + (y[arcs[a].head] - y[arcs[a].tail])^2)

    # convex combination of diameters
    @constraint(model, convexcomb[a=1:narcs], sum(l[a,i] for i=1:ndiams) == 1.0)

    # pressure loss
    @constraint(model, ploss[a=1:narcs], π[arcs[a].tail] - π[arcs[a].head] ==
                C[a] * L[a] * sum(Dm5[i]*l[a,i] for i=1:ndiams))

    # minimize total construction cost
    @constraint(model, objfun,
                z ≥ sum(c[i] * L[a] * l[a,i] for a=1:narcs for i=1:ndiams))
    @objective(model, Min, z)

    model, x, y, L, l, π
end

function PipeLayout.optimize(inst::Instance, topo::Topology, solver::NLP)
    model, x, y, L, l, π = make_nlp(inst, topo, solver)
    JuMP.optimize!(model)
    status = JuMP.termination_status(model)
    if status in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED]
        return Result(status, Solution([], zeros(0,0), []), Inf)
    end
    objval = JuMP.objective_value(model)
    nodes = map(Node, zip(JuMP.value.(x), JuMP.value.(y)))
    sol = Solution(nodes, JuMP.value.(l), JuMP.value.(π))
    Result(status, sol, objval)
end


text:
using SymPy
using Test
using SpecialFunctions

@testset "Math" begin
    ρ = symbols("rho", positive=true)
    ϕ = symbols("phi", real=true)
    x = symbols("x")

    @test simplify(hypot(ρ*cos(ϕ), ρ * sin(ϕ))) == ρ
    @test simplify(hypot(ρ*cos(ϕ), 3)) == sqrt(ρ^2*cos(ϕ)^2 + 9)

    @test atan(Sym(1), 1) == PI/4
    @test atan(Sym(1), -1) == 3PI/4

    @test factorial(Sym(0)) == 1
    @test factorial(Sym(7)) == 5040

    @test sympy.factorial2(Sym(5)) == 15
    @test sympy.factorial2(Sym(-5)) == Sym(1)/3

    @test erf(Sym(0)) == 0
    @test erf(Sym(oo)) == 1
    @test diff(erf(x), x) == 2*exp(-x^2)/sqrt(PI)

    @test sinc(Sym(0)) == 1
    # test consistency with Julia's sinc
    @test sinc(Sym(1)) == 0
    @test N(sinc(Sym(0.2))) ≈ sinc(0.2)

    @test flipsign(Sym(3), 2.0) == 3
    @test flipsign(Sym(3), 0.0) == 3
    @test flipsign(Sym(3), -0.0) == -3
    @test flipsign(Sym(3), -2.0) == -3

    @test eps(Sym) == 0
    #@test rewrite(sinc(x), "jn") == jn(0, PI * x)

    @test abs2(x) == conj(x)*x
    @test abs(x) == sqrt(conj(x)*x)
end