text:
export algebra

add_assert_scope(:AlgAssOrd)
add_verbose_scope(:AlgAssOrd)

elem_type(::AlgAssAbsOrd{S, T}) where {S, T} = AlgAssAbsOrdElem{S, T}

elem_type(::Type{AlgAssAbsOrd{S, T}}) where {S, T} = AlgAssAbsOrdElem{S, T}

ideal_type(::AlgAssAbsOrd{S, T}) where {S, T} = AlgAssAbsOrdIdl{S, T}

ideal_type(::Type{AlgAssAbsOrd{S, T}}) where {S, T} = AlgAssAbsOrdIdl{S, T}

# There is no dedicated type for fractional ideals
fractional_ideal_type(::AlgAssAbsOrd{S, T}) where {S, T} = AlgAssAbsOrdIdl{S, T}

fractional_ideal_type(::Type{AlgAssAbsOrd{S, T}}) where {S, T} = AlgAssAbsOrdIdl{S, T}

@doc Markdown.doc"""
    algebra(O::AlgAssAbsOrd) -> AbsAlgAss

Returns the algebra which contains $O$.
"""
algebra(O::AlgAssAbsOrd) = O.algebra

_algebra(O::AlgAssAbsOrd) = algebra(O)

base_ring(O::AlgAssAbsOrd) = FlintZZ

@doc Markdown.doc"""
    iscommutative(O::AlgAssAbsOrd) -> Bool

Returns `true` if $O$ is a commutative ring and `false` otherwise.
"""
iscommutative(O::AlgAssAbsOrd) = iscommutative(algebra(O))

ismaximal_known(O::AlgAssAbsOrd) = O.ismaximal != 0

@doc Markdown.doc"""
    ismaximal(O::AlgAssAbsOrd) -> Bool

Returns `true` if $O$ is a maximal order and `false` otherwise.
"""
function ismaximal(O::AlgAssAbsOrd)
  if O.ismaximal == 1
    return true
  end
  if O.ismaximal == 2
    return false
  end

  A = algebra(O)
  d = discriminant(O)
  if isdefined(A, :maximal_order)
    if d == discriminant(maximal_order(A))
      O.ismaximal = 1
      return true
    else
      O.ismaximal = 2
      return false
    end
  end

  if typeof(A) <: AlgGrp
    fac = factor(degree(O))
  else
    fac = factor(abs(d))
  end

  for (p, j) in fac
    if j == 1
      continue
    end
    d2 = discriminant(pmaximal_overorder(O, Int(p)))
    if d != d2
      O.ismaximal = 2
      return false
    end
  end
  O.ismaximal = 1
  return true
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    Order(A::AbsAlgAss{fmpq}, B::Vector{<: AbsAlgAssElem{fmpq}}; check::Bool = true,
          isbasis::Bool = false, cached::Bool = true)
      -> AlgAssAbsOrd

Returns the order of $A$ generated by $B$. If `check` is set, it is checked
whether $B$ defines an order. If `isbasis` is set, then the elements are
assumed to form a $\mathbb Z$-basis.
"""
function Order(A::S, B::Vector{T}; check::Bool = true, isbasis::Bool = false, cached::Bool = true) where {S <: AbsAlgAss{fmpq}, T <: AbsAlgAssElem{fmpq}}
  if isbasis
    if check
      b, bmat, bmat_inv, _ = defines_order(A, B)
      if !b
        error("The elements do not define an order")
      else
        return AlgAssAbsOrd{S, elem_type(S)}(A, bmat, bmat_inv, deepcopy(B), cached)
      end
    else
      return AlgAssAbsOrd{S, elem_type(S)}(A, deepcopy(B), cached)
    end
  else
    return _order(A, B; cached = cached, check = check)
  end
end

@doc Markdown.doc"""
    Order(A::AbsAlgAss{fmpq}, M::FakeFmpqMat; check::Bool = true,
          cached::Bool = true)
      -> AlgAssAbsOrd

Returns the order of $A$ with basis matrix $M$. If `check` is set, it is checked
whether $M$ defines an order.
"""
function Order(A::S, M::FakeFmpqMat; check::Bool = true, cached::Bool = true) where {S <: AbsAlgAss{fmpq}}
  if check
    b, Minv, v = defines_order(A, M)
    if !b
      error("The basis matrix does not define an order")
    else
      return AlgAssAbsOrd{S, elem_type(S)}(A, deepcopy(M), Minv, v, cached)
    end
  else
    return AlgAssAbsOrd{S, elem_type(S)}(A, deepcopy(M), cached)
  end
end

function _order(A::S, gens::Vector{T}; cached::Bool = true, check::Bool = true) where {S <: AbsAlgAss, T <: AbsAlgAssElem}
  B_A = basis(A)

  if one(A) in gens
    cur = gens
  else
    cur = append!([one(A)], gens)
  end
  Bmat = basis_matrix(cur, FakeFmpqMat)
  while true
    k = length(cur)
    prods = Vector{elem_type(A)}(undef, k^2)
    for i = 1:k
      ik = (i - 1)*k
      for j = 1:k
        prods[ik + j] = cur[i]*cur[j]
      end
    end
    Ml = hnf(basis_matrix(prods, FakeFmpqMat))
    r = findfirst(i -> !iszero_row(Ml.num, i), 1:k^2)
    nBmat = sub(Ml, r:nrows(Ml), 1:ncols(Ml))
    if nrows(nBmat) == nrows(Bmat) && Bmat == nBmat
      break
    end
    Bmat = nBmat
  end
  if nrows(Bmat) != dim(A)
    error("Elements do not generate an order")
  end

  return Order(A, Bmat, cached = cached, check = check)
end

function _equation_order(A::AbsAlgAss{fmpq})
  @assert iscommutative(A)
  a = primitive_element_via_number_fields(A)
  b = Vector{elem_type(A)}(undef, dim(A))
  b[1] = one(A)
  for i = 2:dim(A)
    b[i] = b[i - 1]*a
  end
  return Order(A, b)
end

################################################################################
#
#  Index
#
################################################################################

function index(O::AlgAssAbsOrd)
  B = basis_mat_inv(O, copy = false)
  n = det(B)
  @assert isinteger(n)
  return FlintZZ(n)
end

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function _assure_has_basis(O::AlgAssAbsOrd)
  if !isdefined(O, :basis)
    B = basis(algebra(O))
    M = basis_matrix(O, copy = false)
    v = Vector{elem_type(O)}(undef, degree(O))
    for i in 1:degree(O)
      w = sum(M.num[i, j]//M.den * B[j] for j in 1:degree(O))
      v[i] = O(w)
    end
    O.basis = v
  end
  return nothing
end

function assure_basis_mat_inv(O::AlgAssAbsOrd)
  if !isdefined(O, :basis_mat_inv)
    O.basis_mat_inv = inv(basis_matrix(O, copy = false))
  end
  return nothing
end

function assure_basis_alg(O::AlgAssAbsOrd)
  if isdefined(O, :basis_alg)
    return nothing
  end

  M = basis_matrix(O, copy = false)
  A = algebra(O)
  O.basis_alg = Vector{elem_type(A)}(undef, dim(A))
  for i = 1:dim(A)
    O.basis_alg[i] = elem_from_mat_row(A, M.num, i, M.den)
  end
  return nothing
end

################################################################################
#
#  Basis
#
################################################################################

@doc Markdown.doc"""
    basis(O::AlgAssAbsOrd; copy::Bool = true) -> Vector{AlgAssAbsOrdElem}

Returns a $\mathbb Z$-basis of $O$.
"""
function basis(O::AlgAssAbsOrd; copy::Bool = true)
  _assure_has_basis(O)
  if copy
    return deepcopy(O.basis)::Vector{elem_type(O)}
  else
    return O.basis::Vector{elem_type(O)}
  end
end

absolute_basis(O::AlgAssAbsOrd) = basis(O)

function basis_alg(O::AlgAssAbsOrd{S, T}; copy::Bool = true) where {S, T}
  assure_basis_alg(O)
  if copy
    return deepcopy(O.basis_alg)::Vector{T}
  else
    return O.basis_alg::Vector{T}
  end
end

################################################################################
#
#  (Inverse) basis matrix
#
################################################################################

@doc Markdown.doc"""
    basis_matrix(O::AlgAssAbsOrd; copy::Bool = true) -> FakeFmpqMat

Returns the basis matrix of $O$.
"""
function basis_matrix(x::AlgAssAbsOrd; copy::Bool = true)
  if copy
    return deepcopy(x.basis_matrix)
  else
    return x.basis_matrix
  end
end

@doc Markdown.doc"""
    basis_mat_inv(O::AlgAssAbsOrd; copy::Bool = true) -> FakeFmpqMat

Returns the inverse of the basis matrix of $O$.
"""
function basis_mat_inv(O::AlgAssAbsOrd; copy::Bool = true)
  assure_basis_mat_inv(O)
  if copy
    return deepcopy(O.basis_mat_inv)::FakeFmpqMat
  else
    return O.basis_mat_inv::FakeFmpqMat
  end
end

################################################################################
#
#  Degree
#
################################################################################

@doc Markdown.doc"""
    degree(O::AlgAssAbsOrd) -> Int

Returns the dimension of the algebra containing $O$.
"""
function degree(O::AlgAssAbsOrd)
  return dim(algebra(O))
end

################################################################################
#
#  Inclusion of algebra elements
#
################################################################################

function _check_elem_in_order(a::T, O::AlgAssAbsOrd{S, T}, short::Type{Val{U}} = Val{false}) where {S, T, U}
  t = zero_matrix(FlintQQ, 1, degree(O))
  elem_to_mat_row!(t, 1, a)
  t = FakeFmpqMat(t)
  t = t*basis_mat_inv(O, copy = false)
  if short == Val{true}
    return isone(t.den)
  elseif short == Val{false}
    if !isone(t.den)
      return false, Vector{fmpz}()
    else
      v = Vector{fmpz}(undef, degree(O))
      for i = 1:degree(O)
        v[i] = deepcopy(t.num[1, i])
      end
      return true, v
    end
  end
end

@doc Markdown.doc"""
    in(x::AbsAlgAssElem, O::AlgAssAbsOrd) -> Bool

Returns `true` if the algebra element $x$ is in $O$ and `false` otherwise.
"""
function in(x::T, O::AlgAssAbsOrd{S, T}) where {S, T}
  return _check_elem_in_order(x, O, Val{true})
end

################################################################################
#
#  Denominator in an order
#
################################################################################

@doc Markdown.doc"""
    denominator(a::AbsAlgAssElem, O::AlgAssAbsOrd) -> fmpz

Returns $d\in \mathbb Z$ such that $d \cdot a \in O$.
"""
function denominator(a::AbsAlgAssElem, O::AlgAssAbsOrd)
  t = zero_matrix(FlintQQ, 1, degree(O))
  elem_to_mat_row!(t, 1, a)
  t = FakeFmpqMat(t)
  t = mul!(t, t, basis_mat_inv(O, copy = false))
  return t.den
end

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

RandomExtensions.maketype(O::AlgAssAbsOrd, R::UnitRange) = elem_type(O)

function rand(rng::AbstractRNG,
              sp::SamplerTrivial{<:Make2{<:RingElem,<:AlgAssAbsOrd,<:UnitRange}})
  O, R = sp[][1:2]
  O(map(fmpz, rand(rng, R, degree(O))))
end

RandomExtensions.make(O::AlgAssAbsOrd, n::Union{Integer, fmpz}) =
  make(O, Integer(-n):Integer(n))

@doc Markdown.doc"""
    rand(O::AlgAssAbsOrd, R::UnitRange) -> AlgAssAbsOrdElem

Returns a random element of $O$ whose coefficients lie in $R$.
"""
rand(O::AlgAssAbsOrd, R::UnitRange) = rand(Random.GLOBAL_RNG, O, R)

@doc Markdown.doc"""
    rand(O::AlgAssAbsOrd, n::Union{Integer, fmpz}) -> AlgAssAbsOrdElem

Returns a random element of $O$ whose coefficients are bounded by $n$.
"""
rand(O::AlgAssAbsOrd, n::fmpz) = rand(Random.GLOBAL_RNG, O, n)
rand(O::AlgAssAbsOrd, n::Integer) = rand(Random.GLOBAL_RNG, O, n)
# these two methods can't be merged with a Union, because of ambiguities

rand(rng::AbstractRNG, O::AlgAssAbsOrd, n::Union{fmpz,UnitRange}) = rand(rng, make(O, n))
rand(rng::AbstractRNG, O::AlgAssAbsOrd, n::Integer) = rand(rng, make(O, n))

################################################################################
#
#  Basis matrices from generators
#
################################################################################

function basis_matrix(A::Array{S, 1}, ::Type{FakeFmpqMat}) where {S <: AbsAlgAssElem{fmpq}}
  @assert length(A) > 0
  n = length(A)
  d = dim(parent(A[1]))

  M = zero_matrix(FlintZZ, n, d)

  dens = [lcm([denominator(coefficients(A[i], copy = false)[j]) for j=1:d]) for i=1:n]
  deno = lcm(dens)

  for i in 1:n
    for j in 1:d
      temp_den = divexact(deno, denominator(coefficients(A[i], copy = false)[j]))
      M[i, j] = numerator(coefficients(A[i], copy = false)[j]) * temp_den
    end
  end
  return FakeFmpqMat(M, deno)
end

function basis_matrix(A::Vector{ <: AbsAlgAssElem{T} }) where T
  @assert length(A) > 0
  n = length(A)
  d = dim(parent(A[1]))
  K = base_ring(parent(A[1]))

  M = zero_matrix(K, n, d)

  for i = 1:n
    elem_to_mat_row!(M, i, A[i])
    #for j = 1:d
    #  M[i, j] = deepcopy(coefficients(A[i], copy = false)[j])
    #end
  end
  return M
end

function basis_matrix(A::Array{AlgAssAbsOrdElem{S, T}, 1}) where S where T
  @assert length(A) > 0
  n = length(A)
  d = degree(parent(A[1]))
  M = zero_matrix(FlintZZ, n, d)

  for i in 1:n
    el = coordinates(A[i])
    for j in 1:d
      M[i, j] = el[j]
    end
  end
  return M
end

################################################################################
#
#  Sum of orders
#
################################################################################

# Be careful!
# To be used only in the case of the construction of a maximal order!
function +(a::AlgAssAbsOrd, b::AlgAssAbsOrd)
  aB = basis_matrix(a, copy = false)
  bB = basis_matrix(b, copy = false)
  d = degree(a)
  c = sub(_hnf(vcat(bB.den*aB.num, aB.den*bB.num), :lowerleft), d + 1:2*d, 1:d)
  return Order(algebra(a), FakeFmpqMat(c, aB.den*bB.den))
end

################################################################################
#
#  Print
#
################################################################################

function show(io::IO, O::AlgAssAbsOrd)
  compact = get(io, :compact, false)
  if compact
    print(io, "Order of ")
    show(IOContext(io, :compact => true), algebra(O))
  else
    print(io, "Order of ")
    print(io, algebra(O))
    println(io, " with basis matrix ")
    print(io, basis_matrix(O))
  end
end

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

@doc Markdown.doc"""
    ==(S::AlgAssAbsOrd, T::AlgAssAbsOrd) -> Bool

Returns `true` if $S$ and $T$ are equal and `false` otherwise.
"""
function ==(S::AlgAssAbsOrd, T::AlgAssAbsOrd)
  return basis_matrix(S, copy = false) == basis_matrix(T, copy = false)
end

################################################################################
#
#  Discriminant and Reduced Trace Matrix
#
################################################################################

@doc Markdown.doc"""
    trred_matrix(O::AlgAssAbsOrd) -> fmpz_mat

Returns the reduced trace matrix $M$ of $O$, i. e. `M[i, j] = trred(b[i]*b[j])`,
where $b$ is a basis of $O$.
"""
function trred_matrix(O::AlgAssAbsOrd)
  if isdefined(O, :trred_matrix)
    return O.trred_matrix
  end
  A=algebra(O)
  x=O.basis_alg
  m=length(x)
  M=zero_matrix(FlintZZ, m, m)
  a=A()
  for i=1:m
    a = mul!(a, x[i], x[i])
    M[i,i] = FlintZZ(trred(a))
  end
  for i = 1:m
    for j = i+1:m
      a = mul!(a, x[i], x[j])
      b = FlintZZ(trred(a))
      M[i,j] = b
      M[j,i] = b
    end
  end
  O.trred_matrix = M
  return M
end

@doc Markdown.doc"""
    discriminant(O::AlgAssAbsOrd) -> fmpz

Returns the discriminant of $O$.
"""
function discriminant(O::AlgAssAbsOrd)
  if isdefined(O, :disc)
    return O.disc
  end
  M = trred_matrix(O)
  O.disc = det(M)
  return O.disc
end

################################################################################
#
#  Schur Index at Infinity
#
################################################################################

#Steel Nebe paper
@doc Markdown.doc"""
    schur_index_at_real_plc(O::AlgAssAbsOrd) -> Int

Given an order $O$, this function returns the schur index
of the algebra over the field of real numbers.
"""
function schur_index_at_real_plc(O::AlgAssAbsOrd)

  x=trace_signature(O)
  n=root(degree(O),2)
  if x[1] == divexact(n*(n+1),2)
    return 1
  else
    return 2
  end
end

function trace_signature(O::AlgAssAbsOrd)

  @vtime :AlgAssOrd 1 M = trred_matrix(O)
  Zx, x = PolynomialRing(FlintZZ, "x", cached = false)
  Qy, y = PolynomialRing(FlintQQ, "y", cached = false)
  @vtime :AlgAssOrd 1 f = charpoly(Zx, M)
  @vtime :AlgAssOrd 1 fac = factor_squarefree(Qy(f))
  npos = 0
  for (t,e) in fac
    @vtime :AlgAssOrd a = number_positive_roots(Zx(t))
    npos += a*e
  end
  return (npos, degree(f) - npos)
end

################################################################################
#
#  Schur Index at p
#
################################################################################

@doc Markdown.doc"""
    schur_index_at_p(O::AlgAssAbsOrd, p::fmpz)

Given a maximal order $O$ and a prime $p$, this function returns the schur index
of the completion of the algebra at $p$.
"""
function schur_index_at_p(O::AlgAssAbsOrd, p::fmpz)
  @assert O.ismaximal==1
  d = discriminant(O)
  v = valuation(d,p)
  if v == 0
    return 1
  end
  n = root(degree(O),2)
  t = n - divexact(v,n)
  return divexact(n,t)
end

################################################################################
#
#  p-maximal overorder
#
################################################################################

function pmaximal_overorder(O::AlgAssAbsOrd{S, T}, p::Union{fmpz, Int}) where S where T
  d = discriminant(O)
  if rem(d, p^2) != 0
    return O
  end

  if p > degree(O)
    @vtime :AlgAssOrd 1 O1 = pmaximal_overorder_tr(O,p)::AlgAssAbsOrd{S, T}
    return O1
  else
    @vtime :AlgAssOrd 1 O1 = pmaximal_overorder_meataxe(O,p)::AlgAssAbsOrd{S, T}
    return O1
  end
end

function pmaximal_overorder_meataxe(O::AlgAssAbsOrd{S, T}, p::Union{fmpz, Int}) where {S, T}

  extend = false
  d = discriminant(O)
  while true
    dd = fmpz(1)
    @vtime :AlgAssOrd 1 max_id =_maximal_ideals(O, p*O, p, strict_containment = true)
    for m in max_id
      @vtime :AlgAssOrd 1 OO = _ring_of_multipliers_integral_ideal(m, fmpz(p))
      dd = discriminant(OO)
      if d != dd
        extend = true
        O = OO
        d = dd
        break
      end
    end

    if extend
      if rem(d, p^2) != 0
        break
      end
      extend = false
      continue
    else
      break
    end

  end
  return O
end

function pmaximal_overorder_tr(O::AlgAssAbsOrd, p::Int)
  #First, the head order by computing the pradical and its ring of multipliers
  d = discriminant(O)
  @vtime :AlgAssOrd 1 I = pradical(O, p)
  @vtime :AlgAssOrd 1 OO = _ring_of_multipliers_integral_ideal(I, fmpz(p))
  dd = discriminant(OO)
  if rem(dd, p^2) != 0
    return OO
  end
  while dd!= d
    d = dd
    O = OO
    @vtime :AlgAssOrd 1 I = pradical(O,p)
    @vtime :AlgAssOrd 1 OO = _ring_of_multipliers_integral_ideal(I, fmpz(p))
    dd = discriminant(OO)
    if rem(dd, p^2) != 0
      return OO
    end
  end
  #Now, we have to check the maximal ideals.

  extend = false
  @vtime :AlgAssOrd 1 max_id = _maximal_ideals(O, I, p, strict_containment = true)
  for m in max_id
    @vtime :AlgAssOrd 1 OO = _ring_of_multipliers_integral_ideal(m, fmpz(p))
    dd = discriminant(OO)
    if d != dd
      extend = true
      O = OO
      d = dd
      break
    end
  end
  if extend
    if rem(dd, p^2) != 0
      return OO
    end
    extend = false
  else
    return OO
  end
  while true
    dd = fmpz(1)
    @vtime :AlgAssOrd 1 max_id = _maximal_ideals(O, p*O, p, strict_containment = true)
    for m in max_id
      OO = _ring_of_multipliers_integral_ideal(m, fmpz(p))
      dd = discriminant(OO)
      if d != dd
        extend = true
        O = OO
        d = dd
        break
      end
    end

    if extend
      if rem(dd, p^2) != 0
        break
      end
      extend = false
      continue
    else
      break
    end

  end
  return O
end

################################################################################
#
#  Maximal Order
#
################################################################################

@doc Markdown.doc"""
    MaximalOrder(O::AlgAssAbsOrd)

Given an order $O$, this function returns a maximal order containing $O$.
"""
function MaximalOrder(O::AlgAssAbsOrd{S, T}) where S where T
  A = algebra(O)

  if isdefined(A, :maximal_order)
    # Check whether O \subseteq OO
    OO = A.maximal_order::order_type(A)
    d = denominator(basis_matrix(O, copy = false)*basis_mat_inv(OO, copy = false))
    if isone(d)
      return OO
    end
  end

  d = discriminant(O)
  @vtime :NfOrd fac = factor(abs(d))

  OO = O
  for (p, j) in fac
    if mod(d, p^2) != 0
      continue
    end
    OO += pmaximal_overorder(O, Int(p))
  end
  OO.ismaximal = 1

  if !isdefined(A, :maximal_order)
    A.maximal_order = OO
  end
  return OO
end

function MaximalOrder(O::AlgAssAbsOrd{S, T}) where { S <: AlgGrp, T <: AlgGrpElem }
  A = algebra(O)

  if isdefined(A, :maximal_order)
    # Check whether O \subseteq OO
    OO = A.maximal_order::order_type(A)
    d = denominator(basis_matrix(O, copy = false)*basis_mat_inv(OO, copy = false))
    if isone(d)
      return OO
    end
  end
  d = discriminant(O)
  fac = factor(degree(O)) # the order of the group

  OO = O
  for (p, j) in fac
    if mod(d, p^2) != 0
      continue
    end
    OO += pmaximal_overorder(O, Int(p))
  end
  OO.ismaximal = 1

  if !isdefined(A, :maximal_order)
    A.maximal_order = OO
  end
  return OO
end

function _denominator_of_mult_table(A::AbsAlgAss{fmpq})
  l = denominator(multiplication_table(A, copy = false)[1, 1, 1])
  for i = 1:dim(A)
    for j = 1:dim(A)
      for k = 1:dim(A)
        l = lcm(l, denominator(multiplication_table(A, copy = false)[i, j, k]))
      end
    end
  end
  return l
end

_denominator_of_mult_table(A::AlgGrp{fmpq}) = fmpz(1)

@doc Markdown.doc"""
    any_order(A::AbsAlgAss{fmpq}) -> AlgAssAbsOrd

Returns any order of $A$.
"""
function any_order(A::AbsAlgAss{fmpq})
  d = _denominator_of_mult_table(A)

  M = vcat(zero_matrix(FlintQQ, 1, dim(A)), d*identity_matrix(FlintQQ, dim(A)))
  oneA = one(A)
  for i = 1:dim(A)
    M[1, i] = deepcopy(coefficients(oneA, copy = false)[i])
  end
  M = FakeFmpqMat(M)
  M = hnf!(M, :lowerleft)
  O = Order(A, sub(M, 2:dim(A) + 1, 1:dim(A)))
  return O
end

@doc Markdown.doc"""
    MaximalOrder(A::AbsAlgAss{fmpq}) -> AlgAssAbsOrd

Returns a maximal order of $A$.
"""
function MaximalOrder(A::AbsAlgAss{S}) where S
  if isdefined(A, :maximal_order)
    return A.maximal_order::AlgAssAbsOrd{typeof(A), elem_type(A)}
  end

  O = any_order(A)
  OO = MaximalOrder(O)
  A.maximal_order = OO
  return OO
end

function maximal_order_via_decomposition(A::AbsAlgAss{fmpq})
  if isdefined(A, :maximal_order)
    return A.maximal_order::AlgAssAbsOrd{AlgAss{fmpq},AlgAssElem{fmpq,AlgAss{fmpq}}}
  end
  fields_and_maps = as_number_fields(A)
  M = zero_matrix(FlintQQ, dim(A), dim(A))
  row = 1
  for i = 1:length(fields_and_maps)
    K = fields_and_maps[i][1]
    AtoK = fields_and_maps[i][2]
    O = maximal_order(K)
    for b in basis(O)
      a = AtoK\K(b)
      elem_to_mat_row!(M, row, a)
      row += 1
    end
  end
  FakeM = FakeFmpqMat(M)
  FakeM = hnf!(FakeM, :lowerleft)
  OO = Order(A, FakeM)
  OO.ismaximal = 1
  A.maximal_order = OO
  return OO
end

# Requires that O is maximal and A = QQ^(n\times n).
# Computes a maximal order of type
#  (O ... O a^(-1))
#  (:     :   :   )
#  (O ... O a^(-1))
#  (a ... a   O   )
# for an ideal a of O.
# See Bley, Johnston "Computing generators of free modules over orders in group
# algebras", Prop. 5.1.
function _simple_maximal_order(O::AlgAssAbsOrd{S1, S2}, with_trafo::Type{Val{T}} = Val{false}) where { S1 <: AlgMat, S2, T }
  A = algebra(O)

  if !(A isa AlgMat)
    throw(ArgumentError("Order must be an order in a matrix algebra"))
  end

  n = degree(A)

  # Build a matrix with the first rows of basis elements of O
  M = zero_matrix(FlintQQ, dim(A), n)
  for i = 1:dim(A)
    for j = 1:n
      M[i, j] = deepcopy(matrix(elem_in_algebra(basis(O, copy = false)[i], copy = false), copy = false)[1, j])
    end
  end
  M = FakeFmpqMat(M)
  M = hnf!(M, :upperright)
  M = fmpq_mat(sub(M, 1:n, 1:n))

  # Compute M * O * M^-1
  iM = inv(M)
  bb = Vector{elem_type(A)}()
  for i = 1:degree(O)
    push!(bb, M*elem_in_algebra(basis(O, copy = false)[i], copy = false)*iM)
  end

  simpleOrder = Order(A, bb)
  simpleOrder.isnice = true

  @assert basis_matrix(simpleOrder) == FakeFmpqMat(identity_matrix(FlintQQ, n^2))

  if with_trafo == Val{true}
    return simpleOrder, A(M)
  else
    return simpleOrder
  end
end

function issimple(O::AlgAssAbsOrd)
  return O.issimple
end

@doc Markdown.doc"""
    nice_order(O::AlgAssAbsOrd) -> AlgAssAbsOrd, AlgElem

Given a maximal order `O` in a full matrix algebra over the rationals, return a
nice maximal order `R` and element `a` such that `a O a^-1 = R`.
"""
function nice_order(O::AlgAssAbsOrd{S, T}) where {S, T}
  if isdefined(O, :nice_order)
    return O.nice_order::Tuple{typeof(O), T}
  else
    sO, A = _simple_maximal_order(O, Val{true})
    O.nice_order = sO, A
    return sO, A
  end
end

################################################################################
#
#  Conductor
#
################################################################################

@doc Markdown.doc"""
    conductor(R::AlgAssAbsOrd, S::AlgAssAbsOrd, action::Symbol) -> AlgAssAbsOrdIdl

Returns the ideal $\{ x \in R \mid xS \subseteq R \}$ if `action == :right` and
$\{ x \in R \mid Sx \subseteq R \}$ if `action == :left`.
It is assumed that $R \subseteq S$.
"""
function conductor(R::AlgAssAbsOrd, S::AlgAssAbsOrd, action::Symbol = :left)
  n = degree(R)
  t = basis_matrix(R, copy = false)*basis_mat_inv(S, copy = false)
  @assert isone(t.den)
  basis_mat_R_in_S_inv_num, d = pseudo_inv(t.num)
  M = zero_matrix(FlintZZ, n^2, n)
  B = basis(S, copy = false)
  for k in 1:n
    a = B[k]
    N = representation_matrix(a, action)*basis_mat_R_in_S_inv_num
    for i in 1:n
      for j in 1:n
        M[(k - 1)*n + i, j] = N[j, i]
      end
    end
  end
  H = sub(hnf(M), 1:n, 1:n)
  Hinv = inv(FakeFmpqMat(transpose(H)))
  Hinv = Hinv*basis_mat_R_in_S_inv_num*basis_matrix(R, copy = false)
  if action == :left
    return ideal(algebra(R), R, Hinv, :right)
  else
    return ideal(algebra(R), R, Hinv, :left)
  end
end

################################################################################
#
#  Units of quotients
#
################################################################################

# Computes a generating system of U in O, where U is a set of representatives of
# the image of the projection map \pi:O^\times -> (O/g*O)^\times.
# Assumes that O is a maximal order in Mat_{n\times n}(QQ).
# See Bley, Johnson: "Computing generators of free modules over orders in
# group algebras", section 6.
function enum_units(O::AlgAssAbsOrd{S, T}, g::fmpz) where { S <: AlgMat, T }
  A = algebra(O)
  @assert degree(A)^2 == dim(A)

  n = degree(A)

  L = _simple_maximal_order(O)
  a = basis_matrix(L, copy = false)[dim(A) - 1, dim(A) - 1]
  ai = basis_matrix(L, copy = false)[n, n]

  result = Vector{elem_type(L)}()
  n1 = n - 1
  # n \nmid i, j or n \mid i, j
  for i = 1:n1
    for j = 1:n1
      if j == i
        continue
      end
      E = identity_matrix(FlintQQ, n)
      E[i, j] = deepcopy(g)
      push!(result, L(A(E)))
    end
  end

  # n \nmid i and n \mid j
  for i = 1:n1
    E = identity_matrix(FlintQQ, n)
    E[i, n] = deepcopy(a)
    push!(result, L(A(E)))
  end

  # n \mid i and n \nmid j
  for j = 1:n1
    E = identity_matrix(FlintQQ, n)
    E[n, j] = deepcopy(ai)
    push!(result, L(A(E)))
  end

  E = identity_matrix(FlintQQ, n)
  E[1, 1] = fmpz(-1)
  push!(result, L(A(E)))
  return result
end

################################################################################
#
#  Trace dual ideal
#
################################################################################

function trace_dual(R::AlgAssAbsOrd)
  t = inv(FakeFmpqMat(trred_matrix(R)))*basis_matrix(R, copy = false)
  return ideal(algebra(R), R, t)
end

################################################################################
#
#  "All" maximal orders
#
################################################################################

# Only works for algebras fulfilling the Eichler condition.
# This is trivial for algebras over QQ, as there is always just one equivalence
# class with respect to conjugation.
representatives_of_maximal_orders(A::AlgAss{fmpq}) = representatives_of_maximal_orders(maximal_order(A))

function representatives_of_maximal_orders(O::AlgAssAbsOrd)
  A = algebra(O)
  @assert issimple(A)
  @assert iseichler(A)
  @assert ismaximal(O)
  return typeof(O)[ O ]
end


text:
module CFTutils

using PyPlot
using PyCall
using PlotTools

immutable CFTField
    w::Float64 #conformal weight
    d::Vector{Int} #Virasoro characters
end

immutable CFTOp
    fh::CFTField
    fa::CFTField
    lbl::AbstractString
end

sdim(oph::CFTField, Ih::Int, opa::CFTField, Ia::Int) = (oph.d[Ih+1] * opa.d[Ia+1], oph.w + Ih + opa.w + Ia, oph.w + Ih - opa.w - Ia)

sdim(op::CFTOp, Ih::Int, Ia::Int) = sdim(op.fh, Ih, op.fa, Ia)

function sop_count(oph::CFTField, opa::CFTField, maxsd::Float64, maxspin::Float64=Inf)
    Imax = floor(Int, maxsd - (oph.w + opa.w))
    numops = 0
    for I in 0:Imax
        for Ih in 0:I
            d, sd, sp = sdim(oph, Ih, opa, I-Ih)
            if abs(sp) <= maxspin
                numops += d
            end
        end
    end
    numops
end

sop_count(op::CFTOp, maxsd::Float64=15.0, maxspin::Float64=Inf) = sop_count(op.fh, op.fa, maxsd, maxspin)

function sops(oph::CFTField, opa::CFTField, maxsd::Float64=15.0)
    Imax = floor(Int, maxsd - (oph.w + opa.w))

    sds = Float64[]
    sps = Float64[]
    ds = Int[]
    for I in 0:Imax
        for Ih in 0:I
            d, sd, sp = sdim(oph, Ih, opa, I-Ih)
            push!(ds, d)
            push!(sds, sd)
            push!(sps, sp)
        end
    end
    sds, sps, ds
end

sops(op::CFTOp, maxsd::Float64=15.0) = sops(op.fh, op.fa, maxsd)

function sops_by_s(op::CFTOp, maxsd::Float64=15.0; maxspin::Float64=Inf)
    sds, sps, ds = sops(op, maxsd)

    ss = unique(round(sps,6))
    ss = ss[abs(ss) .< maxspin + sqrt(eps())]

    haspri = Bool[]
    ress = Vector{Float64}[]
    for s in ss
        flt = abs(sps - s) .< 1e-6

        push!(haspri, flt[1])

        sds_s = sds[flt]
        sps_s = sps[flt]
        ds_s = ds[flt]

        res = Float64[]
        for j in 1:length(sds_s)
            for k in 1:ds_s[j]
                append!(res, sds_s[j])
            end
        end
        push!(ress, res)
    end
    ss, ress, haspri
end

function plot_sops_v3(ops::Vector{CFTOp}, cs::Vector, maxsd::Float64=15.0; minsd::Float64=-Inf, maxspin::Float64=Inf, xsep=0.12, ysep=0.12, szp=40, sz=20, shift_pris::Bool=true)
    sops = [sops_by_s(op, maxsd, maxspin=maxspin) for op in ops]

    ss = unique(round(vcat([sop[1] for sop in sops]...),6))
    sds = Vector{Float64}[]
    cols = Vector[]
    ispri = Vector{Bool}[]
    for s in ss
        sds_s = Float64[]
        cols_s = []
        ispri_s = []
        for (j, (ss_op, sd_op, haspri_op)) in enumerate(sops)
            ind = findfirst(s_op->s_op == s, ss_op)
            if ind > 0
                sd_op = sd_op[ind]
                if length(sd_op) > 0
                    append!(sds_s, sd_op)
                    append!(cols_s, [cs[j] for k in 1:length(sd_op)])
                    is_pri_op = falses(length(sd_op))
                    is_pri_op[1] = haspri_op[ind]
                    append!(ispri_s, is_pri_op)
                end
            end
        end
        srt = sortperm(sds_s)
        sds_s = sds_s[srt]
        cols_s = cols_s[srt]
        ispri_s = ispri_s[srt]

        push!(sds, sds_s)
        push!(cols, cols_s)
        push!(ispri, ispri_s)
    end

    if shift_pris
        for j in 1:length(ss)
            xs = PlotTools.pointexclusion1D(sds[j], ysep, xsep)
            xs += ss[j]

            sdflt = sds[j] .> minsd
            if any(ispri[j][sdflt])
                xs_p = xs[sdflt][ispri[j][sdflt]]
                sds_p = sds[j][sdflt][ispri[j][sdflt]]
                cols_p = cols[j][sdflt][ispri[j][sdflt]]
                scatter(xs_p, sds_p, c=cols_p, marker="D", s=szp, zorder=10)
            end

            scatter(xs[sdflt][!ispri[j][sdflt]], sds[j][sdflt][!ispri[j][sdflt]], c=cols[j][sdflt][!ispri[j][sdflt]], marker="o", s=sz, zorder=10)
        end
    else
        for j in 1:length(ss)
            sdflt = sds[j] .> minsd

            if any(ispri[j][sdflt])
                xs = PlotTools.pointexclusion1D(sds[j][sdflt][ispri[j][sdflt]], ysep, xsep) + ss[j]
                scatter(xs, sds[j][sdflt][ispri[j][sdflt]], c=cols[j][sdflt][ispri[j][sdflt]], marker="D", s=szp, zorder=10)
            end

            xs = PlotTools.pointexclusion1D(sds[j][sdflt][!ispri[j][sdflt]], ysep, xsep)
            xs += ss[j]
            scatter(xs, sds[j][sdflt][!ispri[j][sdflt]], c=cols[j][sdflt][!ispri[j][sdflt]], marker="o", s=sz, zorder=10)
        end
    end

end

const ising_cft_1_1 = CFTField(0.0, Int[1, 0, 1, 1, 2, 2, 3, 3, 5, 5, 7, 8, 11, 12, 16, 18])
const ising_cft_1_2 = CFTField(1/16, Int[1, 1, 1, 2, 2, 3, 4, 5, 6, 8, 10, 12, 15, 18, 22, 27])
const ising_cft_1_3 = CFTField(1/2, Int[1, 1, 1, 1, 2, 2, 3, 4, 5, 6, 8, 9, 12, 14, 17, 20])

const ising_cft_eye = CFTOp(ising_cft_1_1, ising_cft_1_1, "\\mathbb{I}")
const ising_cft_sigma = CFTOp(ising_cft_1_2, ising_cft_1_2, "\\sigma")
const ising_cft_epsilon = CFTOp(ising_cft_1_3, ising_cft_1_3, "\\epsilon")

const potts3_cft_1_1 = CFTField(0.0, Int[1, 0, 1, 1, 2, 2, 4, 4, 7, 8, 12, 14, 21, 24, 34, 41])
const potts3_cft_1_2 = CFTField(1/8, Int[1, 1, 1, 2, 3, 4, 6, 8, 11, 15, 20, 26, 35, 45, 58, 75])
const potts3_cft_1_3 = CFTField(2/3, Int[1, 1, 2, 2, 4, 5, 8, 10, 15, 19, 27, 34, 46, 58, 77, 96])
const potts3_cft_1_4 = CFTField(13/8, Int[1, 1, 2, 3, 4, 6, 9, 12, 16, 22, 29, 38, 50, 64, 82, 105])
const potts3_cft_2_1 = CFTField(2/5, Int[1, 1, 1, 2, 3, 4, 6, 8, 11, 15, 20, 26, 35, 45, 58, 74])
const potts3_cft_2_3 = CFTField(1/15, Int[1, 1, 2, 3, 5, 7, 10, 14, 20, 26, 36, 47, 63, 81, 106, 135])
const potts3_cft_3_1 = CFTField(7/5, Int[1, 1, 2, 2, 4, 5, 8, 10, 15, 19, 26, 33, 45, 56, 74, 92])
const potts3_cft_4_1 = CFTField(3.0, Int[1, 1, 2, 3, 4, 5, 8, 10, 14, 18, 24, 31, 41, 51, 66, 83])

const potts3_cft_eye = CFTOp(potts3_cft_1_1, potts3_cft_1_1, "\\mathbb{I}")
const potts3_cft_sigma = CFTOp(potts3_cft_2_3, potts3_cft_2_3, "\\sigma")
const potts3_cft_epsilon = CFTOp(potts3_cft_2_1, potts3_cft_2_1, "\\epsilon")
const potts3_cft_Z = CFTOp(potts3_cft_1_3, potts3_cft_1_3, "Z")
const potts3_cft_X_eps = CFTOp(potts3_cft_3_1, potts3_cft_2_1, "(X,\\epsilon)")
const potts3_cft_eps_X = CFTOp(potts3_cft_2_1, potts3_cft_3_1, "(\\epsilon,X)")
const potts3_cft_X = CFTOp(potts3_cft_3_1, potts3_cft_3_1, "X")
const potts3_cft_Y = CFTOp(potts3_cft_4_1, potts3_cft_4_1, "Y")
const potts3_cft_Y_eye = CFTOp(potts3_cft_4_1, potts3_cft_1_1, "(Y,\\mathbb{I})")
const potts3_cft_eye_Y = CFTOp(potts3_cft_1_1, potts3_cft_4_1, "(\\mathbb{I},Y)")

end

ext: jl

max_stars_repo_name: amilsted/spinchain-virasoro-demo

max_issues_repo_name: amilsted/spinchain-virasoro-demo

max_forks_repo_name: amilsted/spinchain-virasoro-demo


text:
@testset "NSGA-II" begin
    @testset "nondominated sorting" begin
        ws=[-1, -1]
        let fA = [20, 2.2], fB = [60, 4.4], fD = [15, 4.4]
            @test !dominate(fA, fD, ws)
            @test !dominate(fD, fA, ws)
            @test !dominate(fA, fA, ws)
            @test dominate(fA, fB, ws)
            @test dominate(fD, fB, ws)
            @test !dominate(fB, fD, ws)
        end
        fitness = [20 2.2; 60 4.4; 65 3.5; 15 4.4; 55 4.5; 50 1.8; 80 4.0; 25 4.6]'
        fronts = nondominated_sort(fitness, ws)
        # the order may differ, though we don't care about the order within each front
        @test all(issetequal.(fronts, [[1, 4, 6], [2, 3, 5, 8], [7]]))
    end

    @testset "crowding distance assignment" begin
        fitness = [20 2.2; 60 4.4; 65 3.5; 15 4.4; 55 4.5; 50 1.8; 80 4.0; 25 4.6]'
        ws=[-1, -1]
        fronts = nondominated_sort(fitness, ws)
        cds = assign_crowding_distance(fitness, fronts)
        @test cds[1] == 2.0
        @test cds[4] == Inf
        @test cds[6] == Inf
        @test cds[8] == Inf
        @test cds[3] == Inf
        @test cds[5] ≈ 1.0568181818181819
        @test cds[2] ≈ 1.1590909090909092
        @test cds[7] == Inf
    end
    @testset "crowded comparison" begin
        @test crowded_compare(1, 0.1, 4, Inf) == -1
        @test crowded_compare(5, 0.1, 4, Inf) == 1
        @test crowded_compare(1, 0.1, 1, Inf) == 1
        @test crowded_compare(1, Inf, 1, Inf) == 0
    end
end


text:
module A5P1

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\A5P1.jl"),
  informat="script",
  out_path = "$(pwd())\\A5P1.html",
  doctype = "md2html")
=#
if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTModCopy, JuMP, NLopt, ForwardDiff, Roots

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const DATE_FORMAT_STR = "m/d/y"
const SP500_NAME = "GSPC"
const TBILL_NAME = "t-bill"
const TERM_NAME = "termStructure"
const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const MAX_YEAR = 2016
const IN_SAMPLE_CUTOFF_YEAR = 1975

const DEFAULT_NUM_SIMS = 10000
const DEFAULT_LOGμ_TOL = 10.0^-10.0
const PlotContainer = Union{Plot,Gadfly.Compose.Context}
const CSRegSpec = Tuple{Symbol, Symbol, Int}
const CPRegSpec = Tuple{Symbol, Int}
const TOLERANCE = 10.0^-8.0

const MIN_VALS_FOR_EST = 10

abstract type AbstractARCH end

function preProcessSP500()::Void
  dateFormat::String = DATE_FORMAT_STR
  sp500DF::DataFrame = readtable("$DATA_PATH\\$SP500_NAME.csv")

  #fix the dates
  rename!(sp500DF, :Date, :DATE)
  sp500DF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(sp500DF[:,:DATE])::DataVector{Date}
  sort!(sp500DF, cols = [:DATE])
  sp500DF[:,:returns] = similar(sp500DF[:,:Adj_Close])
  sp500DF[2:end,:returns] = log.(sp500DF[2:end,:Adj_Close] ./  sp500DF[1:(end-1),:Adj_Close])

  sp500DF[:,:dayOfQuarter] = zeros(Int,size(sp500DF,1))
  for i::Int ∈ 2:size(sp500DF,1)
    if Dates.firstdayofquarter(sp500DF[i,:DATE]) ≠ Dates.firstdayofquarter(sp500DF[i-1,:DATE])
      sp500DF[i,:dayOfQuarter] = 1
    else
      sp500DF[i,:dayOfQuarter] = sp500DF[i-1,:dayOfQuarter] + 1
    end
  end
  #showcols(sp500DF)

  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls", "w")
  serialize(stream, sp500DF)
  close(stream)

  return nothing
end

#holds a GARCH parameter object
mutable struct ΘGARCH <: AbstractARCH
  λ::Real
  β::Real
  α::Real
  μ::Real
  valid::Bool

end

#constructor
ΘGARCH(λ::Real, β::Real, α::Real, μ::Real)::ΘGARCH = ΘGARCH(λ, β, α, μ, true)

mutable struct ΘEGARCH <: AbstractARCH
  λ::Real
  β::Real
  α::Real
  μ::Real
  γ::Real
  valid::Bool

end

#constructor
ΘEGARCH(λ::Real, β::Real, α::Real, μ::Real, γ::Real)::ΘEGARCH = ΘEGARCH(λ, β, α, μ, γ, true)

function estimateGARCH(r::Vector{Float64}, verbose::Bool = false)::ΘGARCH

  #for convenience
  T::Int = length(r)

  function obj(λ::Real, β::Real, α::Real, μ::Real)

    #σ²::Vector{Real} = Vector{Real}(T-1)
    σ²tM1::Real = λ/(1.0-β)
    #l::Float64 = -(r[1]-μ)*(r[1]-μ)/(2.0*σ²₀)-0.5*log(σ²₀) #initial likelihood
    l::Real = 0.0

    #calculate the likelihood and get the next volatiltiy in the sequence
    for i::Int ∈ 1:T
      u²::Real = (r[i]-μ)*(r[i]-μ)
      l += -u²/σ²tM1 - 0.5*log(σ²tM1)
      σ²tM1 = λ+β*σ²tM1+α*u²
    end

    #println("λ=$λ, β=$β, α=$α, μ=$μ")

    return l
  end

  meanRet = mean(r)

  #set initial meta parameters
  minVal = 10.0^-10.0
  maxVal = 1.0-minVal

  #Algorithms: :LD_LBFGS, :LN_COBYLA, :LN_BOBYQA
  #  :LD_SLSQP, :LD_TNEWTON, :LD_VAR2/:LD_VAR1,
  #  :LD_MMA, :GN_ESCH, :GN_DIRECT, GN_DIRECT_L, :GD_STOGO
  # GN_ISRES, :LN_NELDERMEAD, :LN_SBPLX,
  #GN_CRS2_LM
  #:GN_CRS2_LM***, :GN_ISRES*** (slow), :LD_MMA**, LN_SBPLX**, LD_TNEWTON** best:62,848
  #MMA: λ=2.991604487710108e-6, β=0.8691809689019686, α=0.2447041114406467, μ=0.0005059861641042039
  #NEWTON: 3.7136249890797323e-6, β=0.8641416287811661, α=0.23581482405030202, μ=0.0005248887137138578
  #CRS2 λ=1.8304983359585568e-6, β=0.9091365756370056, α=0.1664779062477393, μ=0.00048124017464920696
  #SBPLX λ=1.7187359905806954e-6, β=0.9122999784465299, α=0.16106473023211287, μ=0.00029331108088162115
  #LR Vol>0 #:LN_COBYLA(61450), LD_MMA(61144), ISRES(62578),

  mod = Model(solver=NLoptSolver(algorithm=:LD_TNEWTON_PRECOND,maxtime=10))

  @variable(mod, minVal <= λ <= maxVal)
  @variable(mod, minVal <= β <= maxVal)
  @variable(mod, minVal <= α <= maxVal)
  @variable(mod, minVal <= μ <= maxVal)
  #@NLconstraint(mod, 1.0-α-β>=0.0)

  setvalue(λ, 0.5)
  setvalue(β, 0.5)
  setvalue(α, 0.5)
  setvalue(μ, mean(r)) #guess the unconditional mean

  JuMP.register(mod, :obj, 4, obj, autodiff=true)
  @NLobjective(mod, Max, obj(λ, β, α, μ))

  ΘBest::ΘGARCH  = ΘGARCH(0.0, 0.0, 0.0, 0.0, false)
  objBest::Float64 = -9999.0
  try
    output = solve(mod)
    ΘBest  = ΘGARCH(getvalue(λ), getvalue(β), getvalue(α), getvalue(μ))
    objBest= obj(getvalue(λ), getvalue(β), getvalue(α), getvalue(μ))
  end

  if verbose
    impliedAnnualLRVol::Float64 = (ΘBest.λ/(1.0 - ΘBest.β)*255)^0.5
    historicalLRVol::Float64 = std(r)*255.0^0.5

    impliedReturn::Float64 = exp(ΘBest.μ*255)-1.0
    historicalReturn::Float64 = exp(mean(r)*255)-1.0

    println("λ=$(ΘBest.λ), β=$(ΘBest.β), α=$(ΘBest.α), μ=$(ΘBest.μ)")
    println("Vol: Implied=$impliedAnnualLRVol, Hist=$historicalLRVol")
    println("Return: Implied=$impliedReturn, Hist=$historicalReturn")
    println("Objective: $objBest")
    println("For Reference: x̄=$(mean(r))")
  end

  return ΘBest

end

function estimateEGARCH(r::Vector{Float64}, verbose::Bool = false)::ΘEGARCH

  T::Int = length(r)

  function obj(λ::Real, β::Real, α::Real, μ::Real, γ::Real)

    #σ²::Vector{Real} = Vector{Real}(T-1)
    EAbsZ::Float64 = (2.0/π)^0.5
    lσ²tM1::Real = (λ-γ*EAbsZ)/(1.0-β)
    #lσ²tM1::Real = λ/(1.0-β)

    l::Real = 0.0

    #calculate the likelihood and get the next volatiltiy in the sequence
    for i::Int ∈ 1:T
      u²::Real = (r[i]-μ)*(r[i]-μ)
      σtM1::Real = (exp(lσ²tM1))^0.5
      l += -u²/(σtM1*σtM1) - 0.5*lσ²tM1
      lσ²tM1 = λ+β*lσ²tM1+α*u²+γ*(abs((r[i]-μ)/σtM1)-EAbsZ)
    end

    return l
  end

  meanRet = mean(r)

  #set initial meta parameters
  minVal = -10.0
  maxVal = 100.0

  #Algorithms: :LD_LBFGS, :LN_COBYLA, :LN_BOBYQA
  #  :LD_SLSQP, :LD_TNEWTON, :LD_VAR2/:LD_VAR1,
  #  :LD_MMA, :GN_ESCH, :GN_DIRECT, GN_DIRECT_L, :GD_STOGO
  # GN_ISRES, :LN_NELDERMEAD, :LN_SBPLX,
  #GN_CRS2_LM
  #:LD_LBFGS(62774), LD_VAR2(62774), LN_NELDERMEAD(61880), :LN_SBPLX (62791), :LD_TNEWTON(62800)

  mod = Model(solver=NLoptSolver(algorithm=:LD_TNEWTON_PRECOND, maxtime=20))

  @variable(mod, minVal <= λ <= maxVal)
  @variable(mod, minVal <= β <= maxVal)
  @variable(mod, minVal <= α <= maxVal)
  @variable(mod, 0.0 <= μ <= maxVal)
  @variable(mod, minVal <= γ <= maxVal)
  #@NLconstraint(mod, 1.0-α-β>=0.0)

  setvalue(λ, 0.5)
  setvalue(β, 0.5)
  setvalue(α, 0.5)
  setvalue(μ, mean(r)) #guess the unconditional mean
  setvalue(γ, 0.5) #guess the unconditional mean

  JuMP.register(mod, :obj, 5, obj, autodiff=true)
  @NLobjective(mod, Max, obj(λ, β, α, μ, γ))

  ΘBest::ΘEGARCH  = ΘEGARCH(getvalue(λ), getvalue(β), getvalue(α), getvalue(μ), getvalue(γ), false)
  objBest::Float64 = -9999.0
  try
    output = solve(mod)
    ΘBest  = ΘEGARCH(getvalue(λ), getvalue(β), getvalue(α), getvalue(μ), getvalue(γ))
    objBest= obj(getvalue(λ), getvalue(β), getvalue(α), getvalue(μ), getvalue(γ))
  end

  if verbose
    impliedAnnualLRVol::Float64 = (exp((ΘBest.λ-ΘBest.γ*(2/π)^0.5)/(1.0-ΘBest.β))*255)^0.5
    #impliedAnnualLRVol::Float64 = (exp(ΘBest.λ)/(1-ΘBest.β)*255)^0.5
    historicalLRVol::Float64 = std(r)*255.0^0.5

    impliedReturn::Float64 = exp(ΘBest.μ*255)-1.0
    historicalReturn::Float64 = exp(mean(r)*255)-1.0

    println("λ=$(ΘBest.λ), β=$(ΘBest.β), α=$(ΘBest.α), μ=$(ΘBest.μ), γ=$(ΘBest.γ)")
    println("Vol: Implied=$impliedAnnualLRVol, Hist=$historicalLRVol")
    println("Return: Implied=$impliedReturn, Hist=$historicalReturn")
    println("Objective: $objBest")
    println("For Reference: x̄=$(mean(r))")
  end

  return ΘBest

end

function estimateGARCHByQ(sp500DF::DataFrame, estimateModel::Function)#::T where T <: AbstractARCH

  dateDict::Dict = Dict(sp500DF[i,:DATE]=>i for i::Int ∈ 1:size(sp500DF,1))

  minValsForEst::Int = MIN_VALS_FOR_EST
  r::Vector{Float64} = sp500DF[2:end,:returns]

  #get a list of target indices
  quarterlyDates::Vector{Date} = sp500DF[sp500DF[:,:dayOfQuarter].==1,:DATE]
  terminalIndicesToEst::Vector{Int} = ((d::Date)->dateDict[d]).(quarterlyDates)
  #terminalIndicesToEst = terminalIndicesToEst[terminalIndicesToEst[:].>minValsForEst]

  #set up a container of contianers
  estimated::Vector{Vector{T where T<:AbstractARCH}} =
    [Vector{T where T<:AbstractARCH}() for i::Int ∈ 1:length(terminalIndicesToEst)]

  println("beginning $(length(estimated)) simulations")
  fails::Int = 0
  @fastmath for i ∈ 1:length(terminalIndicesToEst)

  #estimate the grabbed indices
    estimated[i] = [estimateModel(r[1:(terminalIndicesToEst[i]-1)], false)]
    if estimated[i][1].valid == false
      fails += 1
    end
  end
  println("simulations complete. $fails failures out of $(length(estimated)) simulations.")

  return [estimated[i][1] for i::Int ∈ 1:length(terminalIndicesToEst)]

end

#script function for running the GARCH problem
function A5P1Script(;refreshData::Bool=true, refreshGARCH::Bool=true, refreshEGARCH::Bool=true,)
  if !isfile("$DATA_PATH\\$SP500_NAME.jls") || refreshData
    preProcessSP500()
  end

  #read the binary file
  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls")
  sp500DF::DataFrame = deserialize(stream)
  close(stream)

  r::Vector{Float64} = sp500DF[2:end,:returns]
  quarterlyDates::Vector{Date} = sp500DF[sp500DF[:,:dayOfQuarter].==1,:DATE]

  #estimate the GARCH model and print the output

  estimateGARCH(r, true)
  estimateEGARCH(r, true)

  #this setup is so that we don't ahve to solve the model each quarter repeatedly
  if refreshGARCH
    estimatedGARCH::Vector{ΘGARCH} = estimateGARCHByQ(sp500DF, estimateGARCH)

    #cache the results on disk
    stream = open("$OUTPUT_PATH\\GARCH_RESULTS.jls", "w")
    serialize(stream, estimatedGARCH)
    close(stream)
  else

    #access the resutls
    stream = open("$OUTPUT_PATH\\GARCH_RESULTS.jls")
    estimatedGARCH = deserialize(stream)
    close(stream)
  end

  if refreshEGARCH
    estimatedEGARCH::Vector{ΘEGARCH} = estimateGARCHByQ(sp500DF, estimateEGARCH)

    #cache the results on disk
    stream = open("$OUTPUT_PATH\\EGARCH_RESULTS.jls", "w")
    serialize(stream, estimatedEGARCH)
    close(stream)
  else

    #access the resutls
    stream = open("$OUTPUT_PATH\\EGARCH_RESULTS.jls")
    estimatedEGARCH = deserialize(stream)
    close(stream)
  end

  ##############Now work on the outputs. First make them into a dataframe
  garchDF::DataFrame = DataFrame(DATE=quarterlyDates, λ=[g.λ for g ∈ estimatedGARCH], β=[g.β for g ∈ estimatedGARCH],
    α=[g.α for g ∈ estimatedGARCH], μ=[g.μ for g ∈ estimatedGARCH],
    μLR=[exp(255*g.μ)-1.0 for g ∈ estimatedEGARCH], σ=[(g.λ/(1.0 - g.β)*255)^0.5 for g ∈ estimatedGARCH])

  egarchDF::DataFrame  = DataFrame(DATE=quarterlyDates, λ=[g.λ for g ∈ estimatedEGARCH], β=[g.β for g ∈ estimatedEGARCH],
    α=[g.α for g ∈ estimatedEGARCH], μ=[g.μ for g ∈ estimatedEGARCH], γ=[g.γ for g ∈ estimatedEGARCH],
    μLR=[exp(255*g.μ)-1.0 for g ∈ estimatedEGARCH], σ=[(exp((g.λ-g.γ*(2/π)^0.5)/(1.0-g.β))*255)^0.5 for g ∈ estimatedEGARCH])

  nQuarters::Int = length(quarterlyDates)

  for i::Int∈1:nQuarters
    if estimatedGARCH[i].valid == false
      garchDF[i,[:λ, :β, :α, :μ, :σ]] = NA
    end

    if estimatedEGARCH[i].valid == false
      egarchDF[i,[:λ, :β, :α, :μ,:γ, :σ]] = NA
    end
  end

  Gadfly.push_theme(:default)
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  push!(plotNames, "GarchParameters")
  push!(plots, vstack(plot(x=garchDF[:,:DATE], y=garchDF[:,:λ],
    Geom.point, #Theme(line_width = 0.25pt),
    Guide.ylabel("Value"), Guide.xlabel("Date"),
    Guide.title("GARCH λ Param")),
    plot(x=garchDF[:,:DATE], y=garchDF[:,:β],
       Geom.point, #Theme(line_width = 0.25pt),
      Guide.ylabel("Value"), Guide.xlabel("Date"),
      Guide.title("GARCH β Param")),
    plot(x=garchDF[:,:DATE], y=garchDF[:,:α],
      Geom.point, #Theme(line_width = 0.25pt),
      Guide.ylabel("Value"), Guide.xlabel("Date"),
      Guide.title("GARCH α Param"))))

  egarchDF = egarchDF[((!isna).(egarchDF[:,:σ])),:]

  push!(plotNames, "EGarchParameters")
  push!(plots, vstack(plot(x=egarchDF[:,:DATE], y=egarchDF[:,:λ],
    Geom.point, #Theme(line_width = 0.25pt),
    Guide.ylabel("Value"), Guide.xlabel("Date"),
    Guide.title("GARCH λ Param")),
    plot(x=egarchDF[:,:DATE], y=egarchDF[:,:β],
       Geom.point, #Theme(line_width = 0.25pt),
      Guide.ylabel("Value"), Guide.xlabel("Date"),
      Guide.title("GARCH β Param")),
    plot(x=egarchDF[:,:DATE], y=egarchDF[:,:α],
      Geom.point, #Theme(line_width = 0.25pt),
      Guide.ylabel("Value"), Guide.xlabel("Date"),
      Guide.title("GARCH α Param"))))

  push!(plotNames, "GarchVol")
  push!(plots, vstack(plot(x=garchDF[:,:DATE], y=garchDF[:,:μLR],
    Geom.point, #Theme(line_width = 0.25pt),
    Guide.ylabel("Value"), Guide.xlabel("Date"), Coord.Cartesian(ymin=0.0,ymax=0.5),
    Guide.title("GARCH μ Param")),
    plot(x=garchDF[:,:DATE], y=garchDF[:,:σ],
       Geom.point, #Theme(line_width = 0.25pt),
      Guide.ylabel("Value"), Guide.xlabel("Date"),
      Guide.title("GARCH σ Param"))))

  push!(plotNames, "EGarchVol")
  push!(plots, vstack(plot(x=egarchDF[:,:DATE], y=egarchDF[:,:μLR],
    Geom.point, #Theme(line_width = 0.25pt),
    Guide.ylabel("Value"), Guide.xlabel("Date"), Coord.Cartesian(ymin=0.0,ymax=0.5),
    Guide.title("GARCH μ Param")),
  plot(x=egarchDF[:,:DATE], y=egarchDF[:,:σ],
     Geom.point, #Theme(line_width = 0.25pt),
    Guide.ylabel("Value"), Guide.xlabel("Date"),
    Guide.title("GARCH σ Param"))))

  for i ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 7.5inch, 9inch),plots[i])
  end

  #graph the GARCH model

end

@time begin
  #A5P1Script(refreshData=false, refreshGARCH=false, refreshEGARCH=false)
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


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

using Markdown
using InteractiveUtils

# ╔═╡ 81cdc1aa-04fd-11eb-3235-bff0e3f52535
using Pkg, DrWatson

# ╔═╡ 81cdee8c-04fd-11eb-29f5-85927396b6dc
begin
	@quickactivate "SR2TuringPluto"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ 5e06a11c-04f6-11eb-0116-1fd51d614f99
md"## clip-04-64-79t"

# ╔═╡ b8b8de5e-0512-11eb-0932-830c3055e00a
#include(srcdir("tools.jl"))

md"### snippet 4.64"

# ╔═╡ 81ce6d32-04fd-11eb-2c4f-77bdf6a3f9e2
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	scale!(df, [:weight])
end;

# ╔═╡ 81e16d54-04fd-11eb-26e6-5b62d66d3d34
f_parabola(weight_s, a, b1, b2) = a + b1 * weight_s + b2 * weight_s^2

# ╔═╡ 81ee40ec-04fd-11eb-1c71-6b4674695b56
@model function parabola(weight_s, heights)
    a ~ Normal(178, 20)
    b1 ~ LogNormal(0, 1)
    b2 ~ Normal(0, 1)
    σ ~ Uniform(0, 50)
    μ = f_parabola.(weight_s, a, b1, b2)
    heights ~ MvNormal(μ, σ)
end

# ╔═╡ 81f015b6-04fd-11eb-23c4-3f3f8a7201c8
begin
	quap4_5t = quap(parabola(df.weight_s, df.height), NelderMead())
	quap4_5t.coef
end

# ╔═╡ 960b36fe-0573-11eb-0e51-9f7bc0d175a5
quap4_5t.vcov

# ╔═╡ 81f6618c-04fd-11eb-27e7-dffa73914969
md"### snippet 4.67"

# ╔═╡ 3228440c-0527-11eb-00f0-9b38ae90372c
begin
	weight_seq_1 = range(-2.2, 2, length = 30)
	post_1 = DataFrame(rand(quap4_5t.distr, 1_000)', quap4_5t.params)
	mu_1 = f_parabola.(weight_seq_1, post_1.a', post_1.b1', post_1.b2')
	sim_1 = rand.(Normal.(mu_1, post_1.σ'))
	mu_1 = meanlowerupper(mu_1)
	sim_1 = meanlowerupper(sim_1)
end;

# ╔═╡ 82285bce-04fd-11eb-0bfd-27cbf60511ee
md"### snippet 4.68"

# ╔═╡ 93badee6-0527-11eb-16bc-b1c6541276ba
begin
	scatter(df.weight_s, df.height, ms = 3, alpha = 0.7, legend = false)
	plot!(weight_seq_1, mu_1.mean, ribbon = (mu_1.mean .- mu_1.lower, mu_1.upper .- mu_1.mean))
	plot!(weight_seq_1, sim_1.lower, fillrange = sim_1.upper, alpha = 0.3, linealpha = 0.0, c = 2)
end

# ╔═╡ 4f386d5c-0573-11eb-084c-1758217d06f9
md"## End of clip-04-64-68t.jl"

# ╔═╡ Cell order:
# ╟─5e06a11c-04f6-11eb-0116-1fd51d614f99
# ╠═81cdc1aa-04fd-11eb-3235-bff0e3f52535
# ╠═81cdee8c-04fd-11eb-29f5-85927396b6dc
# ╟─b8b8de5e-0512-11eb-0932-830c3055e00a
# ╠═81ce6d32-04fd-11eb-2c4f-77bdf6a3f9e2
# ╠═81e16d54-04fd-11eb-26e6-5b62d66d3d34
# ╠═81ee40ec-04fd-11eb-1c71-6b4674695b56
# ╠═81f015b6-04fd-11eb-23c4-3f3f8a7201c8
# ╠═960b36fe-0573-11eb-0e51-9f7bc0d175a5
# ╟─81f6618c-04fd-11eb-27e7-dffa73914969
# ╠═3228440c-0527-11eb-00f0-9b38ae90372c
# ╠═82285bce-04fd-11eb-0bfd-27cbf60511ee
# ╠═93badee6-0527-11eb-16bc-b1c6541276ba
# ╟─4f386d5c-0573-11eb-084c-1758217d06f9


text:
module MyOptim

using LinearAlgebra
using Suppressor

function conjgrad(fn, grad, x0; maxstep=2.0, niters=50, conv_thr = 1e-2, fn_conv = 1e-4, verbosity="low")

    println()
    println("Conj Grad START")
    x = deepcopy(x0)
    storage = zeros(size(x))

    xold = deepcopy(x)

#    println("x0 ", x)

#    println("CG RUN GRAD1")
    #iter 1
    f, g = grad(storage, x)

    fold1 = 1e8
    fold2 = 1e8

#    println("f ", f, " x   ",  x)

    dx = deepcopy(-g)

    dx1 = deepcopy(dx)

    s = deepcopy(dx)

    step_size = maxstep/10.0

#    println("CG RUN LS1")

    old_step=step_size
    x, step_size = linesearch(x, dx, fn, f, step_size, verbosity)
    step_size = min(step_size, maxstep)

    for i = 1:niters

        if sum(abs.(xold - x)) > 1e-7  #only update grad if x changes

            fold2 = fold1
            fold1 = f

            f, g = grad(storage, x)
            xold = deepcopy(x)
#        else
#            println("nochange to grad")
        end

        if (sqrt(sum(g.^2)) < conv_thr && f < 0.1) || ( abs(fold1 - fold2) < fn_conv &&  abs(fold1 - f) < fn_conv)
            println("yes conv sum_grad ", sqrt(sum(g.^2)), " fn_diff_1 ", abs(fold1 - fold2), " fn_diff_2 ", abs(fold1 - f))
            return x, f, g
#        else
#            println("no conv ", sqrt(sum(g.^2)), " ", conv_thr, " f ", f)
        end

#        println("MY CONJGRAD $i ", f, " sg ", sqrt(sum(g.^2)), " step_size $step_size --------------------")

        dx1[:] = dx[:]
        dx[:] = -g[:]

#        println("dx ", dx)
#        println("dx1 ", dx1)

        beta = ( dx' * (dx - dx1)) / (dx1' * dx1)  #PR
#        println("beta $beta")
        beta = max(0, beta)

        s = dx + beta * s

        old_step=step_size

        x, step_size, good = linesearch(x, s, fn, f, step_size, verbosity)
        if good
            step_size =min(step_size, maxstep)
        else
            beta = 0.0
        end
        println("MY CG quadratic linesearch iter $i fn val: $f | rms grad: ", sqrt(sum(g.^2)), " within_stepsize: $good , old step: $old_step , new step_size: $step_size ")

    end

    println("WARNING, conv not achieved")
    return x, f, g

end

function linesearch(x, dx, fn, f0, step_size, verbosity)
#    println("CG FN0 $f0")
#    println("CG dx $dx")
#    println("ls $x $dx $f0 $step_size")
#    println("CG FN1 ss $step_size")

    if verbosity=="low"
        @suppress f0r, flag0r = fn( x + dx * step_size * 0.0)
    else
        f0r, flag0r = fn( x + dx * step_size * 0.0)
    end

#    println("CG real $f0r $flag0r")

    if verbosity=="low"
        @suppress f1, flag1 = fn( x + dx * step_size * 0.5)
    else
        f1, flag1 = fn( x + dx * step_size * 0.5)
    end

#    println("CG FN1 $f1 $flag1")
    if f1 > f0 || flag1 == false
#        println("MY LS uphill $f1 > $f0 flag $flag1, reduce step")
        return x, step_size/3.1, false
    end
#    println("CG FN2")

    if verbosity == "low"
        @suppress f2, flag2 = fn( x + dx * step_size)
    else
        f2, flag2 = fn( x + dx * step_size)
    end

#    println("CG FN2 $f2 $flag2")

    if flag2 == false
#        println("MY LS flag2 $flag2")
        return x + dx * step_size * 0.5, step_size/1.75, false
    end

    c = f0
    a = 2 * (f2-2*f1+f0)
    b = f2 - a - c

#    println("a b c $a $b $c ", -b/(2*a))

    if a < 0
#        println("MY LS a negative $a, take max step increase stepsize")
        return x + dx * step_size * 1.0, step_size * 2.0, true
    else
        step = min(-b/(2*a), 1.0)
#        println("MY LS NORMAL $step")
        return x + dx * step * step_size, step * step_size * 1.5, true
    end

end

end #end module

#=
function myf(x)

    return -0.5*x'*A*x + 0.1*sum(x.^4)

end

using ForwardDiff

function grad(storage, x)
    sleep(0.1)

    return myf(x), ForwardDiff.gradient(myf, x)

end

function grad_only(storage, x)
    storage[:] = ForwardDiff.gradient(myf, x)
    sleep(0.1)
    return storage
end
=#


text:
"""
  localMass2D!()

Computes local mass matrix for 2D geometries, e.g. for use in discretized time derivative
"""
function localMass2D!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                       At,phi,dphidx,dphidy,dphids,dphidt,order,parameter...)
  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(Float64))

  for gpt=1:nGaussNodes
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
    jac = jacCalc(xN,yN,dphids,dphidt)

    for tj = 1:nNodesPerElm, ti = 1:nNodesPerElm
      integrand = phi[ti]*phi[tj]*jac
      At[ti,tj] += integrand*w[gpt]
    end
  end

  nothing
end

"""
  localMassWithScalar2D!()

Computes local mass matrix for 2D geometries, e.g. for use in discretized time derivative
"""
function localMassWithScalar2D!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                       At,phi,dphidx,dphidy,dphids,dphidt,order,param::Vector{T}) where T<:Real
  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(Float64))

  αNodes = zeros(T,nGaussNodes)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(T))
  for i=1:nGaussNodes
    αNodes[i] = param[mesh.cm[el].NodeList[i]]
  end

  for gpt=1:nGaussNodes
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
    jac = jacCalc(xN,yN,dphids,dphidt)
    αg  = shapeEval(αNodes,phi)

    for tj = 1:nNodesPerElm, ti = 1:nNodesPerElm
      integrand = αg*phi[ti]*phi[tj]*jac
      At[ti,tj] += integrand*w[gpt]
    end
  end

  nothing
end

"""
  localLaplace2D!()

Computes local laplace matrix for 2D geometries.

ASSUMPTIONS:
- doesn't use parameter! Assumes any parameter information is put into forcing function

Should really be called poisson2D, but whatever
"""
function localLaplace2D!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                       At,phi,dphidx,dphidy,dphids,dphidt,order,param)
  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(Float64))

  for gpt=1:nGaussNodes
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
    jac = jacCalc(xN,yN,dphids,dphidt)

    for ti = 1:nNodesPerElm, tj = 1:nNodesPerElm
      integrand = (dphidx[ti]*dphidx[tj] + dphidy[ti]*dphidy[tj])/jac
      At[ti,tj] += integrand*w[gpt]
    end
  end

  nothing
end

"""
localDarcy2D!()

Computes local Darcy matrix for 2D geometries.
"""
function localDarcy2D!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
  At,phi,dphidx,dphidy,dphids,dphidt,order,param::Vector{T}) where T<:Real
  αNodes = zeros(T,nGaussNodes)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(T))
  for i=1:nGaussNodes
    αNodes[i] = param[mesh.cm[el].NodeList[i]]
  end

  for gpt=1:nGaussNodes
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
    jac = jacCalc(xN,yN,dphids,dphidt)
    αg  = shapeEval(αNodes,phi)

    for tj = 1:nNodesPerElm, ti = 1:nNodesPerElm
    integrand = (dphidx[ti]*dphidx[tj] + dphidy[ti]*dphidy[tj])/(jac*αg)
    At[ti,tj] += integrand*w[gpt]
    end
  end

  nothing
end

# Advection part of Advection-Diffusion Equation
function localAdvDiff2D!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                       At,phi,dphidx,dphidy,dphids,dphidt,order,parameter)
  uNodes = zeros(Float64,nGaussNodes)
  vNodes = zeros(Float64,nGaussNodes)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(Float64))
  for i=1:nNodesPerElm
    uNodes[i] = parameter.u[mesh.cm[el].NodeList[i]]
    vNodes[i] = parameter.v[mesh.cm[el].NodeList[i]]
  end

  for gpt=1:nGaussNodes
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
    jac = jacCalc(xN,yN,dphids,dphidt)
    ug  = shapeEval(uNodes,phi)
    vg  = shapeEval(vNodes,phi)

    for tj = 1:nNodesPerElm, ti = 1:nNodesPerElm
	    integrand =  ug*phi[ti]*dphidx[tj] + vg*phi[ti]*dphidy[tj]
      At[ti,tj] += integrand*w[gpt]
    end
  end

  nothing
end

# 2D SUPG stablization term for 2D Advection-diffusion equation
function localAdvDiff2D_SUPGmod!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                                 At,phi,dphidx,dphidy,dphids,dphidt,order,parameter)
  uNodes = zeros(Float64,nGaussNodes)
  vNodes = zeros(Float64,nGaussNodes)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(Float64))
  for i=1:nNodesPerElm
    uNodes[i] = parameter.u[mesh.cm[el].NodeList[i]]
    vNodes[i] = parameter.v[mesh.cm[el].NodeList[i]]
  end

  # compute h (max element length)
  hk = sqrt(areaCalc(xN,yN))

  # compute element Peclet number
  wmag = 0.0
  elPe = 0.0
  if mesh.order == :Linear
    uavg = sum(uNodes)/4
    vavg = sum(vNodes)/4
    wmag = sqrt(uavg^2 + vavg^2)

    elPe = 0.5*wmag*hk/parameter.κ
  elseif mesh.order == :Quadratic
    wmag = sqrt(uNodes[9]^2 + vNodes[9]^2)

    elPe = 0.5*wmag*hk/parameter.κ
  else
    error("mesh.order not specified correctly")
  end

  # compute delta for SUPG
  δ = 0.5*hk/wmag*(1.0 - 1.0/elPe)

  if mesh.order == :Linear
    for gpt=1:nGaussNodes
      shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
      derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
      jac = jacCalc(xN,yN,dphids,dphidt)
      ug  = shapeEval(uNodes,phi)
      vg  = shapeEval(vNodes,phi)

      for tj = 1:nNodesPerElm, ti = 1:nNodesPerElm
        integrand =  δ*(ug*dphidx[tj] + vg*dphidy[tj])*(ug*dphidx[ti] + vg*dphidy[ti])
        At[ti,tj] += integrand*w[gpt]
      end
    end
  elseif mesh.order == :Quadratic

  end
  nothing
end

# Stokes Constant viscosity
function localStokesConst2D!(mesh,el,xN,yN,w,s,t,Axt,Ayt,Bxt,Byt,
                        phi,dphidx,dphidy,dphids,dphidt,
                        psi,dpsidx,dpsidy,dpsids,dpsidt,
                        parameter...)
  μ0 = parameter[1]
  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)

    for ti=1:9
      for tj=1:9
		      integrandX = μ0*dphidx[ti]*dphidx[tj]/jac +
                       μ0*dphidy[ti]*dphidy[tj]/jac
		      Axt[ti,tj] += integrandX*w[gpt]

		      integrandY = μ0*dphidx[ti]*dphidx[tj]/jac +
                       μ0*dphidy[ti]*dphidy[tj]/jac
		      Ayt[ti,tj] += integrandY*w[gpt]
	    end

	    for tk=1:4
	      integrandX = -psi[tk]*dphidx[ti]
	      Bxt[tk,ti] += integrandX*w[gpt]

	      integrandY = -psi[tk]*dphidy[ti]
	      Byt[tk,ti] += integrandY*w[gpt]
	    end
    end
  end
end

# Stokes -- also correct!! -- Variable Viscosity
function localStokesVar2D_WRONG!(mesh,el,xN,yN,w,s,t,Axt,Ayt,A12t,A21t,Bxt,Byt,
                            phi,dphidx,dphidy,dphids,dphidt,
                            psi,dpsidx,dpsidy,dpsids,dpsidt,
                            parameter...)
  μ = parameter[1]

  μNodes = zeros(Float64,9)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  for i=1:9
    μNodes[i] = μ[mesh.cm[el].NodeList[i]]
  end
  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(A12t,zero(Float64))
  fill!(A21t,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)
    μg  = shapeEval(μNodes,phi)

    for ti=1:9
      for tj=1:9
        integrandX = 2*μg*dphidx[ti]*dphidx[tj]/jac +
                       μg*dphidy[ti]*dphidy[tj]/jac
        Axt[ti,tj] += integrandX*w[gpt]

        integrand12 = μg*dphidx[ti]*dphidy[tj]/jac
        A12t[ti,tj] += integrand12*w[gpt]

        integrand21 = μg*dphidy[ti]*dphidx[tj]/jac
        A21t[ti,tj] += integrand21*w[gpt]

        integrandY =   μg*dphidx[ti]*dphidx[tj]/jac +
                     2*μg*dphidy[ti]*dphidy[tj]/jac
        Ayt[ti,tj] += integrandY*w[gpt]
    end

    for tk=1:4
      integrandX = psi[tk]*dphidx[ti]
      Bxt[tk,ti] -= integrandX*w[gpt]

      integrandY = psi[tk]*dphidy[ti]
      Byt[tk,ti] -= integrandY*w[gpt]
    end
    end
  end

  nothing
end
# Stokes -- correct!! -- Variable Viscosity
function localStokesVar2D!(mesh,el,xN,yN,w,s,t,
                           At,Bt,Ct,Dt,Et,Ft,Gt,Ht,
                           phi,dphidx,dphidy,dphids,dphidt,
                           psi,dpsidx,dpsidy,dpsids,dpsidt,
                           parameter...)
  μ = parameter[1]

  μNodes = zeros(Float64,9)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  for i=1:9
    μNodes[i] = μ[mesh.cm[el].NodeList[i]]
  end

  fill!(At,zero(Float64))
  fill!(Bt,zero(Float64))
  fill!(Ct,zero(Float64))
  fill!(Et,zero(Float64))
  fill!(Ft,zero(Float64))
  fill!(Gt,zero(Float64))
  fill!(Ht,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)
    μg   = shapeEval(μNodes,phi)

    for ti=1:9
      # vector-valued test function
      for tj=1:9
        # integrand for A
        iA = 2*μg*dphidx[ti]*dphidx[tj]/jac +
               μg*dphidy[ti]*dphidy[tj]/jac
        At[ti,tj] += iA*w[gpt]

        # integrand for B
        iB = μg*dphidx[ti]*dphidy[tj]/jac
        Bt[ti,tj] += iB*w[gpt]

        # integrand for D -- transpose of B
        iD = μg*dphidx[tj]*dphidy[ti]/jac # iB
        Dt[ti,tj] += iD*w[gpt]        # transpose of B so i and j are switched

        # integrand for E
        iE =   μg*dphidx[ti]*dphidx[tj]/jac +
             2*μg*dphidy[ti]*dphidy[tj]/jac
        Et[ti,tj] += iE*w[gpt]
      end

      # scalar-valued test function
      for tj=1:4
        # integrand for C
        iC = psi[tj]*dphidx[ti] # C is Nu x Np block matrix, so i and j are correct
        Ct[ti,tj] -= iC*w[gpt]

        # integrand for F
        iF = psi[tj]*dphidy[ti] # F is Nu x Np block matrix, so i and j are correct
        Ft[ti,tj] -= iF*w[gpt]

        # integrand for G -- transpose of C
        iG = iC
        Gt[tj,ti] -= iG*w[gpt]    # G is Np x Nu block matrix, so i and j are switched

        # integrand for H -- transpose of H
        iH = iF
        Ht[tj,ti] -= iH*w[gpt]    # H is Np x Nu block matrix, so i and j are switched
      end
    end
  end

  nothing
end
# Stokes -- correct!! -- Variable Viscosity
function localStokesAlternateVar2D!(mesh,el,xN,yN,w,s,t,
                                    At,Bt,Ct,Dt,Et,Ft,Gt,Ht,
                                    phi,dphidx,dphidy,dphids,dphidt,
                                    psi,dpsidx,dpsidy,dpsids,dpsidt,
                                    parameter...)
  μ = parameter[1]

  μNodes = zeros(Float64,9)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  for i=1:9
    μNodes[i] = μ[mesh.cm[el].NodeList[i]]
  end

  fill!(At,zero(Float64))
  fill!(Bt,zero(Float64))
  fill!(Ct,zero(Float64))
  fill!(Et,zero(Float64))
  fill!(Ft,zero(Float64))
  fill!(Gt,zero(Float64))
  fill!(Ht,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)

    μg    = shapeEval(μNodes,phi)
    dμgdx = shapeEval(μNodes,dphidx)
    dμgdy = shapeEval(μNodes,dphidy)

    for ti=1:9
      # vector-valued test function
      for tj=1:9
        # integrand for A
        iA = μg*(dphidx[ti]*dphidx[tj] + dphidy[ti]*dphidy[tj])/jac -
             phi[ti]*(dμgdx*dphidx[tj] + dμgdy*dphidy[tj])/jac
        At[ti,tj] += iA*w[gpt]

        # integrand for B
        iB = 0.0
        Bt[ti,tj] += iB*w[gpt]

        # integrand for D -- transpose of B
        iD = 0.0
        Dt[ti,tj] += iD*w[gpt]

        # integrand for E
        iE = iA
        Et[ti,tj] += iE*w[gpt]
      end

      # scalar-valued test function
      for tj=1:4
        # integrand for C
        iC = -psi[tj]*dphidx[ti] # C is Nu x Np block matrix, so i and j are correct
        Ct[ti,tj] += iC*w[gpt]

        # integrand for F
        iF = -psi[tj]*dphidy[ti] # F is Nu x Np block matrix, so i and j are correct
        Ft[ti,tj] += iF*w[gpt]

        # integrand for G -- transpose of C
        iG = iC
        Gt[tj,ti] += iG*w[gpt]    # G is Np x Nu block matrix, so i and j are switched

        # integrand for H -- transpose of H
        iH = iF
        Ht[tj,ti] += iH*w[gpt]    # H is Np x Nu block matrix, so i and j are switched
      end
    end
  end

nothing
end

# Brinkman Constant parameter
function localBrinkmanConst2D!(mesh,el,xN,yN,w,s,t,Axt,Ayt,Bxt,Byt,
                               phi,dphidx,dphidy,dphids,dphidt,
                               psi,dpsidx,dpsidy,dpsids,dpsidt,
                               parameter...)
  μ = parameter[1]
  α = parameter[2]

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)

    for ti=1:9
      for tj=1:9
        integrandX = μ*dphidx[ti]*dphidx[tj]/jac +
                     μ*dphidy[ti]*dphidy[tj]/jac +
                     α*phi[ti]*phi[tj]*jac
        Axt[ti,tj] += integrandX*w[gpt]

        integrandY = μ*dphidx[ti]*dphidx[tj]/jac +
                     μ*dphidy[ti]*dphidy[tj]/jac +
                     α*phi[ti]*phi[tj]*jac
        Ayt[ti,tj] += integrandY*w[gpt]
      end

      for tk=1:4
        integrandX = psi[tk]*dphidx[ti]
        Bxt[tk,ti] -= integrandX*w[gpt]

        integrandY = psi[tk]*dphidy[ti]
        Byt[tk,ti] -= integrandY*w[gpt]
      end
    end
  end
end

# Brinkman Variable parameter
function localBrinkmanVar2D!(mesh,el,xN,yN,w,s,t,Axt,Ayt,A12t,A21t,Bxt,Byt,
                             phi,dphidx,dphidy,dphids,dphidt,
                             psi,dpsidx,dpsidy,dpsids,dpsidt,
                             parameter...)
  μ = parameter[1]
  α = parameter[2]

  μNodes = zeros(Float64,9)
  αNodes  = zeros(Float64,9)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  for i=1:9
    μNodes[i] = μ[mesh.cm[el].NodeList[i]]
    αNodes[i] = α[mesh.cm[el].NodeList[i]]
  end
  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(A12t,zero(Float64))
  fill!(A21t,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)

    μg = shapeEval(μNodes,phi)
    αg = shapeEval(αNodes,phi)

    for ti=1:9
      for tj=1:9
        integrandX  = 2*μg*dphidx[ti]*dphidx[tj]/jac +
                        μg*dphidy[ti]*dphidy[tj]/jac +
                        αg*phi[ti]*phi[tj]*jac
        Axt[ti,tj] += integrandX*w[gpt]

        integrand12  = μg*dphidx[tj]*dphidy[ti]/jac
        A12t[ti,tj] += integrand12*w[gpt]

        integrand21  = μg*dphidy[tj]*dphidx[ti]/jac
        A21t[ti,tj] += integrand21*w[gpt]

        integrandY   =   μg*dphidx[ti]*dphidx[tj]/jac +
                       2*μg*dphidy[ti]*dphidy[tj]/jac +
                         αg*phi[ti]*phi[tj]*jac
        Ayt[ti,tj]  += integrandY*w[gpt]
      end

      for tk=1:4
        integrandX = psi[tk]*dphidx[ti]
        Bxt[tk,ti] -= integrandX*w[gpt]

        integrandY = psi[tk]*dphidy[ti]
        Byt[tk,ti] -= integrandY*w[gpt]
      end
    end
  end
end

# Brinkman Multiphase 2D Constant Parameter
function localBrinkmanMPConst2D!(mesh,el,xN,yN,w,s,t,Axt,Ayt,Bxt,Byt,
                                  phi,dphidx,dphidy,dphids,dphidt,
                                  psi,dpsidx,dpsidy,dpsids,dpsidt,
                                  parameter...)
  α1 = parameter[1]
  α2 = parameter[2]
  α3 = parameter[3]

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)

  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)

    for ti=1:9
      for tj=1:9
        integrandX = α1*dphidx[ti]*dphidx[tj]/jac +
                     α1*dphidy[ti]*dphidy[tj]/jac +
                     α2*phi[ti]*phi[tj]*jac
        Axt[ti,tj] += integrandX*w[gpt]

        integrandY = α1*dphidx[ti]*dphidx[tj]/jac +
                     α1*dphidy[ti]*dphidy[tj]/jac +
                     α2*phi[ti]*phi[tj]*jac
        Ayt[ti,tj] += integrandY*w[gpt]
      end

      for tk=1:4
        integrandX = α3*psi[tk]*dphidx[ti]
        Bxt[tk,ti] -= integrandX*w[gpt]

        integrandY = α3*psi[tk]*dphidy[ti]
        Byt[tk,ti] -= integrandY*w[gpt]
      end
    end
  end
end

# Brinkman Multiphase -- wrong!! -- 2D Variable Parameter
function localBrinkmanMPVar2D!(mesh,el,xN,yN,w,s,t,Axt,Ayt,A12t,A21t,Bxt,Byt,
                               phi,dphidx,dphidy,dphids,dphidt,
                               psi,dpsidx,dpsidy,dpsids,dpsidt,
                               parameter...)
  α1 = parameter[1]
  α2 = parameter[2]
  α3 = parameter[3]

  α1Nodes = zeros(Float64,9)
  α2Nodes = zeros(Float64,9)
  α3Nodes = zeros(Float64,9)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  for i=1:9
    α1Nodes[i] = α1[mesh.cm[el].NodeList[i]]
    α2Nodes[i] = α2[mesh.cm[el].NodeList[i]]
    α3Nodes[i] = α3[mesh.cm[el].NodeList[i]]
  end

  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(A12t,zero(Float64))
  fill!(A21t,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

    for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)

    α1g = shapeEval(α1Nodes,phi)
    α2g = shapeEval(α2Nodes,phi)
    α3g = shapeEval(α3Nodes,phi)
    yg  = shapeEval(yN,phi)

    for ti=1:9
      for tj=1:9
        integrandX = 2*α1g*dphidx[ti]*dphidx[tj]/jac +
                       α1g*dphidy[ti]*dphidy[tj]/jac +
                       α2g*phi[ti]*phi[tj]*jac
        Axt[ti,tj] += integrandX*w[gpt]

        integrand12  = α1g*dphidx[tj]*dphidy[ti]/jac
        A12t[ti,tj] += integrand12*w[gpt]

        integrand21  = α1g*dphidy[tj]*dphidx[ti]/jac
        A21t[ti,tj] += integrand21*w[gpt]

        integrandY =   α1g*dphidx[ti]*dphidx[tj]/jac +
                     2*α1g*dphidy[ti]*dphidy[tj]/jac +
                       α2g*phi[ti]*phi[tj]*jac
        Ayt[ti,tj] += integrandY*w[gpt]
      end

      for tk=1:4
        integrandX = α3g*psi[tk]*dphidx[ti]
        Bxt[tk,ti] -= integrandX*w[gpt]

        integrandY = α3g*psi[tk]*dphidy[ti]
        Byt[tk,ti] -= integrandY*w[gpt]
      end
    end
  end
end
# Multiphase Brinkman (MPB) -- correct!! -- 2D Variable Parameter
function localMPBVar2D!(mesh,el,xN,yN,w,s,t,
                        At,Ct,Et,Ft,Gt,Ht,
                        phi,dphidx,dphidy,dphids,dphidt,
                        psi,dpsidx,dpsidy,dpsids,dpsidt,
                        parameter...)
  # make parameters easier to deal with
  α1 = parameter[1]
  α2 = parameter[2]
  α3 = parameter[3]

  # temporary arrays for variable-coefficients
  α1Nodes = zeros(Float64,9)
  α2Nodes = zeros(Float64,9)
  α3Nodes = zeros(Float64,9)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  for i=1:9
    α1Nodes[i] = α1[mesh.cm[el].NodeList[i]]
    α2Nodes[i] = α2[mesh.cm[el].NodeList[i]]
    α3Nodes[i] = α3[mesh.cm[el].NodeList[i]]
  end

  fill!(At,zero(Float64))
  fill!(Ct,zero(Float64))
  fill!(Et,zero(Float64))
  fill!(Ft,zero(Float64))
  fill!(Gt,zero(Float64))
  fill!(Ht,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)

    α1g = shapeEval(α1Nodes,phi)
    α2g = shapeEval(α2Nodes,phi)
    α3g = shapeEval(α3Nodes,phi)

    dα1gdx = shapeEval(α1Nodes,dphidx)
    dα1gdy = shapeEval(α1Nodes,dphidy)
    dα3gdx = shapeEval(α3Nodes,dphidx)
    dα3gdy = shapeEval(α3Nodes,dphidy)

    for ti=1:9
      # vector-valued test function
      for tj=1:9
        # integrand for A
        iA = α1g*(dphidx[ti]*dphidx[tj] + dphidy[ti]*dphidy[tj])/jac +
             phi[ti]*(dα1gdx*dphidx[tj] + dα1gdy*dphidy[tj])/jac +
             α2g*phi[ti]*phi[tj]*jac
        At[ti,tj] += iA*w[gpt]

        # integrand for E -- same as A
        iE = iA
        Et[ti,tj] += iE*w[gpt]    # A and E are same block matrices
      end

      # scalar-valued test function
      for tj=1:4
        # integrand for C
        iC = -psi[tj]*(α3g*dphidx[ti] + dα3gdx*phi[ti]) # C is Nu x Np block matrix, so i and j are correct
        Ct[ti,tj] += iC*w[gpt]

        # integrand for F
        iF = -psi[tj]*(α3g*dphidy[ti] + dα3gdy*phi[ti]) # F is Nu x Np block matrix, so i and j are correct
        Ft[ti,tj] += iF*w[gpt]

        # integrand for G
        iG = -psi[tj]*dphidx[ti]
        Gt[tj,ti] += iG*w[gpt]    # G is Np x Nu block matrix, so i and j are switched

        # integrand for H
        iH = -psi[tj]*dphidy[ti]
        Ht[tj,ti] += iH*w[gpt]    # H is Np x Nu block matrix, so i and j are switched
      end
    end
  end
end

# Constant diffusion Laplace Axisymmetric
function localLaplaceConstAS!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                         At,phi,dphidx,dphidy,dphids,dphidt,order,parameter...)
  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  # zero out At
  fill!(At,zero(Float64))

  for gpt=1:nGaussNodes
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
    jac = jacCalc(xN,yN,dphids,dphidt)
    yg = shapeEval(yN,phi)

    for ti = 1:nNodesPerElm, tj = 1:nNodesPerElm
      integrand = (dphidx[ti]*dphidx[tj] + dphidy[ti]*dphidy[tj])/jac
      At[ti,tj] += integrand*w[gpt]*yg
    end
  end

  nothing
end

# Constant diffusion Laplace Axisymmetric variable coefficient
function localLaplaceVarAS!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                            At,phi,dphidx,dphidy,dphids,dphidt,order,parameter)

  PeNodes = zeros(Float64,nGaussNodes)
  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  # zero out At
  fill!(At,zero(Float64))
  for i=1:nGaussNodes
    PeNodes[i] = parameter[mesh.cm[el].NodeList[i]]
  end
  for gpt=1:nGaussNodes
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
    jac = jacCalc(xN,yN,dphids,dphidt)
    Peg = shapeEval(PeNodes,phi)
    yg = shapeEval(yN,phi)

    for ti = 1:nNodesPerElm, tj = 1:nNodesPerElm
      integrand = (dphidx[ti]*dphidx[tj] + dphidy[ti]*dphidy[tj])/jac/Peg
      At[ti,tj] += integrand*w[gpt]*yg
    end
  end

  nothing
end

function localAdvDiffAS!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                         At,phi,dphidx,dphidy,dphids,dphidt,order,parameter)
  U = parameter.u
  V = parameter.v

  UNodes = zeros(Float64,nGaussNodes)
  VNodes = zeros(Float64,nGaussNodes)
  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(Float64))

  for i=1:nGaussNodes
    UNodes[i] = U[mesh.cm[el].NodeList[i]]
    VNodes[i] = V[mesh.cm[el].NodeList[i]]
  end

  for gpt=1:nGaussNodes
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
    jac = jacCalc(xN,yN,dphids,dphidt)
    ug  = shapeEval(UNodes,phi)
    vg  = shapeEval(VNodes,phi)
    yg  = shapeEval(yN,phi)

    for ti = 1:nNodesPerElm, tj = 1:nNodesPerElm
      integrand =  ug*phi[ti]*dphidx[tj] + vg*phi[ti]*dphidy[tj]
      At[ti,tj] += integrand*w[gpt]*yg
    end
  end

  nothing
end

function localAdvDiffAS_SUPGmod!(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
                                 At,phi,dphidx,dphidy,dphids,dphidt,order,parameter)
  U = parameter.u
  V = parameter.v

  UNodes = zeros(Float64,nGaussNodes)
  VNodes = zeros(Float64,nGaussNodes)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,nGaussNodes)
  fill!(At,zero(Float64))

  for i=1:nGaussNodes
    UNodes[i] = U[mesh.cm[el].NodeList[i]]
    VNodes[i] = V[mesh.cm[el].NodeList[i]]
  end

  # compute h (max element length)
  hk = sqrt(areaCalc(xN,yN))

  # compute element Peclet number
  wmag = 0.0
  elPe = 0.0
  if mesh.order == :Linear
    uavg = sum(UNodes)/4
    vavg = sum(VNodes)/4
    wmag = sqrt(uavg^2 + vavg^2)

    elPe = 0.5*wmag*hk/parameter.κ
  elseif mesh.order == :Quadratic
    wmag = sqrt(UNodes[9]^2 + VNodes[9]^2)

    elPe = 0.5*wmag*hk/parameter.κ
  else
    error("mesh.order not specified correctly")
  end

  # compute delta for SUPG
  δ = 0.5*hk/wmag*(1.0 - 1.0/elPe)

  if mesh.order == :Linear
    for gpt=1:nGaussNodes
      shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],order)
      derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,order)
      jac = jacCalc(xN,yN,dphids,dphidt)
      ug  = shapeEval(UNodes,phi)
      vg  = shapeEval(VNodes,phi)
      yg  = shapeEval(yN,phi)

      for ti = 1:nNodesPerElm, tj = 1:nNodesPerElm
        integrand = δ*(ug*dphidx[tj] + vg*dphidy[tj])*(ug*dphidx[ti] + vg*dphidy[ti]) -
                    parameter.κ*δ*(dphidy[tj]/yg)*(ug*dphidx[ti] + vg*dphidy[ti])         # might need to switch ti & tj
        At[ti,tj] += integrand*w[gpt]*yg
      end
    end

  elseif mesh.order == :Quadratic

  else
    error("mesh.order not specified correctly")
end

  nothing
end

function localStokesConstAS!(mesh,el,xN,yN,w,s,t,Axt,Ayt,Bxt,Byt,
                             phi,dphidx,dphidy,dphids,dphidt,
                             psi,dpsidx,dpsidy,dpsids,dpsidt,
                             parameter...)
  μ0 = parameter[1]
  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)

  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

  for gpt=1:9
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)
    yg  = shapeEval(yN,phi)

    for ti=1:9
      for tj=1:9
        integrandX = μ0*dphidx[ti]*dphidx[tj]*yg/jac +
                     μ0*dphidy[ti]*dphidy[tj]*yg/jac
        Axt[ti,tj] += integrandX*w[gpt]

        integrandY = μ0*dphidx[ti]*dphidx[tj]*yg/jac +
                     μ0*dphidy[ti]*dphidy[tj]*yg/jac +
                     μ0*phi[ti]*phi[tj]*jac/yg
        Ayt[ti,tj] += integrandY*w[gpt]
      end

      for tk=1:4
        integrandX = psi[tk]*dphidx[ti]*yg
        Bxt[tk,ti] -= integrandX*w[gpt]

        integrandY = psi[tk]*dphidy[ti]*yg + psi[tk]*phi[ti]*jac
        Byt[tk,ti] -= integrandY*w[gpt]
      end
    end
  end
end

function localStokesVarAS!(mesh,el,xN,yN,w,s,t,Axt,Ayt,A12t,A21t,Bxt,Byt,
                            phi,dphidx,dphidy,dphids,dphidt,
                            psi,dpsidx,dpsidy,dpsids,dpsidt,
                            parameter...)
  μ = parameter[1]
  μNodes = zeros(Float64,9)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  for i=1:9
    μNodes[i] = μ[mesh.cm[el].NodeList[i]]
  end

  nGpt = length(s)

  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(A12t,zero(Float64))
  fill!(A21t,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

  for gpt=1:nGpt
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)
    yg  = shapeEval(yN,phi)
    μg  = shapeEval(μNodes,phi)

    for ti=1:9
      for tj=1:9
        integrandX = 2.0*μg*dphidx[ti]*dphidx[tj]*yg/jac +
                         μg*dphidy[ti]*dphidy[tj]*yg/jac
        Axt[ti,tj] += integrandX*w[gpt]

        # NOTE: I'm pretty sure the A21t is wrong, but there's something later on in the assembly process
        #       that fixes this bug and so for the love of god I know it's wrong but please don't change it
        #       because it's fixed later on by another bug...will have to go through and check carefully
        #       at a later date.
        #
        # integrand21 should be = μg*dphidx[tj]*dphidy[ti]*yg/jac. Note indices
        integrand12 = μg*dphidx[tj]*dphidy[ti]*yg/jac
        A12t[ti,tj] += integrand12*w[gpt]

        integrand21 = μg*dphidx[ti]*dphidy[tj]*yg/jac
        A21t[ti,tj] += integrand21*w[gpt]

        integrandY =     μg*dphidx[ti]*dphidx[tj]*yg/jac +
                     2.0*μg*dphidy[ti]*dphidy[tj]*yg/jac +
                     2.0*μg*phi[ti]*phi[tj]*jac/yg
        Ayt[ti,tj] += integrandY*w[gpt]
      end

      for tk=1:4
        integrandX = psi[tk]*dphidx[ti]*yg
        Bxt[tk,ti] -= integrandX*w[gpt]

        integrandY = psi[tk]*dphidy[ti]*yg + psi[tk]*phi[ti]*jac
        Byt[tk,ti] -= integrandY*w[gpt]
      end
    end
  end
end
function localStokesVarASAlternate!(mesh,el,xN,yN,w,s,t,Axt,Ayt,Bxt,Byt,
                                    phi,dphidx,dphidy,dphids,dphidt,
                                    psi,dpsidx,dpsidy,dpsids,dpsidt,
                                    parameter...)
  μ = parameter[1]
  μNodes = zeros(Float64,9)

  nGpt = length(s)

  # generate local stiffness matrices
  getNodes!(xN,yN,mesh,el,9)
  for i=1:9
    μNodes[i] = μ[mesh.cm[el].NodeList[i]]
  end

  fill!(Axt,zero(Float64))
  fill!(Ayt,zero(Float64))
  fill!(Bxt,zero(Float64))
  fill!(Byt,zero(Float64))

  for gpt=1:nGpt
    shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],2)
    shape2D!(psi,dpsids,dpsidt,s[gpt],t[gpt],1)
    derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xN,yN,2)
    derivShape2D!(psi,dpsidx,dpsidy,dpsids,dpsidt,s[gpt],t[gpt],xN[1:4],yN[1:4],1)
    jac = jacCalc(xN,yN,dphids,dphidt)
    yg  = shapeEval(yN,phi)
    μg  = shapeEval(μNodes,phi)

    dμgdx = shapeEval(μNodes,dphidx)
    dμgdy = shapeEval(μNodes,dphidy)

    for ti=1:9
      for tj=1:9
        integrandX = μg*dphidx[ti]*dphidx[tj]*yg/jac +
                     μg*dphidy[ti]*dphidy[tj]*yg/jac -
                     dμgdx*phi[tj]*dphidx[ti]*yg/jac -
                     dμgdy*phi[tj]*dphidy[ti]*yg/jac
        Axt[ti,tj] += integrandX*w[gpt]

        integrandY = integrandX + μg*phi[ti]*phi[tj]/yg*jac
        Ayt[ti,tj] += integrandY*w[gpt]
      end

      for tk=1:4
        integrandX = psi[tk]*dphidx[ti]*yg
        Bxt[tk,ti] -= integrandX*w[gpt]

        integrandY = psi[tk]*dphidy[ti]*yg + psi[tk]*phi[ti]*jac
        Byt[tk,ti] -= integrandY*w[gpt]
      end
    end
  end
end

function getNodes(mesh,el,nGaussNodes)
  xNodes = Array{Float64}(nGaussNodes); yNodes = Array{Float64}(nGaussNodes)
  for i=1:nGaussNodes
    xNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].x
    yNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].y
  end
  return xNodes,yNodes
end

function getNodes!(xNodes,yNodes,mesh,el,nGaussNodes)
  for i=1:nGaussNodes
    xNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].x
    yNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].y
  end
end

# given 4 nodes, computes area
function areaCalc(xNodes,yNodes)
  area1 = xNodes[1]*(yNodes[2] - yNodes[3]) +
          xNodes[2]*(yNodes[3] - yNodes[1]) +
          xNodes[3]*(yNodes[1] - yNodes[2])
  area2 = xNodes[1]*(yNodes[4] - yNodes[3]) +
          xNodes[4]*(yNodes[3] - yNodes[1]) +
          xNodes[3]*(yNodes[1] - yNodes[4])

  return 0.5*(abs(area1) + abs(area2))
end


text:
include("ToricVarieties.jl")

######################
# 1: The Julia type for ToricDivisors
######################

struct toric_divisor
           GapToricDivisor::GapObj
end
export toric_divisor

######################
# 2: Generic constructors
######################

function create_divisor( coeffs::Vector{Int}, v::toric_variety )
    # create the divisor
    gap_coeffs = GapObj( coeffs )
    gap_divisor = GAP.Globals.CreateDivisor( gap_coeffs, v.GapToricVariety )

    # wrap and return
    return toric_divisor( gap_divisor )
end
export create_divisor

function divisor_of_character( character::Vector{Int}, v::toric_variety )
    # create the divisor
    gap_character = GapObj( character )
    gap_divisor = GAP.Globals.DivisorOfCharacter( gap_character, v.GapToricVariety )

    # wrap and return
    return toric_divisor( gap_divisor )
end
export divisor_of_character

function divisor_of_class( v::toric_variety, class::Vector{Int} )
    # create the divisor
    gap_class = GapObj( class )
    gap_divisor = GAP.Globals.DivisorOfGivenClass( v.GapToricVariety, gap_class )

    # wrap and return
    return toric_divisor( gap_divisor )
end
export divisor_of_class

######################
# 3: Properties
######################

function is_cartier( d::toric_divisor )
    return GAP.Globals.IsCartier( d.GapToricDivisor )::Bool
end
export is_cartier

function is_principal( d::toric_divisor )
    return GAP.Globals.IsPrincipal( d.GapToricDivisor )::Bool
end
export is_principal

function is_primedivisor( d::toric_divisor )
    return GAP.Globals.IsPrimedivisor( d.GapToricDivisor )::Bool
end
export is_primedivisor

function is_basepoint_free( d::toric_divisor )
    return GAP.Globals.IsBasepointFree( d.GapToricDivisor )::Bool
end
export is_basepoint_free

function is_ample( d::toric_divisor )
    return GAP.Globals.IsAmple( d.GapToricDivisor )::Bool
end
export is_ample

function is_very_ample( d::toric_divisor )
    if ! is_ample( d )
        @warn "Can (current) only tell for ample toric divisors if they are very ample."
        return "fail"
    end

    return GAP.Globals.IsVeryAmple( d.GapToricDivisor )::Bool
end
export is_very_ample

function is_numerically_effective( d::toric_divisor )
    return GAP.Globals.IsNumericallyEffective( d.GapToricDivisor )::Bool
end
export is_numerically_effective

max_stars_repo_name: fingolfin/JToric

max_issues_repo_name: fingolfin/JToric

max_forks_repo_name: fingolfin/JToric


text:
using DashBootstrapComponents, DashHtmlComponents

cards = dbc_card(
    dbc_cardbody([
        html_h5("Custom CSS", className="card-title"),
        html_p(
            "This card has inline styles applied controlling the width. " *
            "You could also apply the same styles with a custom CSS class.",
        ),
    ]),
    style=Dict("width" => "18rem"),
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# Autogenerated wrapper script for GeographicLib_jll for i686-w64-mingw32-cxx03
export CartConvert, ConicProj, GeoConvert, GeodSolve, GeodesicProj, GeoidEval, Gravity, MagneticField, Planimeter, RhumbSolve, TransverseMercatorProj, libGeographic

JLLWrappers.@generate_wrapper_header("GeographicLib")
JLLWrappers.@declare_executable_product(CartConvert)
JLLWrappers.@declare_executable_product(ConicProj)
JLLWrappers.@declare_executable_product(GeoConvert)
JLLWrappers.@declare_executable_product(GeodSolve)
JLLWrappers.@declare_executable_product(GeodesicProj)
JLLWrappers.@declare_executable_product(GeoidEval)
JLLWrappers.@declare_executable_product(Gravity)
JLLWrappers.@declare_executable_product(MagneticField)
JLLWrappers.@declare_executable_product(Planimeter)
JLLWrappers.@declare_executable_product(RhumbSolve)
JLLWrappers.@declare_executable_product(TransverseMercatorProj)
JLLWrappers.@declare_library_product(libGeographic, "libGeographic.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        CartConvert,
        "bin\\CartConvert.exe",
    )

    JLLWrappers.@init_executable_product(
        ConicProj,
        "bin\\ConicProj.exe",
    )

    JLLWrappers.@init_executable_product(
        GeoConvert,
        "bin\\GeoConvert.exe",
    )

    JLLWrappers.@init_executable_product(
        GeodSolve,
        "bin\\GeodSolve.exe",
    )

    JLLWrappers.@init_executable_product(
        GeodesicProj,
        "bin\\GeodesicProj.exe",
    )

    JLLWrappers.@init_executable_product(
        GeoidEval,
        "bin\\GeoidEval.exe",
    )

    JLLWrappers.@init_executable_product(
        Gravity,
        "bin\\Gravity.exe",
    )

    JLLWrappers.@init_executable_product(
        MagneticField,
        "bin\\MagneticField.exe",
    )

    JLLWrappers.@init_executable_product(
        Planimeter,
        "bin\\Planimeter.exe",
    )

    JLLWrappers.@init_executable_product(
        RhumbSolve,
        "bin\\RhumbSolve.exe",
    )

    JLLWrappers.@init_executable_product(
        TransverseMercatorProj,
        "bin\\TransverseMercatorProj.exe",
    )

    JLLWrappers.@init_library_product(
        libGeographic,
        "bin\\libGeographic.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()