text:
export isinvertible, contract

@doc Markdown.doc"""
    order(I::AlgAssAbsOrdIdl) -> AlgAssAbsOrd

> Returns the order containing $I$.
"""
@inline order(I::AlgAssAbsOrdIdl) = I.order

iszero(I::AlgAssAbsOrdIdl) = (I.iszero == 1)
isone(I::AlgAssAbsOrdIdl) = isone(abs(det(basis_matrix(I, copy = false))))

function Base.one(S::AlgAssAbsOrdIdlSet)
  O = order(S)
  M = identity_matrix(FlintZZ, degree(O))
  return ideal(O, M, :twosided, true)
end

function one(I::AlgAssAbsOrdIdl)
  return ideal(order(I), one(order(I)))
end

function Base.copy(I::AlgAssAbsOrdIdl)
  return I
end

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

function show(io::IO, a::AlgAssAbsOrdIdl)
  print(io, "Ideal of ")
  show(IOContext(io, :compact => true), order(a))
  println(io, " with basis matrix")
  print(io, a.basis_matrix)
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AlgAssAbsOrdIdl, dict::IdDict)
  b = typeof(a)(order(a))
  for i in fieldnames(typeof(a))
    if isdefined(a, i)
      if i != :order
        setfield!(b, i, Base.deepcopy_internal(getfield(a, i), dict))
      end
    end
  end
  return b
end

################################################################################
#
#  Ideal Set
#
################################################################################

function show(io::IO, a::AlgAssAbsOrdIdlSet)
  print(io, "Set of ideals of $(order(a))\n")
end

order(a::AlgAssAbsOrdIdlSet) = a.order

parent(I::AlgAssAbsOrdIdl) = AlgAssAbsOrdIdlSet(order(I))

function IdealSet(O::AlgAssAbsOrd)
   return AlgAssAbsOrdIdlSet(O)
end

elem_type(::Type{AlgAssAbsOrdIdlSet{S, T}}) where {S, T} = AlgAssAbsOrdIdl{S, T}

elem_type(::AlgAssAbsOrdIdlSet{S, T}) where {S, T} = AlgAssAbsOrdIdl{S, T}

parent_type(::Type{AlgAssAbsOrdIdl{S, T}}) where {S, T} = AlgAssAbsOrdIdlSet{S, T}

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

function assure_has_basis(a::AlgAssAbsOrdIdl{S, T}) where {S, T}
  if isdefined(a, :basis)
    return nothing
  end

  O = order(a)
  a.basis = Vector{AlgAssAbsOrdElem{S, T}}(undef, degree(O))
  for i = 1:degree(O)
    a.basis[i] = elem_from_mat_row(O, basis_matrix(a, copy = false), i)
  end
  return nothing
end

function assure_has_basis_mat_inv(a::AlgAssAbsOrdIdl)
  if isdefined(a, :basis_mat_inv)
    return nothing
  else
    a.basis_mat_inv = FakeFmpqMat(pseudo_inv(basis_matrix(a, copy = false)))
    return nothing
  end
end

################################################################################
#
#  Basis (matrices)
#
################################################################################

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

> Returns the basis of $a$.
"""
function basis(a::AlgAssAbsOrdIdl; copy::Bool = true)
  assure_has_basis(a)
  if copy
    return deepcopy(a.basis)
  else
    return a.basis
  end
end

@doc Markdown.doc"""
    basis_matrix(a::AlgAssAbsOrdIdl; copy::Bool = true) -> fmpz_mat

> Returns the basis matrix of $a$ with respect to the basis of the order.
"""
function basis_matrix(a::AlgAssAbsOrdIdl; copy::Bool = true)
  if copy
    return deepcopy(a.basis_matrix)
  else
    return a.basis_matrix
  end
end

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

> Returns the inverse of the basis matrix of $a$.
"""
function basis_mat_inv(a::AlgAssAbsOrdIdl; copy::Bool = true)
  assure_has_basis_mat_inv(a)
  if copy
    return deepcopy(a.basis_mat_inv)
  else
    return a.basis_mat_inv
  end
end

function det_of_basis_matrix(a::AlgAssAbsOrdIdl)
  d = fmpz(1)
  M = basis_matrix(a, copy = false)
  for i = 1:nrows(M)
    d = mul!(d, d, M[i, i])
  end
  return d
end

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

@doc Markdown.doc"""
    +(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

> Returns $a + b$.
"""
function +(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  if iszero(a)
    return deepcopy(b)
  elseif iszero(b)
    return deepcopy(a)
  end

  d = degree(order(a))
  M = vcat(basis_matrix(a), basis_matrix(b))
  deta = det_of_basis_matrix(a)
  detb = det_of_basis_matrix(b)
  M = hnf_modular_eldiv!(M, gcd(deta, detb), :lowerleft)
  M = sub(M, (d + 1):2*d, 1:d)
  return ideal(order(a), M, :nothing, true)
end

@doc Markdown.doc"""
    *(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

> Returns $a \cdot b$.
"""
function *(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  if iszero(a)
    return deepcopy(a)
  elseif iszero(b)
    return deepcopy(b)
  end

  d = degree(order(a))
  ba = basis(a, copy = false)
  bb = basis(b, copy = false)
  d2 = d^2

  M = zero_matrix(FlintZZ, d2, d)
  t = one(order(a))
  for i = 1:d
    for j = 1:d
      t = mul!(t, ba[i], bb[j])
      for k = 1:d
        M[(i - 1)*d + j, k] = coordinates(t; copy = false)[k]
      end
    end
  end

  deta = det_of_basis_matrix(a)
  detb = det_of_basis_matrix(b)
  M = hnf_modular_eldiv!(M, deta*detb, :lowerleft)
  M = sub(M, (d2 - d + 1):d2, 1:d)
  return ideal(order(a), M, :nothing, true)
end

@doc Markdown.doc"""
    ^(a::AlgAssAbsOrdIdl, e::Int) -> AlgAssAbsOrdIdl
    ^(a::AlgAssAbsOrdIdl, e::fmpz) -> AlgAssAbsOrdIdl

> Returns $a^e$.
"""
^(A::AlgAssAbsOrdIdl, e::Int) = Base.power_by_squaring(A, e)
^(A::AlgAssAbsOrdIdl, e::fmpz) = Base.power_by_squaring(A, BigInt(e))

@doc Markdown.doc"""
    intersect(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

> Returns $a \cap b$.
"""
function intersect(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  d = degree(order(a))
  H = vcat(basis_matrix(a), basis_matrix(b))
  K = left_kernel(H)[2]
  return ideal(order(a), _hnf(view(K, 1:d, 1:d)*basis_matrix(a, copy = false), :lowerleft), :nothing, true)
end

################################################################################
#
#  Ad hoc multiplication
#
################################################################################

*(x::Union{ Int, fmpz }, a::AlgAssAbsOrdIdl) = ideal(order(a), x*basis_matrix(a, copy = false), :nothing, true)
*(a::AlgAssAbsOrdIdl, x::Union{ Int, fmpz }) = ideal(order(a), basis_matrix(a, copy = false)*x, :nothing, true)

function *(x::AlgAssAbsOrdElem, a::AlgAssAbsOrdIdl)
  @assert parent(x) === order(a)
  @assert isleft_ideal(a)
  O = order(a)
  return ideal(O, x)*a
end

function *(a::AlgAssAbsOrdIdl, x::AlgAssAbsOrdElem)
  @assert parent(x) === order(a)
  @assert isright_ideal(a)
  O = order(a)
  return a*ideal(O, x)
end

################################################################################
#
#  Division
#
################################################################################

function divexact(a::AlgAssAbsOrdIdl, x::fmpz)
  O = order(a)
  M = FakeFmpqMat(basis_matrix(a, copy = false), deepcopy(x))
  if M.den != 1
    error("Ideal not divisible by x")
  end
  return ideal(O, M.num)
end

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

# Assumes that B_I[i] == elem_from_mat_row(O, M, i)
function defines_ideal(O::Union{ NfAbsOrd, AlgAssAbsOrd }, M::fmpz_mat, B_I::Vector{ <: Union{ NfAbsOrdElem, AlgAssAbsOrdElem } })
  if nrows(M) != degree(O) || ncols(M) != degree(O)
    return false, FakeFmpqMat()
  end
  local Minv
  try
    Minv = FakeFmpqMat(pseudo_inv(M))
  catch
    return false, FakeFmpqMat()
  end
  n = degree(O)
  B_O = basis(O, copy = false)

  N = zero_matrix(FlintZZ, n, n)
  t = O()
  for i = 1:n
    for j = 1:n
      t = mul!(t, B_O[i], B_I[j])
      for k = 1:n
        N[j, k] = deepcopy(coordinates(t, copy = false)[k])
      end
    end
    if !isone((N*Minv).den)
      return false, Minv
    end
    if iscommutative(O)
      continue
    end
    for j = 1:n
      t = mul!(t, B_I[i], B_O[j])
      for k = 1:n
        N[j, k] = deepcopy(coordinates(t, copy = false)[k])
      end
    end
    if !isone((N*Minv).den)
      return false, Minv
    end
  end
  return true, Minv
end

function defines_ideal(O::Union{ NfAbsOrd, AlgAssAbsOrd }, M::fmpz_mat)
  if nrows(M) != degree(O) || ncols(M) != degree(O)
    return false, FakeFmpqMat(M), elem_type(O)[]
  end
  B_I = Vector{elem_type(O)}(undef, degree(O))
  for i = 1:degree(O)
    B_I[i] = elem_from_mat_row(O, M, i)
  end
  b, Minv = defines_ideal(O, M, B_I)
  return b, Minv, B_I
end

function defines_ideal(O::Union{ NfAbsOrd, AlgAssAbsOrd }, B_I::Vector{ <: Union{ NfAbsOrdElem, AlgAssAbsOrdElem } })
  if length(B_I) != degree(O)
    return false, zero_matrix(O, degree(O), degree(O)), FakeFmpqMat()
  end
  M = zero_matrix(FlintZZ, degree(O), degree(O))
  for i = 1:degree(O)
    for j = 1:degree(O)
      M[i, j] = deepcopy(coordinates(B_I[i], copy = false)[j])
    end
  end
  b, Minv = defines_ideal(O, M, B_I)
  return b, M, Minv
end

@doc Markdown.doc"""
    ideal(O::AlgAssAbsOrd, M::fmpz_mat, side::Symbol = :nothing,
          M_in_hnf::Bool = false)
      -> AlgAssAbsOrdIdl

> Returns the ideal of $O$ with basis matrix $M$.
> If the ideal is known to be a right/left/twosided ideal of $O$, `side` may be
> set to `:right`/`:left`/`:twosided` respectively.
> If `M_in_hnf == true` it is assumed that $M$ is already in lower left HNF.
"""
function ideal(O::AlgAssAbsOrd{S, T}, M::fmpz_mat, side::Symbol = :nothing, M_in_hnf::Bool = false) where {S, T}
  !M_in_hnf ? M = _hnf(M, :lowerleft) : nothing
  a = AlgAssAbsOrdIdl{S, T}(O, M)
  _set_sidedness(a, side)
  return a
end

@doc Markdown.doc"""
    ideal(O::AlgAssAbsOrd, x::AlgAssAbsOrdElem) -> AlgAssAbsOrdIdl

> Returns the twosided principal ideal of $O$ generated by $x$.
"""
function ideal(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}) where {S, T}
  @assert parent(x) == O

  if iscommutative(O)
    a = ideal(O, x, :left)
    a.isright = 1
    return a
  end

  t1 = O()
  t2 = O()
  M = zero_matrix(FlintZZ, degree(O)^2, degree(O))
  b = basis(O, copy = false)
  for i = 1:degree(O)
    t1 = mul!(t1, b[i], x)
    ii = (i - 1)*degree(O)
    for j = 1:degree(O)
      t2 = mul!(t2, t1, b[j])
      elem_to_mat_row!(M, ii + j, t2)
    end
  end
  M = sub(_hnf(M, :lowerleft), nrows(M) - degree(O) + 1:nrows(M), 1:ncols(M))

  return ideal(O, M, :twosided, true)
end

@doc Markdown.doc"""
    ideal(O::AlgAssAbsOrd, x::AlgAssAbsOrdElem, action::Symbol) -> AlgAssAbsOrdIdl

> Returns the ideal $x \cdot O$, if `action == :left`, and $O \cdot x$, if
> `action == :right`.
"""
function ideal(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}, action::Symbol) where {S, T}
  @assert parent(x) == O
  M = representation_matrix(x, action)
  a = ideal(O, M)
  if iszero(x)
    a.iszero = 1
  end

  if action == :left
    a.isright = 1
  elseif action == :right
    a.isleft = 1
  end

  return a
end

@doc Markdown.doc"""
    *(O::AlgAssAbsOrd, x::AlgAssAbsOrdElem) -> AlgAssAbsOrdIdl
    *(O::AlgAssAbsOrd, x::Int) -> AlgAssAbsOrdIdl
    *(O::AlgAssAbsOrd, x::fmpz) -> AlgAssAbsOrdIdl
    *(x::AlgAssAbsOrdElem, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    *(x::Int, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    *(x::fmpz, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl

> Returns the ideal $O \cdot x$ or $x \cdot O$ respectively.
"""
*(O::AlgAssAbsOrd{S, T}, x::AlgAssAbsOrdElem{S, T}) where {S, T} = ideal(O, x, :right)
*(x::AlgAssAbsOrdElem{S, T}, O::AlgAssAbsOrd{S, T}) where {S, T} = ideal(O, x, :left)

*(O::AlgAssAbsOrd, x::Union{ Int, fmpz }) = O*O(x)
*(x::Union{ Int, fmpz }, O::AlgAssAbsOrd) = O(x)*O

@doc Markdown.doc"""
    ideal_from_z_gens(O::AlgAssAbsOrd, b::Vector{ <: AlgAssAbsOrdElem},
                      side::Symbol = :nothing)
      -> AlgAssAbsOrdIdl

> Returns the ideal of $O$ generated by the elements of `b` as a lattice over
> $\mathbb Z$.
"""
function ideal_from_z_gens(O::AlgAssAbsOrd, b::Vector{T}, side::Symbol = :nothing) where { T <: AlgAssAbsOrdElem }
  d = degree(O)
  @assert length(b) >= d

  M = zero_matrix(FlintZZ, length(b), d)
  for i = 1:length(b)
    elem_to_mat_row!(M, i, b[i])
  end
  M = _hnf(M, :lowerleft)
  if d < length(b)
    M = sub(M, (length(b) - d + 1):length(b), 1:d)
  end
  return ideal(O, M, side, true)
end

################################################################################
#
#  Extend/contract
#
################################################################################

@doc Markdown.doc"""
    extend(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl

> Returns $O \cdot a \cdot O$.
> It is assumed that the order of $a$ is contained in $O$.
"""
function extend(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd, action::Symbol = :twosided)
  # Assumes order(A) \subseteq O

  d = degree(O)
  X = map(O, basis(A, copy = false))
  Y = basis(O, copy = false)
  if action == :left || action == :twosided
    M = zero_matrix(FlintZZ, d^2, d)
    t = O()
    for i = 1:d
      ii = (i - 1)*d
      for j = 1:d
        t = Y[i]*X[j]
        elem_to_mat_row!(M, ii + j, t)
      end
    end
    M = sub(_hnf(M, :lowerleft), nrows(M) - d + 1: nrows(M), 1:d)
    if action == :left
      return ideal(O, M, :left, true)
    end
    for i = 1:d
      X[i] = elem_from_mat_row(O, M, i)
    end
  end
  if action == :right || action == :twosided
    M = zero_matrix(FlintZZ, d^2, d)
    t = O()
    for i = 1:d
      ii = (i - 1)*d
      for j = 1:d
        t = X[i]*Y[j]
        elem_to_mat_row!(M, ii + j, t)
      end
    end
    M = sub(_hnf(M, :lowerleft), nrows(M) - d + 1: nrows(M), 1:d)
    return ideal(O, M, action, true)
  end
  error("action should be :left, :right or :twosided")
end

@doc Markdown.doc"""
    *(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    *(O::AlgAssAbsOrd, a::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

> Returns $a \cdot O$ or $O \cdot a$ respectively.
> It is assumed that the order of $a$ is contained in $O$.
"""
*(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) = extend(A, O, :right)
*(O::AlgAssAbsOrd, A::AlgAssAbsOrdIdl) = extend(A, O, :left)

@doc Markdown.doc"""
    contract(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    intersect(a::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) -> AlgAssAbsOrdIdl
    intersect(O::AlgAssAbsOrd, a::AlgAssAbsOrdIdl) -> AlgAssAbsOrdIdl

> Returns $a \cap O$.
> It is assumed that the order of $a$ contains $O$.
"""
function contract(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd)
  # Assumes O \subseteq order(A)

  d = degree(O)
  M = basis_matrix(O, copy = false)*basis_mat_inv(order(A), copy = false)
  @assert M.den == 1
  H = vcat(basis_matrix(A), M.num)
  K = left_kernel(H)[2]
  M = sub(K, 1:d, 1:d)*basis_matrix(A, copy = false)
  M = M*basis_matrix(order(A), copy = false)*basis_mat_inv(O, copy = false)
  @assert M.den == 1
  M = _hnf(M.num, :lowerleft)
  return ideal(O, M, :nothing, true)
end

intersect(O::AlgAssAbsOrd, A::AlgAssAbsOrdIdl) = contract(A, O)
intersect(A::AlgAssAbsOrdIdl, O::AlgAssAbsOrd) = contract(A, O)

################################################################################
#
#  Move ideals to another algebra
#
################################################################################

# Let m: A -> B be an injection and I an ideal in B. Then this computes the
# ideal O_A\cap I, where O_A is the maximal order of A.
# OB is assumed to be a maximal order of B containing order(I).
function _as_ideal_of_smaller_algebra(m::AbsAlgAssMor, I::AlgAssAbsOrdIdl, OB::AlgAssAbsOrd)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(order(I)) == B
  OA = maximal_order(A)
  IB = extend(I, OB)
  # Transport OA to B
  M = zero_matrix(FlintZZ, dim(B), dim(B))
  for i = 1:dim(A)
    t = OB(m(elem_in_algebra(basis(OA, copy = false)[i])))
    elem_to_mat_row!(M, i, t)
  end
  # Compute the intersection of M and IB
  H = vcat(M, basis_matrix(IB))
  K = left_kernel(H)[2]
  N = view(K, 1:dim(B), 1:dim(B))*M
  # Map the basis to A
  basis_in_A = Vector{elem_type(OA)}(undef, dim(B))
  for i = 1:dim(B)
    t = elem_from_mat_row(OB, N, i)
    b, s = haspreimage(m, elem_in_algebra(t, copy = false))
    @assert b
    basis_in_A[i] = OA(s)
  end
  J = ideal_from_z_gens(OA, basis_in_A)
  return J
end

# Let m: A -> B be an injection and O an order in B. Then this computes the
# order O_A\cap O, where O_A is the maximal order of A.
# OB is assumed to be a maximal order of B containing O.
function _as_order_of_smaller_algebra(m::AbsAlgAssMor, O::AlgAssAbsOrd, OB::AlgAssAbsOrd)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(O) == B
  OA = maximal_order(A)
  # Transport OA to B
  M = zero_matrix(FlintZZ, dim(B), dim(B))
  for i = 1:dim(A)
    t = OB(m(elem_in_algebra(basis(OA, copy = false)[i])))
    elem_to_mat_row!(M, i, t)
  end
  # Compute the intersection of M and O (in OB)
  N = basis_matrix(O, copy = false)*basis_mat_inv(OB, copy = false)
  @assert N.den == 1
  H = vcat(M, N.num)
  K = left_kernel(H)[2]
  N = sub(K, 1:dim(B), 1:dim(B))*M
  # Map the basis to A
  basis_in_A = Vector{elem_type(A)}(undef, dim(B))
  for i = 1:dim(B)
    t = elem_from_mat_row(OB, N, i)
    b, s = haspreimage(m, elem_in_algebra(t, copy = false))
    @assert b
    basis_in_A[i] = s
  end
  M = zero_matrix(base_ring(A), length(basis_in_A), dim(A))
  for i = 1:length(basis_in_A)
    elem_to_mat_row!(M, i, basis_in_A[i])
  end
  M = FakeFmpqMat(M)
  MM = sub(hnf!(M, :lowerleft), (dim(B) - dim(A) + 1):dim(B), 1:dim(A))
  OO = Order(A, MM)
  return OO
end

# Let m: A -> B be an injection, I an ideal in A and O an order of B. Then this
# computes I*O.
# WE ASSUME THAT m(A) LIES IN THE CENTRE OF B!
function _as_ideal_of_larger_algebra(m::AbsAlgAssMor, I::AlgAssAbsOrdIdl, O::AlgAssAbsOrd)
  A = domain(m)
  B = codomain(m)
  @assert dim(A) <= dim(B)
  @assert algebra(order(I)) == A
  @assert algebra(O) == B

  M = zero_matrix(FlintZZ, dim(A)*dim(B), dim(B))
  X = Vector{elem_type(O)}(undef, dim(A))
  for i = 1:dim(A)
    t = elem_in_algebra(basis(I, copy = false)[i])
    X[i] = O(m(t))
  end
  for i = 1:dim(A)
    for j = 1:dim(B)
      t = X[i]*basis(O, copy = false)[j]
      elem_to_mat_row!(M, (i - 1)*dim(B) + j, t)
    end
  end
  M = sub(_hnf(M, :lowerleft), dim(B)*(dim(A) - 1) + 1:dim(A)*dim(B), 1:dim(B))
  return ideal(O, M, :nothing, true)
end

################################################################################
#
#  Inclusion of elements in ideals
#
################################################################################

@doc Markdown.doc"""
    in(x::AlgAssAbsOrdElem, a::AlgAssAbsOrdIdl) -> Bool

> Returns `true` if the order element $x$ is in $a$ and `false` otherwise.
"""
function in(x::AlgAssAbsOrdElem, I::AlgAssAbsOrdIdl)
  el = coordinates(x)
  y = matrix(FlintZZ, 1, length(el), el)
  M1, d = pseudo_inv(basis_matrix(I, copy = false))
  if FakeFmpqMat(y*M1, d).den == 1
    return true
  else
    return false
  end
end

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

@doc Markdown.doc"""
    ==(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl) -> Bool

> Returns `true` if $a$ and $b$ are equal and `false` otherwise.
"""
function ==(A::AlgAssAbsOrdIdl, B::AlgAssAbsOrdIdl)
  order(A) !== order(B) && return false
  return basis_matrix(A, copy = false) == basis_matrix(B, copy = false)
end

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

function Base.hash(A::AlgAssAbsOrdIdl, h::UInt)
  return Base.hash(basis_matrix(A, copy = false), h)
end

################################################################################
#
#  Idempotents
#
################################################################################

@doc Markdown.doc"""
    idempotents(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl)
      -> AlgAssAbsOrdElem, AlgAssAbsOrdElem

> Given coprime ideals $a$ and $b$, this function returns elements $x \in a$
> and $y \in b$ such that $x + y = 1$.
"""
function idempotents(a::AlgAssAbsOrdIdl, b::AlgAssAbsOrdIdl)

  !(order(a) === order(b)) && error("Parent mismatch")
  O = order(a)
  d = degree(O)

  V = zero_matrix(FlintZZ, 1 + 2*degree(O), 1 + 2*degree(O))
  V[1, 1] = 1
  u = coordinates(one(O))
  for i = 1:d
    V[1, i + 1] = u[i]
  end

  _copy_matrix_into_matrix(V, 2, 2, basis_matrix(a, copy = false))
  _copy_matrix_into_matrix(V, 2 + d, 2, basis_matrix(b, copy = false))
  for i = 2:d + 1
    V[i, d + i] = 1
  end

  H = hnf!(V)
  for i = 2:d + 1
    if H[1, i] != 0
      error("Ideals are not coprime")
    end
  end

  t = O()
  z = basis(a, copy = false)[1]*H[1, d + 2]
  for i = 2:d
    t = mul!(t, basis(a, copy = false)[i], H[1, d + 1 + i])
    z = add!(z, z, t)
  end

  @assert -z in a
  @assert one(O) + z in b
  return -z, one(O) + z
end

################################################################################
#
#  From algebra to number field
#
################################################################################

function _as_ideal_of_number_field(I::AlgAssAbsOrdIdl, m::AbsAlgAssToNfAbsMor)
  @assert algebra(order(I)) == domain(m)
  K = codomain(m)
  OK = maximal_order(K)

  b = Vector{elem_type(OK)}()
  for i = 1:dim(domain(m))
    push!(b, OK(m(elem_in_algebra(basis(I, copy = false)[i]))))
  end
  return ideal_from_z_gens(OK, b)
end

function _as_ideal_of_number_field(I::FacElem{<: AlgAssAbsOrdIdl, <: AlgAssAbsOrdIdlSet}, m::AbsAlgAssToNfAbsMor)
  @assert algebra(order(base_ring(parent(I)))) == domain(m)
  K = codomain(m)
  OK = maximal_order(K)

  if isempty(I)
    return FacElemMon(IdealSet(OK))()
  end

  bases = Vector{ideal_type(OK)}()
  exps = Vector{fmpz}()
  for (b, e) in I
    push!(bases, _as_ideal_of_number_field(b, m))
    push!(exps, e)
  end
  return FacElem(bases, exps)
end

function _as_ideal_of_algebra(I::NfAbsOrdIdl, i::Int, O::AlgAssAbsOrd, one_ideals::Vector{Vector{T}}) where { T <: AlgAssAbsOrdElem }
  A = algebra(O)
  fields_and_maps = as_number_fields(A)
  b = Vector{elem_type(O)}()
  for j = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[j]
    if j == i
      @assert nf(order(I)) == K
      append!(b, [ O(AtoK\K(bb)) for bb in basis(I, copy = false) ])
    else
      append!(b, one_ideals[j])
    end
  end
  return ideal_from_z_gens(O, b)
end

function _as_ideal_of_algebra(I::FacElem{NfOrdIdl, NfOrdIdlSet}, i::Int, O::AlgAssAbsOrd, one_ideals::Vector{Vector{T}}) where { T <: AlgAssAbsOrdElem }
  if isempty(I)
    return FacElemMon(IdealSet(O))()
  end

  base = Vector{ideal_type(O)}()
  exp = Vector{fmpz}()
  for (b, e) in I
    push!(base, _as_ideal_of_algebra(b, i, O, one_ideals))
    push!(exp, e)
  end
  return FacElem(base, exp)
end

# Returns an array of bases of the ideals O_i(1)*O_i lifted to O, where O_i
# are the maximal orders of the number fields in the decomposition of
# algebra(O).
function _lift_one_ideals(O::AlgAssAbsOrd)
  A = algebra(O)
  fields_and_maps = as_number_fields(A)

  one_ideals = Vector{Vector{elem_type(O)}}()
  for i = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[i]
    OK = maximal_order(K)
    a = OK(1)*OK
    b = Vector{elem_type(O)}()
    for i = 1:degree(K)
      push!(b, O(AtoK\K(basis(a, copy = false)[i])))
    end
    push!(one_ideals, b)
  end
  return one_ideals
end

################################################################################
#
#  Factorization
#
################################################################################

@doc Markdown.doc"""
    factor(I::AlgAssAbsOrdIdl) -> Dict{AlgAssAbsOrdIdl, Int}

> Returns the prime ideal factorization of $I$ as a dictionary.
"""
function factor(I::AlgAssAbsOrdIdl)
  O = order(I)
  A = algebra(O)
  fields_and_maps = as_number_fields(A)

  one_ideals = _lift_one_ideals(O)

  fac = Dict{ideal_type(O), Int}()
  for i = 1:length(fields_and_maps)
    K, AtoK = fields_and_maps[i]
    J = _as_ideal_of_number_field(I, AtoK)
    fac2 = factor(J)
    for (p, e) in fac2
      P = _as_ideal_of_algebra(p, i, O, one_ideals)
      fac[P] = e
    end
  end
  return fac
end

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

# If I is not coprime to the conductor of O in the maximal order, then this might
# not be an inverse.
function inv(a::AlgAssAbsOrdIdl)
  @assert !iszero(a)
  O = order(a)
  return colon(ideal(O, one(O)), a)
end

# Mostly the same as colon in NfOrd/Ideal/Ideal.jl
function colon(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  O = order(a)
  n = degree(O)
  B = basis(b)

  bmatinv = basis_mat_inv(a, copy = false)

  n = representation_matrix(B[1])*bmatinv
  m = numerator(n)
  d = denominator(n)
  for i in 2:length(B)
    n = representation_matrix(B[i])*bmatinv
    mm = n.num
    dd = n.den
    l = lcm(dd, d)
    if l == d && l == dd
      m = hcat(m, mm)
    elseif l == d
      m = hcat(m, div(d, dd)*mm)
    elseif l == dd
      m = hcat(div(dd, d)*m, mm)
      d = dd
    else
      m = hcat(m*div(l, d), mm*div(l, dd))
      d = l
    end
  end
  m = hnf(transpose(m))
  # n is upper right HNF
  m = transpose(sub(m, 1:degree(O), 1:degree(O)))
  b = inv(FakeFmpqMat(m, d))
  return frac_ideal_type(O)(O, ideal(O, b.num), b.den)
end

function isinvertible(a::AlgAssAbsOrdIdl)
  if iszero(a)
    return false, a
  end

  bmata = basis_matrix(a, copy = false)

  if any(iszero(bmata[i, i]) for i in 1:degree(order(a)))
    return false, a
  end

  if order(a).ismaximal == 1
    return true, inv(a)
  end

  b = inv(a)
  c = a*b
  return isone(c), b
end

################################################################################
#
#  Normal ideals
#
################################################################################

@doc Markdown.doc"""
    isnormal(a::AlgAssAbsOrdIdl) -> Bool

> Returns `true` if $a$ is a normal ideal and `false` otherwise.
"""
isnormal(a::AlgAssAbsOrdIdl) = ismaximal(left_order(a))

################################################################################
#
#  Locally free basis
#
################################################################################

@doc Markdown.doc"""
    locally_free_basis(a::AlgAssAbsOrdIdl, p::Union{ Int, fmpz }) -> AlgAssAbsOrdElem

> Returns an element $x$ of the order $O$ of $a$ such that $a_p = O_p x$
> where $p$ is a prime of `base_ring(O)`.
> See also `islocally_free`.
"""
function locally_free_basis(I::AlgAssAbsOrdIdl, p::Union{ Int, fmpz })
  b, x = islocally_free(I, p)
  if !b
    error("The ideal is not locally free at the prime")
  end
  return x
end

# See Bley, Wilson "Computations in relative algebraic K-groups", section 4.2
@doc Markdown.doc"""
    islocally_free(a::AlgAssAbsOrdIdl, p::Union{ Int, fmpz })
      -> Bool, AlgAssAbsOrdElem

> Returns a tuple `(true, x)` with an element $x$ of the order $O$ of $a$ such
> that $a_p = O_p x$ if $a$ is locally free at $p$, and `(false, 0)` otherwise.
> $p$ is a prime of `base_ring(O)`.
> See also `locally_free_basis`.
"""
function islocally_free(I::AlgAssAbsOrdIdl, p::Union{ Int, fmpz })
  O = order(I)
  OpO, toOpO = AlgAss(O, p*O, p)
  J = radical(OpO)
  OJ, toOJ = quo(OpO, J)
  decOJ = decompose(OJ)

  pI = p*I
  IpI, toIpI = AlgAss(I, p*I, p)
  gensJ = Vector{elem_type(IpI)}()
  for b in basis(J, copy = false)
    bb = toOpO\b
    for c in basis(I, copy = false)
      push!(gensJ, toIpI(bb*c))
    end
  end
  JinIpI = ideal_from_gens(IpI, gensJ)
  IJ, toIJ = quo(IpI, JinIpI)

  a = O()
  for i = 1:length(decOJ)
    A, AtoOJ = decOJ[i]
    B, AtoB, BtoA = _as_algebra_over_center(A)
    C, BtoC = _as_matrix_algebra(B)
    e = toOpO\(toOJ\(AtoOJ(BtoA(BtoC\C[1]))))
    basiseIJ = Vector{elem_type(IJ)}()
    for b in basis(I, copy = false)
      bb = toIJ(toIpI(e*b))
      if !iszero(bb)
        push!(basiseIJ, bb)
      end
    end

    # Construct an Fq-basis for e*IJ where Fq \cong centre(A)
    Z, ZtoA = center(A)
    basisZ = [ toOpO\(toOJ\(AtoOJ(ZtoA(Z[i])))) for i = 1:dim(Z) ]

    basiseIJoverZ = Vector{elem_type(O)}()
    M = zero_matrix(base_ring(IJ), dim(Z), dim(IJ))
    MM = zero_matrix(base_ring(IJ), 0, dim(IJ))
    r = 0
    for i = 1:length(basiseIJ)
      b = toIpI\(toIJ\basiseIJ[i])

      for j = 1:dim(Z)
        bb = toIJ(toIpI(basisZ[j]*b))
        elem_to_mat_row!(M, j, bb)
      end

      N = vcat(MM, M)
      s = rank(N)
      if s > r
        push!(basiseIJoverZ, b)
        MM = N
        r = s
      end
      if r == length(basiseIJ)
        break
      end
    end

    if length(basiseIJoverZ) != degree(C)
      # I is not locally free
      return false, O()
    end

    for i = 1:length(basiseIJoverZ)
      a += mod(toOpO\(toOJ\(AtoOJ(BtoA(BtoC\C[i]))))*basiseIJoverZ[i], pI)
    end
  end

  return true, mod(a, pI)
end

################################################################################
#
#  Ring of multipliers, left and right order
#
################################################################################

@doc Markdown.doc"""
    ring_of_multipliers(I::AlgAssAbsOrdIdl) -> AlgAssAbsOrd

> Given an ideal $I$, it returns the ring $(I : I)$.
"""
function ring_of_multipliers(I::AlgAssAbsOrdIdl, p::fmpz=fmpz(1), action::Symbol = :left)
  O = order(I)
  @hassert :AlgAssOrd 1 Hecke.check_associativity(algebra(O))
  @hassert :AlgAssOrd 1 Hecke.check_distributivity(algebra(O))
  @hassert :AlgAssOrd 1 defines_ideal(O, basis_matrix(I, copy = false))[1]
  bmatinv = basis_mat_inv(I, copy = false)
  if isdefined(I, :gens) && length(I.gens) < degree(O)
    B = I.gens
  else
    B = basis(I, copy = false)
  end
  m = zero_matrix(FlintZZ, degree(O)*length(B), degree(O))
  for i = 1:length(B)
    M = representation_matrix(B[i], action)
    mul!(M, M, bmatinv.num)
    if bmatinv.den == 1
      for s = 1:degree(O)
        for t = 1:degree(O)
          m[t + (i - 1)*degree(O), s] = M[s, t]
        end
      end
    else
      for s = 1:degree(O)
        for t = 1:degree(O)
          @hassert :AlgAssOrd 1 divisible(M[s, t], bmatinv.den)
          m[t + (i - 1)*degree(O), s] = divexact(M[s, t], bmatinv.den)
        end
      end
    end
  end
  #In the case of the p-radical, it is important to do this modulo p
  if p == 1
    m = hnf(m)
  else
    hnf_modular_eldiv!(m, p)
  end
  s = prod(m[i,i] for i=1:ncols(m))
  if s==1
    return O
  end
  # n is upper right HNF
  n = transpose(view(m, 1:degree(O), 1:degree(O)))
  b = FakeFmpqMat(pseudo_inv(n))
  mul!(b, b, basis_matrix(O, copy = false))
  @hassert :AlgAssOrd 1 defines_order(algebra(O), b)[1]
  O1 = Order(algebra(O), b)
  O1.disc = divexact(discriminant(O), s^2)
  return O1
end

@doc Markdown.doc"""
    left_order(a::AlgAssAbsOrdIdl) -> AlgAssAbsOrd

> Returns the order $\{ x \in A \mid xa \subseteq a\}$, where $A$ is the algebra
> containing $a$.
"""
left_order(a::AlgAssAbsOrdIdl) = ring_of_multipliers(a, fmpz(1), :right)

@doc Markdown.doc"""
    right_order(a::AlgAssAbsOrdIdl) -> AlgAssAbsOrd

> Returns the order $\{ x \in A \mid ax \subseteq a\}$, where $A$ is the algebra
> containing $a$.
"""
right_order(a::AlgAssAbsOrdIdl) = ring_of_multipliers(a, fmpz(1), :left)

################################################################################
#
#  p-radical
#
################################################################################

function pradical_meataxe(O::AlgAssAbsOrd, p::Int)

  A1, OtoA1  = quo(O, p*O, p)
  #@show dim(A1)
  @vtime :AlgAssOrd 1 lg = Hecke.gens(A1)
  #@show length(lg)
  lM = dense_matrix_type(base_ring(A1))[transpose(representation_matrix(lg[i])) for i=1:length(lg)]
  M = ModAlgAss(lM)
  ls = minimal_submodules(M)
  l = sum(nrows(x) for x in ls)
  M1 = zero_matrix(base_ring(A1), l, degree(O))
  i=1
  for x in ls
    for j=1:nrows(x)
      for k=1:degree(O)
        M1[i,k] = x[j,k]
      end
      i+=1
    end
  end
  r = rref!(M1)
  if r == degree(O)
    J = ideal(O, scalar_matrix(FlintZZ, degree(O), p))
    J.gens = elem_type(O)[O(p*one(O.algebra))]
    return J
  end
  M1 = view(M1, 1:r, 1:degree(O))
  dM = transpose(nullspace(M1)[2])
  gens=Vector{elem_type(O)}(undef, nrows(dM)+1)
  m = zero_matrix(FlintZZ, degree(O), degree(O))
  for i=1:nrows(dM)
    for j=1:ncols(dM)
      m[i,j] = lift(dM[i,j])
    end
    gens[i]= elem_from_mat_row(O,m,i)
  end
  hnf_modular_eldiv!(m, fmpz(p))
  gens[nrows(dM)+1]=O(p*one(algebra(O)))
  J=ideal(O,m)
  J.gens=gens
  return J

end

@doc Markdown.doc"""
    pradical(O::AlgAssAbsOrd, p::Int) -> AlgAssAbsOrdIdl

> Given an order $O$ and a prime $p$, it returns the radical of the ideal generated by $p$.
"""
function pradical(O::AlgAssAbsOrd, p::Int)

  #See the paper from Ronyai, Structure of finite algebra
  l = clog(fmpz(degree(O)),p)
  if l > 1
    return pradical_meataxe(O,p)
  end
  n = root(degree(O),2)
  F = GF(p, cached=false)

  #First step: kernel of the trace matrix mod p

  I = change_base_ring(n*trred_matrix(O), F)
  k, B = nullspace(I)
  # The columns of B give the coordinates of the elements in the order.
  if k==0
    J= ideal(O, scalar_matrix(FlintZZ, degree(O), p))
    J.gens = AlgAssAbsOrdElem[O(p*one(algebra(O)))]
    return J
  end
  if l==1
    #In this case, we can output I: it is the standard p-trace method.
    M=zero_matrix(FlintZZ, degree(O), degree(O))
    for i=1:ncols(B)
      for j=1:degree(O)
        M[i,j]=lift(B[j,i])
      end
    end
    M1 = hnf_modular_eldiv!(M, fmpz(p))
    res = ideal(O, view(M1,1:degree(O),1:degree(O)))
    B1 = lift(B')
    res.gens = Vector{elem_type(O)}(undef, k+1)
    for i=1:k
      res.gens[i] = elem_from_mat_row(O, B1, i)
    end
    res.gens[k+1] = O(p*one(algebra(O)))
    @hassert :AlgAssOrd 1 check_pradical(res,p)
    return res
  end

  Ide = transpose(lift(B))
  #Now, iterate: we need to find the kernel of tr((xy)^(p^i))/p^i mod p
  #on the subspace generated by I
  #Hard to believe, but this is linear on I!!!!
  for i = 1:l
    N = zero_matrix(F, degree(O), nrows(Ide))
    a = algebra(O)()
    for t = 1:nrows(Ide)
      elm = elem_from_mat_row(O,Ide,t)
      for s = 1:degree(O)
        mul!(a, elem_in_algebra(elm, copy = false), O.basis_alg[s])
        bbb = (a^(p^i))
        trel = tr(bbb)
        el = divexact(numerator(trel),p^i)
        N[s,t] = F(el)
      end
    end
    k, B2 = nullspace(N)
    if k == 0
      J = ideal(O, scalar_matrix(FlintZZ, degree(O), p))
      J.gens = AlgAssAbsOrdElem[O(p*one(algebra(O)))]
      return J
    end
    Ide = lift(transpose(B2))*Ide
  end
  gens = Vector{AlgAssAbsOrdElem}(undef, nrows(Ide)+1)
  for i in 1:nrows(Ide)
    gens[i] = elem_from_mat_row(O, Ide, i)
  end
  gens[nrows(Ide)+1]= O(p*one(algebra(O)))
  #Now, I have to give a basis for the ideal.
  m=zero_matrix(FlintZZ, degree(O), degree(O))
  for i=1:nrows(Ide)
    for j=1:ncols(Ide)
      m[i,j]=Ide[i,j]
    end
  end
  hnf_modular_eldiv!(m, fmpz(p))
  res = ideal(O, m)
  res.gens=gens
  @hassert :AlgAssOrd 1 check_pradical(res,p)
  return res

end

################################################################################
#
#  Some tests
#
################################################################################

function check_pradical(I::AlgAssAbsOrdIdl, p::Int)

  O= order(I)
  for i=1:degree(O)
    x=elem_from_mat_row(O,basis_matrix(I, copy = false), i)
    for j=1:degree(O)
      @assert divisible(numerator(tr(elem_in_algebra(x; copy = false)*O.basis_alg[j])), p)
    end
  end
  #=
  if p==2
    for i=1:O.dim
      x=elem_from_mat_row(O,I.basis_matrix, i)
      for j=1:O.dim
        for k=1:clog(fmpz(O.dim), p)
          @assert divisible(numerator(tr((x.elem_in_algebra*O.basis_alg[j])^(p^k))), p^(k+1))
        end
      end
    end
  end
  =#
  return true
end

################################################################################
#
#  Primes lying over a prime
#
################################################################################

# Returns the twosided maximal ideals of O containing I, where p*O \subseteq I.
# If strict_containment == true and I is already prime, we return an empty array.
function _maximal_ideals(O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, p::Union{Int, fmpz}; strict_containment::Bool = false)

  A1, OtoA1 = quo(O, I, p)
  @vtime :AlgAssOrd 1 lg = gens(A1)
  lM = dense_matrix_type(base_ring(A1))[representation_matrix(lg[i]) for i=1:length(lg)]
  append!(lM, dense_matrix_type(base_ring(A1))[representation_matrix(lg[i], :right) for i=1:length(lg)])
  M = ModAlgAss(lM)
  ls = maximal_submodules(M)
  if strict_containment && isone(length(ls)) && iszero(nrows(ls[1]))
    ls = typeof(ls)[]
  end
  poneO = O(p*one(algebra(O)))
  return ( _from_submodules_to_ideals(M, O, I, x, fmpz(p), poneO, A1, OtoA1) for x in ls )

end

function _from_submodules_to_ideals(M::ModAlgAss, O::AlgAssAbsOrd, I::AlgAssAbsOrdIdl, x::Union{Zmodn_mat, Generic.Mat{Generic.ResF{fmpz}}}, p::fmpz, poneO::AlgAssAbsOrdElem, A1::AlgAss, OtoA1::AbsOrdToAlgAssMor)
  @hassert :AlgAssOrd 1 begin r = rref(x)[1]; closure(x, M.action) == sub(rref(x)[2], 1:r, 1:ncols(x)) end
  m = zero_matrix(FlintZZ, nrows(x)+degree(O) , degree(O))
  gens = Vector{AlgAssAbsOrdElem}(undef, nrows(x))
  for i = 1:nrows(x)
    el = OtoA1\(elem_from_mat_row(A1, x, i))
    for j = 1:degree(O)
      m[i,j] = coordinates(el, copy = false)[j]
    end
    gens[i] = elem_from_mat_row(O, m, i)
  end
  for i=nrows(x)+1:nrows(m)
    for j=1:degree(O)
      m[i,j] = basis_matrix(I, copy = false)[i-nrows(x), j]
    end
  end
  hnf_modular_eldiv!(m, fmpz(p))
  J = ideal(O, view(m, 1:degree(O), 1:degree(O)))
  if isdefined(I, :gens)
    append!(gens, I.gens)
    J.gens = gens
  else
    append!(gens, basis(I, copy = false))
  end
  return J

end

@doc Markdown.doc"""
    prime_ideals_over(O::AlgAssAbsOrd, p::Union{ Int, fmpz })
      -> Vector{AlgAssAbsOrdIdl}

> Returns all prime ideal of $O$ lying over the prime number $p$.
"""
function prime_ideals_over(O::AlgAssAbsOrd, p::Union{ Int, fmpz })
  @vtime :AlgAssOrd 1 max_id = collect(_maximal_ideals(O, p*O, Int(p)))
  return max_id
end

@doc Markdown.doc"""
    prime_ideals_over(O::AlgAssAbsOrd, P::AlgAssAbsOrdIdl) -> Vector{AlgAssAbsOrdIdl}

> Given a prime ideal $P$ in an order contained in $O$, this function returns
> the prime ideals of $O$ lying over $P$.
"""
function prime_ideals_over(O::AlgAssAbsOrd, P::AlgAssAbsOrdIdl)
  O1 = order(P)
  if O1 == O
    return ideal_type(O)[P]
  end
  M = maximal_order(O)
  lp = prime_ideals_over(M, minimum(P))
  p_critical_primes = Vector{ideal_type(O)}()
  for Q in lp
    c = contract(Q, O1)
    if c == P
      c1 = contract(Q, O)
      if !(c1 in p_critical_primes)
        push!(p_critical_primes, c1)
      end
    end
  end
  return p_critical_primes
end

################################################################################
#
#  Minimum
#
################################################################################

function minimum(P::AlgAssAbsOrdIdl)
  N = norm(P, copy = false)
  f, p = ispower(N)
  @assert isprime(p)
  return p
end

################################################################################
#
#  Norm
#
################################################################################

function assure_has_norm(a::AlgAssAbsOrdIdl)
  if isdefined(a, :norm)
    return nothing
  end

  a.norm = abs(det(basis_matrix(a, copy = false)))
  return nothing
end

@doc Markdown.doc"""
    norm(a::AlgAssAbsordIdl; copy::Bool = true) -> fmpz

> Returns the norm of $a$.
"""
function norm(a::AlgAssAbsOrdIdl; copy::Bool = true)
  assure_has_norm(a)
  if copy
    return deepcopy(a.norm)
  else
    return a.norm
  end
end

function assure_has_normred(a::AlgAssAbsOrdIdl)
  if isdefined(a, :normred)
    return nothing
  end

  A = algebra(order(a))
  m = isqrt(dim(A))
  @assert m^2 == dim(A)
  N = norm(a, copy = false)
  b, n = ispower(N, m)
  @assert b "Cannot compute reduced norm. Maybe the algebra is not simple and central?"
  a.normred = n
  return nothing
end

@doc Markdown.doc"""
    normred(a::AlgAssAbsOrdIdl; copy::Bool = true) -> fmpz

> Returns the reduced norm of $a$.
> It is assumed that the algebra containing $a$ is simple and central.
"""
function normred(a::AlgAssAbsOrdIdl; copy::Bool = true)
  @assert dimension_of_center(algebra(order(a))) == 1
  @assert algebra(order(a)).issimple == 1
  assure_has_normred(a)
  if copy
    return deepcopy(a.normred)
  else
    return a.normred
  end
end


text:
module RayTrace # main model that calls on all the other modules

using Shapes, Vector3D

export trace! # nice, just one funnction to export

type Ray # the ray type
    pos::Vec # has a position, a location
    dir::Vec # and a direction, which has to be unitized, norm(dir) = 1
end

type Bull # this type is really just a container for the following fields. It can be thought of as a ray in the gradient lens.
    ray::Ray # the ray, note its origin is now the center of the lens
    N::Vec # the normal to the refractive surface
    dis::Float64 # the distance to the lens center
    n1::Float64 # the refractive index at the previous position
end

# functions:
# initiates the Bull type.
function Bull(ray::Ray,lens::Lens)
    ray.pos -= lens.c # set ray's origin on the lens' center
    return Bull(ray,Vec(0.,0.,0.),0.,0.)
end

# initiates the array of Rays. It starts them off from directly above the lens at distance L. Their directions are uniformly and randomly distributed across the lens.
function letBlight(lens::Lens,L::Float64)
    if isinf(L) # if the light source is infinitely far away
        # just ditribute the rays uniformly on the disk that the sphere makes:
        sqrtr = sqrt(rand()) # this is taken from http://mathworld.wolfram.com/DiskPointPicking.html
        θ = rand()*2π
        a = Vec([sqrtr*cos(θ),sqrtr*sin(θ),2lens.r]) # initiate the Ray with those starting points at distance 2*rad (could be any distance larger than rad really
        pos = lens.c + a # adjust in case the lens center is not at 0,0,0
        dir = Vec([0.,0.,-1.]) # just point them directly down
    else
        cosα = cos(asin(lens.r/L)) # theta is the angle between the optical axis and the ray that touches the surface of the sphere

        # z and α are randomly sampled from those ranges to generate the uniformly spaced points on the sphere
        z = rand()*(cosα - 1.) - cosα
        θ = rand()*2π

        # the x, y, and z that lie on the sphere's surface
        x = sqrt(1. - z*z)*cos(θ)
        y = sqrt(1. - z*z)*sin(θ)

        # this is the starting position of all the rays -- point source
        a = Vec([0.,0.,L])
        pos = lens.c + a
        dir = unitize(Vec([x,y,z]))
    end
    return Ray(pos,dir) # the ray type
end

# This function finds the intersection point of the light as it starts from its initial point. It's a basic vector-sphere intersection thing. It works. I copied it from the ray_sphere function in https://github.com/JuliaLang/julia/blob/master/test/perf/kernel/raytracer.jl
function intersect!(ray::Ray,s::Sphere)
    v = s.c - ray.pos
    b = dot(v, ray.dir)
    disc = b*b - dot(v, v) + s.r*s.r
    if disc >= 0.
        d = sqrt(disc)
        t2 = b + d
        if t2 >= 0.
            t1 = b - d
            ray.pos += (t1 > 0. ? t1 : t2)*ray.dir
        end
    end
end

# the refraction/reflection function, ray is ray, N is the normal to the refractive durface, n is the ratio between the refractive index before and afer the surface.
function refract!(ray::Ray,N::Vec,n::Float64)
    # this is the 'r' in page 74 in that pdf I just mentioned
    a = -dot(ray.dir,N)

    # this is the 'c' in page 74
    c = 1. - n*n*(1. - a*a)

    # and finally the 't' in page 74
    # added reflection as well... See page 54 (v is -ray.dir)
    t = c>0 ? n*ray.dir + (n*a - sqrt(c))*N : 2* dot(N,-ray.dir) *N + ray.dir

    # unitize the direction vector for the ray
    ray.dir = unitize(t)
end

# a step function. moves the ray forward one ε. Notice that the ray is now packed inside Bull, and is normalized to the lens center
function step!(bull::Bull,lens::Lens)
    bull.n1 = lens.ri(norm(bull.ray.pos)) # first I find the refractive index at the current position

    bull.ray.pos += lens.ε*bull.ray.dir # then I move the ray one ε forward with its original direction

    bull.dis = norm(bull.ray.pos) # get the distance to the lens center from the new location
    bull.N = bull.ray.pos/bull.dis # find the new location's normal
end

# bend the ray accordingly to the refraction/reflection.
function bend!(bull::Bull,lens::Lens)
    n2 = lens.ri(bull.dis) # new refractive index here

    # here I find out if the ray is on it's way in or way out. If its distance to the lens center is diminishing then the normal to the refractive surface is going outwards. But if the ray is leaving the lens the normal to the surface is pointing inwards.
    costheta = dot(bull.ray.dir,bull.ray.pos)

    out = -sign(costheta) # is it going out? If so make out equal to +1, otherwise make it -1

    bull.N = out*bull.ray.pos/bull.dis # this is the normal to the refractive surface

    refract!(bull.ray,bull.N,bull.n1/n2) # refraction

end

# main tracing function. p is a matrix where all the ray's locations are saved. lens is lens, retina is retina, L is the distance to the light source, n_medium you know, maxiterATION, diagnose is a Dict that contains some diagnostics.
function trace!(p::Array{Float64,2},lens::Lens,retina::Retina,L::Float64,n_medium::Float64,maxiter::Int,diagnose::Dict)
    ray = letBlight(lens,L) # initiate all the rays

    p[:,1] = myconvert(ray.pos) # populate with the start position

    intersect!(ray,lens) # find the intersection point of the ray with the lens and update ray's position

    p[:,2] = myconvert(ray.pos) # add said intersection to the position matrix

    refract!(ray,unitize(ray.pos - lens.c),n_medium/lens.ri(lens.r)) # refract at the surface of the lens due to the medium's refractive index (that might be different from that of the lens' periphery)
    bull = Bull(ray,lens) # make a bull instance
    for i = 3:maxiter-1 # go for maximum maxiterATIONS
        step!(bull,lens) # advance one ε
        if bull.dis > lens.r # is the ray exiting the lens? If so, intersect it with the lens (which is therefore closer)
            ray = bull.ray # extract the ray from bull
            ray.pos -= lens.ε*ray.dir-lens.c # correct for the mistake of exiting the lens AND the lens center
            intersect!(ray,lens) # intersect with the lens
            p[:,i] = myconvert(ray.pos) # add position to the position matrix
            diagnose[:Δ] = sqrt(sum((diagnose[:goal]-vec(p[:,i])).^2)) # find distance to the goal

            break # and break the loop, cause we're done with the lens
        end # if the ray didn't exit the lens:
        p[:,i] = myconvert(bull.ray.pos + lens.c) # save the new position in the positions matrix
        bend!(bull,lens) # and refract
    end
    refract!(ray,-unitize(ray.pos - lens.c),lens.ri(lens.r)/n_medium) # refract when exiting the lens and entering the surrounding medium
    intersect!(ray,retina) # intersection point with the retina
    i += 1 # just for the position matrix
    p[:,i] = myconvert(ray.pos) # add said position to the position matrix
    diagnose[:len] = i # and set the length of the ray refraction points in the positions matrix
end

end

ext: jl

max_stars_repo_name: JuliaTagBot/borat

max_issues_repo_name: JuliaTagBot/borat

max_forks_repo_name: JuliaTagBot/borat


text:
using Documenter, TestPkgLH

makedocs(
    modules = [TestPkgLH],
    format = :html,
    checkdocs = :exports,
    sitename = "TestPkgLH.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/hendri54/TestPkgLH.jl.git",
)


text:
include("mul.jl")
include("muladd.jl")
include("inv.jl")
include("lazymul.jl")
include("add.jl")

mulapplystyle(::TriangularLayout, ::AbstractStridedLayout) = LmulStyle()
mulapplystyle(::AbstractStridedLayout, ::TriangularLayout) = RmulStyle()


text:
#= note: to run this in parallel, from the system
prompt, and before starting julia, do
export JULIA_NUM_THREADS=X
where X is the desired number of threads. On
my 32 core system, 10 is good (cache contention issues)
Remember to re-set threads to 1 before using MPI
=#
using Statistics, Econometrics, LinearAlgebra, Plots
include("QIVmodel.jl")
function main()
LNW, X, Z = getdata()
n = size(LNW,1)
mcmcreps = 100000
burnin = 10000
# use st. errs. from ordinary QR as quide to setting tuning.
basetuning = [0.1, 0.05, 0.03, 0.005, 0.3, 0.2, 0.2]
# adjust this up or down to achieve decent acceptance rate
scale = 0.1 # higher for higher taus
tuning = scale*basetuning
EducEffect = zeros(9,3)
Constant = zeros(9,3)
for i = 1:9
    τ = round(i/10,digits=1)
    Σ = τ*(1.0-τ)*Z'Z/n
    Σinv = inv(Σ)
    θ = X\LNW  # OLS start values
    # to do ordinary QR via MCMC, set Z=X (just to verify that MCMC works!)
    lnL = θ -> likelihood(θ, LNW, X, Z, τ, Σinv)
    Proposal = θ -> proposal(θ, tuning)
    Prior = θ -> prior(θ)
    chain = mcmc(θ, mcmcreps, burnin, Prior, lnL, Proposal)
    V = cov(chain[:,1:end-1])
    d = dstats(chain)
    Constant[i,:] = d[1,[1,7,8]] # save posterior mean, 5% and 95% quantile
    EducEffect[i,:] = d[2,[1,7,8]] # save posterior mean, 5% and 95% quantile
    #scatter(chain[:,1],chain[:,2], show=true, reuse=false)
end

# plot results
# educ
τ = (1:9)/10
plot(reuse=false)
educ = EducEffect[:,1]
lb = EducEffect[:,2]
ub = EducEffect[:,3]
plot!(τ, [educ educ], fillrange=[lb ub], fillalpha=0.3, c=:green, legend=:none)
xticks!((1:9)/10)
#savefig("EducJUNK.png")
# constant
τ = (1:9)/10
plot(reuse=false)
c = Constant[:,1]
lb = Constant[:,2]
ub = Constant[:,3]
plot!(τ, [c c], fillrange=[lb ub], fillalpha=0.3, c=:green, legend=:none)
xticks!((1:9)/10)
#savefig("ConstantJUNK.png")
return
end
main()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
import Base.Iterators.Repeated

mutable struct Buffer{T} <: GPUArray{T, 1}
    id        ::GLuint
    len       ::Int
    buffertype::GLenum
    usage     ::GLenum
    context   ::Context
    function Buffer{T}(cpu_data_ptr::Ptr{T}, buff_length::Int, buffertype::GLenum, usage::GLenum) where T
        id = glGenBuffers()
	    # size of 0 can segfault it seems
	    buff_length = buff_length == 0 ? 1 : buff_length
        upload_data(id, cpu_data_ptr, buff_length * sizeof(T), buffertype, usage)
        obj = new{T}(id, buff_length, buffertype, usage, current_context())
        finalizer(free!, obj)
        obj
    end
end
#Function to deal with any Immutable type with Real as Subtype
function Buffer(
        buffer::DenseVector{T};
        buffertype::GLenum = GL_ARRAY_BUFFER, usage::GLenum = GL_STATIC_DRAW
    ) where T <: Real
    Buffer{T}(pointer(buffer), length(buffer), buffertype, usage)
end
function Buffer(
        buffer::DenseVector{T};
        buffertype::GLenum = GL_ARRAY_BUFFER, usage::GLenum = GL_STATIC_DRAW
    ) where T
    glasserteltype(T)
    Buffer{T}(pointer(buffer), length(buffer), buffertype, usage)
end
function Buffer(
        ::Type{T}, len::Int;
        buffertype::GLenum = GL_ARRAY_BUFFER, usage::GLenum = GL_STATIC_DRAW
    ) where T
    glasserteltype(T)
    Buffer{T}(Ptr{T}(C_NULL), len, buffertype, usage)
end

free!(x::Buffer) =
        context_command(x.context, () -> glDeleteBuffers(1, [x.id]))

function indexbuffer(
        buffer::Vector{T};
        usage::GLenum = GL_STATIC_DRAW
    ) where T
    glasserteltype(T)
    Buffer(buffer, buffertype = GL_ELEMENT_ARRAY_BUFFER, usage=usage)
end
indexbuffer(x::Buffer) = x.buffertype == GL_ELEMENT_ARRAY_BUFFER ? x : @error "Indexbuffer must be of enum GL_ELEMENT_ARRAY_BUFFER!"
#-------------------------------- END CONSTRUCTORS -------------------------------------#

Base.length(x::Buffer) = x.len

function upload_data(id::GLuint, cpu_data_ptr::Ptr, data_length::Int, buffertype::GLenum, usage::GLenum)
    glBindBuffer(buffertype, id)
    glBufferData(buffertype, data_length, cpu_data_ptr, usage)
    glBindBuffer(buffertype, 0)
end

function upload_data!(buf::Buffer{T}, new_data::Vector{T}) where {T}
	bind(buf)
	if length(new_data) == length(buf)
        glBufferSubData(buf.buffertype, 0, length(new_data) * sizeof(T), pointer(new_data))
    else
        glBufferData(buf.buffertype, length(new_data) * sizeof(T), pointer(new_data), buf.usage)
    end
    unbind(buf)
end

function reupload_data!(buf::Buffer{T}, new_data::Vector{T}, offset=0) where {T}
	bind(buf)
    glBufferSubData(buf.buffertype, offset, length(new_data) * sizeof(T), pointer(new_data))
    unbind(buf)
end

cardinality(::Buffer{T}) where {T} = cardinality(T)

bind(buffer::Buffer) = glBindBuffer(buffer.buffertype, buffer.id)
#used to reset buffer target
bind(buffer::Buffer, other_target) = glBindBuffer(buffer.buffertype, other_target)
unbind(buffer::Buffer) = glBindBuffer(buffer.buffertype, buffer.id)

Base.convert(::Type{Buffer}, x::Buffer)   = x
Base.convert(::Type{Buffer}, x::Array)    = Buffer(x)
Base.convert(::Type{Buffer}, x::Repeated) = convert(Buffer, x.xs.x)

#TODO: implement iteration
function Base.similar(x::Buffer{T}, buff_length::Int) where T
    Buffer{T}(Ptr{T}(C_NULL), buff_length, x.buffertype, x.usage)
end

# copy between two buffers
# could be a setindex! operation, with subarrays for buffers
function Base.unsafe_copyto!(a::Buffer{T}, readoffset::Int, b::Buffer{T}, writeoffset::Int, len::Int) where T
    multiplicator = sizeof(T)
    glBindBuffer(GL_COPY_READ_BUFFER, a.id)
    glBindBuffer(GL_COPY_WRITE_BUFFER, b.id)
    glCopyBufferSubData(
        GL_COPY_READ_BUFFER, GL_COPY_WRITE_BUFFER,
        multiplicator*(readoffset-1),
        multiplicator*(writeoffset-1),
        multiplicator*len
    )
    glBindBuffer(GL_COPY_READ_BUFFER, 0)
    glBindBuffer(GL_COPY_WRITE_BUFFER, 0)
    return nothing
end
#copy inside one buffer
function Base.unsafe_copyto!(buffer::Buffer{T}, readoffset::Int, writeoffset::Int, len::Int) where T
    len <= 0 && return nothing
    bind(buffer)
    ptr = Ptr{T}(glMapBuffer(buffer.buffertype, GL_READ_WRITE))
    for i=1:len+1
        unsafe_store!(ptr, unsafe_load(ptr, i+readoffset-1), i+writeoffset-1)
    end
    glUnmapBuffer(buffer.buffertype)
    bind(buffer,0)
    return nothing
end
function Base.unsafe_copyto!(a::Vector{T}, readoffset::Int, b::Buffer{T}, writeoffset::Int, len::Int) where T
    bind(b)
    ptr = Ptr{T}(glMapBuffer(b.buffertype, GL_WRITE_ONLY))
    for i=1:len
        unsafe_store!(ptr, a[i+readoffset-1], i+writeoffset-1)
    end
    glUnmapBuffer(b.buffertype)
    bind(b,0)
end
function Base.unsafe_copyto!(a::Buffer{T}, readoffset::Int, b::Vector{T}, writeoffset::Int, len::Int) where T
    bind(a)
    ptr = Ptr{T}(glMapBuffer(a.buffertype, GL_READ_ONLY))
    for i=1:len
        b[i+writeoffset-1] = unsafe_load(ptr, i+readoffset-2) #-2 => -1 to zero offset, -1 gl indexing starts at 0
    end
    glUnmapBuffer(a.buffertype)
    bind(a,0)
end

#--------------------------------- END BASEOVERLOADS-----------------------------------#

# GPUArray interface
function gpu_data(b::Buffer{T}) where T
    data = Vector{T}(length(b))
    bind(b)
    glGetBufferSubData(b.buffertype, 0, sizeof(data), data)
    bind(b, 0)
    data
end

# Resize buffer
function gpu_resize!(buffer::Buffer{T}, newdims::NTuple{1, Int}) where T
    #TODO make this safe!
    newlength = newdims[1]
    oldlen    = length(buffer)
    if oldlen > 0
        old_data = gpu_data(buffer)
    end
    bind(buffer)
    glBufferData(buffer.buffertype, newlength*sizeof(T), C_NULL, buffer.usage)
    bind(buffer, 0)
    buffer.length = newdims
    if oldlen>0
        max_len = min(length(old_data), newlength) #might also shrink
        buffer[1:max_len] = old_data[1:max_len]
    end
    #probably faster, but changes the buffer ID
    # newbuff     = similar(buffer, newdims...)
    # unsafe_copy!(buffer, 1, newbuff, 1, length(buffer))
    # buffer.id   = newbuff.id
    # buffer.length = newbuff.length
    nothing
end

function gpu_setindex!(b::Buffer{T}, value::Vector{T}, offset::Integer) where T
    multiplicator = sizeof(T)
    bind(b)
    glBufferSubData(b.buffertype, multiplicator*offset-1, sizeof(value), value)
    bind(b, 0)
end
function gpu_setindex!(b::Buffer{T}, value::Vector{T}, offset::UnitRange{Int}) where T
    multiplicator = sizeof(T)
    bind(b)
    glBufferSubData(b.buffertype, multiplicator*(first(offset)-1), sizeof(value), value)
    bind(b, 0)
    return nothing
end

function gpu_getindex(b::Buffer{T}, range::UnitRange) where T
    multiplicator = sizeof(T)
    offset        = first(range)-1
    value         = Vector{T}(undef, length(range))
    bind(b)
    glGetBufferSubData(b.buffertype, multiplicator*offset, sizeof(value), value)
    bind(b, 0)
    value
end


text:
"""
    LatinHypercube(gens = 1,
                   popsize = 100,
                   ntour = 2,
                   ptour = 0.8.,
                   interSampleWeight = 1.0,
                   ae_power = 2,
                   periodic_ae = false,
                   rng=Random.GLOBAL_RNG)

Instantiate grid-based hyperparameter tuning strategy using the
library [LatinHypercubeSampling.jl](https://github.com/MrUrq/LatinHypercubeSampling.jl).

An optimised Latin Hypercube sampling plan is created using a genetic
based optimization algorithm based on the inverse of the Audze-Eglais
function.  The optimization is run for `nGenerations` and creates `n`
models for evaluation, where `n` is specified by a corresponding `TunedModel`
instance, as in

    tuned_model = TunedModel(model=...,
                             tuning=LatinHypercube(...),
                             range=...,
                             measures=...,
                             n=...)

(See [`TunedModel`](@ref) for complete options.)

To use a periodic version of the Audze-Eglais function (to reduce
clustering along the boundaries) specify `periodic_ae = true`.

### Supported ranges:

A single one-dimensional range or vector of one-dimensioinal ranges
can be specified. Specifically, in `LatinHypercubeSampling` search,
the `range` field of a `TunedModel` instance can be:

- A single one-dimensional range - ie, `ParamRange` object - `r`, constructed
using the `range` method.

- Any vector of objects of the above form

Both `NumericRange`s and `NominalRange`s are supported, and
hyper-parameter values are sampled on a scale specified by the range
(eg, `r.scale = :log`).

"""
mutable struct LatinHypercube <: TuningStrategy
    gens::Int
    popsize::Int
    ntour::Int
    ptour::Number
    interSampleWeight::Number
    ae_power::Number
    periodic_ae::Bool
    rng::Random.AbstractRNG
end

function LatinHypercube(; gens = 1,
                        popsize = 100,
                        ntour = 2,
                        ptour = 0.8,
                        interSampleWeight = 1.0,
                        ae_power = 2,
                        periodic_ae = false,
                        rng=Random.GLOBAL_RNG)

    _rng = rng isa Integer ? Random.MersenneTwister(rng) : rng

    return LatinHypercube(gens, popsize, ntour,
                          ptour, interSampleWeight, ae_power,
                          periodic_ae, _rng)

end

function _create_bounds_and_dims_type(d,r)
    bounds = []
    dims_type = Array{LatinHypercubeSampling.LHCDimension}(undef,0)
    for i = 1:d
        if r[i] isa NumericRange
            if !(r[i].scale isa Symbol)
                throw(ArgumentError("Callable scale not supported."))
            end
            push!(dims_type,LatinHypercubeSampling.Continuous())
            if isfinite(r[i].lower) && isfinite(r[i].upper)
                push!(bounds,
                      Float64.([transform(MLJBase.Scale,
                                          MLJBase.scale(r[i].scale),
                                          r[i].lower),
                                transform(MLJBase.Scale,
                                          MLJBase.scale(r[i].scale),
                                          r[i].upper)]))
            elseif !isfinite(r[i].lower) && isfinite(r[i].upper)
                push!(bounds,
                      Float64.([transform(MLJBase.Scale,
                                          MLJBase.scale(r[i].scale),
                                          r[i].upper - 2*r[i].unit),
                                transform(MLJBase.Scale,
                                          MLJBase.scale(r[i].scale),
                                          r[i].upper)]))
            elseif isfinite(r[i].lower) && !isfinite(r[i].upper)
                push!(bounds, Float64.([transform(MLJBase.Scale,
                                                 MLJBase.scale(r[i].scale),
                                                 r[i].lower),
                                       transform(MLJBase.Scale,
                                                 MLJBase.scale(r[i].scale),
                                                 r[i].lower + 2*r[i].unit)]))
            else
                push!(bounds, Float64.([transform(MLJBase.Scale,
                                                 MLJBase.scale(r[i].scale),
                                                 r[i].origin - r[i].unit),
                                       transform(MLJBase.Scale,
                                                 MLJBase.scale(r[i].scale),
                                                 r[i].origin + r[i].unit)]))
            end
        else
            push!(dims_type,
                  LatinHypercubeSampling.Categorical(length(r[i].values),
                                                     1.0))
            push!(bounds,Float64.([1,length(r[i].values)]))
        end
    end
    return Tuple.(bounds), dims_type
end

# Hyper-param values generated by LatinHypercubeSampling library need
# further rescalings, and possibly rounding, depending on range type

# nominal ranges:
_transform(r::NominalRange, x) = r.values[round(Int, x)]

# numeric ranges:
_transform_noround(r, x) =
    inverse_transform(MLJBase.Scale, MLJBase.scale(r.scale), x)
_transform(r::NumericRange, x) = _transform_noround(r, x)
_transform(r::NumericRange{T}, x) where T<:Integer =
    round(T, _transform_noround(r, x))

function setup(tuning::LatinHypercube, model, range, n, verbosity)
    ranges = range isa AbstractVector ? range : [range, ]
    d = length(ranges)
    bounds, dims_type = _create_bounds_and_dims_type(d, ranges)
    plan, _ = LatinHypercubeSampling.LHCoptim(n, d, tuning.gens,
                    rng = tuning.rng,
                    popsize = tuning.popsize,
                    ntour = tuning.ntour,
                    ptour = tuning.ptour,
                    dims = dims_type,
                    interSampleWeight = tuning.interSampleWeight,
                    periodic_ae = tuning.periodic_ae,
                    ae_power = tuning.ae_power)
    scaled_plan = LatinHypercubeSampling.scaleLHC(plan, bounds)
    rescaled_plan = map(tuple(1:d...)) do k
        broadcast(x -> _transform(ranges[k], x), scaled_plan[:,k])
    end

    fields = map(r -> r.field, ranges)
    parameter_scales = scale.(ranges)
    models = makeLatinHypercube(model, fields, rescaled_plan)
    state = (models=models,
             fields=fields,
             parameter_scales=parameter_scales)
    return state
end

function MLJTuning.models(tuning::LatinHypercube,
                          model,
                          history,
                          state,
                          n_remaining,
                          verbosity)
     return state.models[_length(history) + 1:end], state
end

tuning_report(tuning::LatinHypercube, history, state) =
    (plotting = plotting_report(state.fields, state.parameter_scales, history),)

function makeLatinHypercube(prototype::Model, fields, plan)
    N = length(first(plan))
    map(1:N) do i
        clone = deepcopy(prototype)
        for k in eachindex(fields)
            recursive_setproperty!(clone,fields[k], plan[k][i])
        end
        clone
    end
end

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
@inline _select(::StaticInt{M}, ::StaticInt{M}) where {M} = StaticInt{M}()
@noinline _select(::StaticInt{M}, ::StaticInt{N}) where {M,N} = throw("$M ≠ $N")
@inline _select(::StaticInt{M}, _) where {M} = StaticInt{M}()
@inline _select(_, ::StaticInt{M}) where {M} = StaticInt{M}()
@inline _select(x, _) = x

"""
Checks sizes for compatibility, and preserves the static size information if
given a mix of static and dynamic sizes.
"""
@inline function matmul_sizes(C,A,B)
    MC, NC = size(C)
    MA, KA = size(A)
    KB, NB = size(B)
    @assert ((MC == MA) & (KA == KB) & (NC == NB)) "Size mismatch."
    (_select(MA, MC), _select(KA, KB), _select(NB, NC))
end

function unsafe_copyto_turbo!(pB, pA, M, N)
    LoopVectorization.@turbo for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
        pB[m,n] = pA[m,n]
    end
end

function default_stridedpointer_quote(::Type{T}, N, Ot) where {T}
  C = 1
  B = 0
  R = Expr(:tuple)
  o = Expr(:tuple)
  xt = Expr(:tuple)
  st = StaticInt(sizeof(T))
  gf = GlobalRef(Core,:getfield)
  for n ∈ 1:N
    push!(xt.args, Expr(:call, :*, :st, Expr(:call, gf, :x, n, false)))
    push!(R.args, n)
    push!(o.args, Expr(:call, Ot))
  end
  quote
    $(Expr(:meta,:inline))
    st = $st
    si = StrideIndex{$N,$R,$C}($xt, $o)
    stridedpointer(ptr, si, StaticInt{$B}())
  end
end

@generated function default_stridedpointer(ptr::Ptr{T}, x::X) where {T, N, X <: Tuple{Vararg{Integer,N}}}
    default_stridedpointer_quote(T, N, :One)
end
@generated function default_zerobased_stridedpointer(ptr::Ptr{T}, x::X) where {T, N, X <: Tuple{Vararg{Integer,N}}}
    default_stridedpointer_quote(T, N, :Zero)
end

@generated function splitfirstdim(sp::StridedPointer{T,N,C,B,R}, ::StaticInt{M}) where {T,N,C,B,R,M}
  gf = GlobalRef(Core,:getfield)
  fx = M*sizeof(T)
  xt = Expr(:tuple, :x0, Expr(:call, *, :x0, StaticInt(M*sizeof(T))))
  ot = Expr(:tuple, Zero(), Zero())
  R₁ = (R[1])::Int
  Rn = Expr(:tuple, R₁, R₁+1)
  for n ∈ 2:N
    push!(xt.args, Expr(:call, gf, :x, n, false))
    push!(ot.args, Zero())
    Rₙ = (R[n])::Int
    push!(Rn.args, Core.ifelse(Rₙ > R₁, Rₙ + 1, Rₙ))
  end
  Cn = Core.ifelse(C > 1, C+1, C)
  Bn = Core.ifelse(B > 1, B+1, B)
  quote
    $(Expr(:meta,:inline))
    x = strides(sp)
    x0 = $gf(x, 1, false)
    si = StrideIndex{$(N+1),$Rn,$Cn}($xt, $ot)
    stridedpointer($gf(sp,:p), si, StaticInt{$Bn}())
  end
end
@generated function droplastdim(sp::StridedPointer{T,N,C,B,R}) where {T,N,C,B,R}
  Cn = Core.ifelse(C == N, -1, C)
  Bn = Core.ifelse(B == N, -1, B)
  xt = Expr(:tuple)
  ot = Expr(:tuple)
  rt = Expr(:tuple)
  gf = GlobalRef(Core,:getfield)
  for n ∈ 1:N-1
    push!(xt.args, Expr(:call, gf, :x, n, false))
    push!(ot.args, Expr(:call, gf, :o, n, false))
    push!(rt.args, R[n])
  end
  quote
    $(Expr(:meta,:inline))
    x = strides(sp)
    o = offsets(sp)
    si = StrideIndex{$(N-1),$rt,$Cn}($xt, $ot)
    stridedpointer($gf(sp,:p), si, StaticInt{$Bn}())
  end
end


text:
"""
Authors : Adrien BLASSIAU Corentin JUVIGNY
TP RODD - 2019-2020
"""

@doc """
Module contenant 3 fonctions de génération d'instances du problème G_ULS.
  Il contient 3 fonctions
  - `generate(T, M, R, Emax, dmin, dmax)` génère une instance.
  - `generate(T, M, R, Emax, dmin, dmax, filename)` génère une instance et écrit le tout dans le fichier filename.
  - `generate(filename) renvoie l'instance décrite dans le fichier filename (généré avec la fonction précédente)

Utilisez using .Generator pour utiliser le module et son contenu dans votre code.

Utiliser ?NOM dans l'interface en ligne de commande Julia pour avoir de la documentation sur NOM.
Par exemple ?generate affiche la documetation sur les fonctions generate
""" ->
module Generator

  export generate

  include("./problem.jl")
  import .G_ULS: Instance
  using Dates

  @doc """
  Prérequis : installation des paquets suivants :
  (using Pkg ; Pkg.add("Distributions"))
  (using Pkg ; Pkg.add("Dates"))

  Fonction de génération d'instance.
  """ ->
  function generate(T::Int, M::Int, R::Int, Emax::Int, f::Array{Int}, e::Array{Int}, h::Int, p::Array{Int}, d::Array{Int})

    inst = Instance()

    if (length(f) != M) || (length(e) != M) || (length(p) != M)
      println("Impossible de générer une instance, f ou e de maucvaise taille")
      return inst
    end

    # Paramètres de base
    inst.T = T
    inst.M = M
    inst.R = R

    # Emax
    for i in 1:inst.T
        push!(inst.Emax, Emax)
    end

    # f
    for i in 1:inst.T
      push!(inst.f, f)
    end

    # e
    for i in 1:inst.T
      push!(inst.e, e)
    end

    # h
    for i in 1:inst.T
      push!(inst.h, h)
    end

    # p
    for i in 1:inst.T
      push!(inst.p, [])
      for m in 1:inst.M
        push!(inst.p[i], p[m])
      end
    end

    # d
    for i in 1:inst.T
      push!(inst.d, d[i])
    end

    return inst
  end

  @doc """
  Prérequis : installation des paquets suivants :
  (using Pkg ; Pkg.add("Distributions"))
  (using Pkg ; Pkg.add("Dates"))

  Utilise la fonction `generate(w, h, n, pw, ph, bounds_ω, bounds_ma)` pour construction une instance puis écriture de cette instance dans le fichier dont le chemin est `filename`.
  """ ->
  function generate(T::Int, M::Int, R::Int, Emax::Int, f::Array{Int}, e::Array{Int}, h::Int, p::Array{Int}, d::Array{Int}, filename::String)

    inst = generate(T, M, R, Emax, f, e, h, p, d)

    datetimenow = Dates.now()
    date = Dates.format(datetimenow, "e, dd u yyyy HH:MM:SS")

    open(filename, "w") do file
      write(file,
        "/**********************************************\n* Author: Adrien BLASSIAU Corentin JUVIGNY\n* Creation Date: $date\n*********************************************/")
      write(file, "\n\n")

      write(file, "T=$(inst.T);\n")
      write(file, "M=$(inst.M);\n")
      write(file, "R=$(inst.R);\n")
      write(file, "Emax=$(inst.Emax);\n")

      write(file, "f=[")
      for i in 1:inst.T-1
        write(file, "$(inst.f[i])")
        write(file, ",")
      end
      write(file, "$(inst.f[T])];\n")

      write(file, "e=[")
      for i in 1:inst.T-1
        write(file, "$(inst.e[i])")
        write(file, ",")
      end
      write(file, "$(inst.e[T])];\n")

      write(file, "h=$(inst.h);\n")

      write(file, "p=[")
      for i in 1:inst.T-1
        write(file, "$(inst.p[i])")
        write(file, ",")
      end
      write(file, "$(inst.p[T])];\n")

      write(file, "d=$(inst.d);\n")

    end
    return inst
  end

  @doc """
  generate(filename::String)

  Lecture fichier de données.

  Cette fonction lit un fichier dont le chemin est `filename` et renvoie un objet de type Instance contenant les informations écrites dans le fichier.
  """ ->
  function generate(filename::String)

    open(filename) do file
      lines = readlines(file)

      inst = Instance()

      line = lines[6]
      inst.T = parse(Int64, SubString(line,3,length(line)-1))

      line = lines[7]
      inst.M = parse(Int64, SubString(line,3,length(line)-1))

      line = lines[8]
      inst.R = parse(Int64, SubString(line,3,length(line)-1))

      line = lines[9]
      inst.Emax = Meta.parse(SubString(line,6,length(line)-1)) |> eval

      line = lines[10]
      inst.f = Meta.parse(SubString(line,3,length(line)-1)) |> eval

      line = lines[11]
      inst.e = Meta.parse(SubString(line,3,length(line)-1)) |> eval

      line = lines[12]
      inst.h = Meta.parse(SubString(line,3,length(line)-1)) |> eval

      line = lines[13]
      inst.p = Meta.parse(SubString(line,3,length(line)-1)) |> eval

      line = lines[14]
      inst.d = Meta.parse(SubString(line,3,length(line)-1)) |> eval

      return inst
    end

  end
end

max_stars_repo_name: AdrienBlassiau/RODD

max_issues_repo_name: AdrienBlassiau/RODD

max_forks_repo_name: AdrienBlassiau/RODD


text:
import Revise
import GLRenderer as GL
import ThreeDP3 as T
import Images as I
import MiniGSG as S
import Rotations as R
import PoseComposition: Pose
using Gen
import Serialization
import MeshCatViz

import MeshCatViz
MeshCatViz.setup_visualizer()

YCB_DIR = "/home/nishadg/mcs/ThreeDVision.jl/data/ycbv2"
world_scaling_factor = 100.0
id_to_cloud, id_to_shift, id_to_box  = T.load_ycbv_models_adjusted(YCB_DIR, world_scaling_factor);
all_ids = sort(collect(keys(id_to_cloud)));

# +
# IDX =6
IDX =301
# IDX =629
IDX = 539
IDX = 1228
IDX = 6
IDX = 913
IDX = 1898
IDX = 539

@show T.get_ycb_scene_frame_id_from_idx(YCB_DIR,IDX)

# +
gt_poses, ids, gt_rgb_image, gt_depth_image, cam_pose, camera = T.load_ycbv_scene_adjusted(
    YCB_DIR, IDX, world_scaling_factor, id_to_shift
);
dense_poses, dense_ids = T.load_ycbv_dense_fusion_predictions_adjusted(YCB_DIR, IDX, world_scaling_factor, id_to_shift);
dense_poses_reordered = [dense_poses[findfirst(dense_ids .== id)]   for id in ids];
inf_poses = Serialization.deserialize("../assets/results/$(lpad(IDX,4,"0")).poses").poses_after
inf_poses = T.adjust_poses(inf_poses, ids, world_scaling_factor, id_to_shift)

img = I.colorview(I.Gray, gt_depth_image ./ maximum(gt_depth_image))
I.colorview(I.RGB, permutedims(Float64.(gt_rgb_image)./255.0, (3,1,2)))

# +
resolution = 0.5

renderer_color = GL.setup_renderer(camera, GL.RGBMode())
obj_paths = T.load_ycb_model_obj_file_paths(YCB_DIR)
for id in all_ids
    v,n,f,_ = renderer_color.gl_instance.load_obj_parameters(obj_paths[id])
    v = v * world_scaling_factor
    v .-= id_to_shift[id]'
    GL.load_object!(renderer_color, v, n, f)
end
# colors = I.distinguishable_colors(length(ids), I.colorant"green")
colors = [
    I.colorant"yellow", I.colorant"cyan", I.colorant"lightgreen",
    I.colorant"red", I.colorant"purple", I.colorant"orange"
]

# +
focus_idx = 2

rgb_image_dense, depth_image = GL.gl_render(
    renderer_color, [ids[focus_idx]], [dense_poses_reordered[focus_idx]], [I.colorant"cyan"], T.IDENTITY_POSE)
rgb_image_dense= I.colorview(I.RGBA, permutedims(rgb_image_dense,(3,1,2)))

new_gt_rgb_image = cat(Float64.(gt_rgb_image) ./ 255.0, ones(size(gt_rgb_image)[1:2]),dims=3)
new_gt_rgb_image = I.colorview(I.RGBA, permutedims(new_gt_rgb_image,(3,1,2)))
@show size(new_gt_rgb_image)
mask = (rgb_image_dense .!= I.RGBA(1.0, 1.0, 1.0, 1.0))
new_gt_rgb_image[mask] .= alpha.*new_gt_rgb_image[mask] +(1-alpha).*rgb_image_dense[mask]
new_gt_rgb_image
# +

rgb_image_3dp3, depth_image = GL.gl_render(
    renderer_color, [ids[focus_idx]], [inf_poses[focus_idx]], [I.colorant"limegreen"], T.IDENTITY_POSE)
rgb_image_3dp3= I.colorview(I.RGBA, permutedims(rgb_image_3dp3,(3,1,2)))

new_gt_rgb_image = cat(Float64.(gt_rgb_image) ./ 255.0, ones(size(gt_rgb_image)[1:2]),dims=3)
new_gt_rgb_image = I.colorview(I.RGBA, permutedims(new_gt_rgb_image,(3,1,2)))
@show size(new_gt_rgb_image)
@show size(rgb_image_3dp3)
alpha = 0.4
mask = (rgb_image_3dp3 .!= I.RGBA(1.0, 1.0, 1.0, 1.0))
new_gt_rgb_image[mask] .= alpha.*new_gt_rgb_image[mask] +(1-alpha).*rgb_image_3dp3[mask]
new_gt_rgb_image
# -

# +

@gen function slam_model(T, room_bounds, wall_colors, cov)
    prev_data = nothing
    for t in 1:T
        if t==1
            pose ~ pose_uniform(room_bounds[1,:],room_bounds[2,:])
        else
            pose ~ pose_gaussian(prev_pose, [1.0 0.0;0.0 1.0] * 0.1, deg2rad(20.0))
        end
        prev_pose = pose
    end

    rgb, depth = GL.gl_render(renderer, wall_colors, pose)
    sense_depth ~ mvnormal(depth[1,:], cov)
    sense_rgb ~ color_distribution(rgb, 0.999)
end

# +

corners = get_corners(get_depth(tr_gt,t))
likely_poses = []
for c in corners
    for c2 in gt_corners
        p = c2 * inv(c)
        push!(likely_poses, p)
    end
end

for _ in 1:10
    PF.pf_move_accept!(pf_state, Gen.metropolis_hastings, (pose_mixture_proposal,
            (likely_poses, t, [1.0 0.0;0.0 1.0] * 0.5, deg2rad(5.0))), 10);
end

# +

@gen function generative_model()
    floor = Plane([0., 0., 0.], [0., 0., 1.])
    room_height = @trace(uniform(2.5, 3.0), :room_height)
    ceiling = Plane([0., 0., room_height], [0., 0., -1.])
    objects = [floor, ceiling]

    camera_z = @trace(uniform(0.2, 2.0), :z)
    camera_location = [0., 0., camera_z]

    camera_pitch = @trace(uniform(pi/2 - pi/4, pi/2 + pi/4), :pitch)
    camera_roll = @trace(uniform(-pi/4, pi/4), :roll)
    camera_rotation = make_rotation_matrix(camera_pitch, 0., camera_roll)

    depths = render(objects, camera_location, camera_rotation)

    noise = 0.1
    @trace(realsense_sensor(depths, noise), :observation)
end

# initialize trace with first observation
frame = get_frame(depth_camera)
constraints = Gen.choicemap((:observation, frame))
trace, = Gen.generate(generative_model, (), constraints)

while true

    trace, = Gen.mh(trace, Gen.select(:pitch, :roll, :z, :room_height))
    trace, = Gen.mh(trace, random_walk, (pi/64, :pitch))
    trace, = Gen.mh(trace, random_walk, (pi/64, :roll))
    trace, = Gen.mh(trace, random_walk, (0.05, :z))
    trace, = Gen.mh(trace, random_walk, (0.05, :room_height))
    trace, = Gen.mh(trace, Gen.select(:noise))

    # update trace with new observation
    frame = get_frame(depth_camera)
    constraints = Gen.choicemap((:observation, frame))
    trace, = Gen.update(trace, (), (), constraints)
end


text:
# using DifferentialEquations
using LinearAlgebra
using Statistics
using Distributions
using Plots
using Random

function Lorenz63(u, θ)
    σ,β,ρ = θ
    x,y,z = u

    du = [
        σ * (y-x);
        x*(ρ-z) - y;
        x*y - β*z
    ]
    return du
end

function JLorenz63(u, θ)
    σ,β,ρ = θ
    x,y,z = u

    J = [
        -σ σ 0;
        (ρ-z) -1 -x;
        y x -β
    ]

    return J
end

function RK4(f, u, dt, θ)
    k1 = f(u,θ)
    k2 = f(u+k1*dt/2, θ)
    k3 = f(u+k2*dt/2, θ)
    k4 = f(u+k3*dt, θ)

    u_new = u + (dt/6)*(k1 + 2k2 + 2k3 + k4)
    return u_new
end

function JRK4(f, Jf, u, dt, θ)
    n = size(u,1)
    k1 = f(u, θ)
    k2 = f(u+k1*dt/2, θ)
    k3 = f(u+k2*dt/2, θ)
    #k4 = rhs(state+k3*dt,*args)

    dk1 = Jf(u, θ)
    dk2 = Jf(u+k1*dt/2, θ) * (I(n) + dk1*dt/2)
    dk3 = Jf(u+k2*dt/2, θ) * (I(n) + dk2*dt/2)
    dk4 = Jf(u+k3*dt, θ) * (I(n) + dk3*dt)

    return I(n) + (dt/6)*(dk1 + 2dk2 + 2dk3 + dk4)
end

function EKF(ub, w, ObsOp, JObsOp, R, B)
    n = size(ub, 1)

    Dh = JObsOp(ub)

    D = Dh*B*Dh' + R
    K = B * Dh' * inv(D)

    ua = ub + K * (w - ObsOp(ub))
    P = (I(n) - K*Dh)*B

    return ua, P
end

function h(u)
    w = u
    return w
end

function Dh(u)
    n = size(u, 1)
    D = I(n)
    return D
end

θ = [10.0, (8/3), 28.0]
dt = 0.01
tm = 10
nt = Int(tm/dt)
t = 0:dt:tm

u0True = [1.0; 1.0; 1.0]
Random.seed!(1)
σ_m = 0.15
R = σ_m^2 .* I(3)

dt_m = 0.2
tm_m = 2
nt_m = Int(tm_m/dt_m)

ind_m = (Int(dt_m/dt)+1):Int(dt_m/dt):(Int(tm_m/dt)+1)
t_m = t[ind_m]

uTrue = zeros(3,nt+1)
uTrue[:,1] = u0True
km = 1
w = zeros(3, nt_m)
for k ∈ 1:nt
    uTrue[:, k+1] = RK4(Lorenz63, uTrue[:,k], dt, θ)

    if (km <= nt_m) && (k+1==ind_m[km])
        w[:,km] = h(uTrue[:,k+1]) + rand(Normal(0,σ_m), 3)
        km = km+1
    end
end

u0b = [2.0; 3.0; 4.0]
σ_b = 0.1
B = σ_b^2 .* I(3)
Q = 0.0 * I(3)

ub = zeros(3, nt+1)
ub[:,1] = u0b

ua = zeros(3, nt+1)
ua[:,1] = u0b
km = 1

for k ∈ 1:nt
    ub[:, k+1] = RK4(Lorenz63, ub[:,k], dt, θ)
    ua[:, k+1] = RK4(Lorenz63, ua[:,k], dt, θ)

    DM = JRK4(Lorenz63, JLorenz63, ua[:,k], dt, θ)

    B = DM * B * DM' +  Q

    if (km <= nt_m) && (k+1 == ind_m[km])
        ua[:, k+1], B = EKF(ua[:, k+1], w[:, km], h, Dh, R, B)
        km = km+1
    end

end

var_names = ["x(t)", "y(t)", "z(t)"]
plots = []
for i ∈ 1:3
    p = plot(t, uTrue[i, :], color=:black, alpha=0.75, linewidth=2, label="truth")
    plot!(p, t_m, w[i, :], seriestype=:scatter, color="light green", label="observation" )
    plot!(p, t, ub[i, :], linestyle=:dashdot, color="purple", label="forecast")
    plot!(p, t, ua[i, :], color="royal blue", alpha=1, linewidth = 1, label="analysis")
    xlabel!(p, "t")
    ylabel!(p, var_names[i])
    push!(plots, p)
end

plot(plots..., layout=(3,1))


text:
using MappedArrays
using Test

@test isempty(detect_ambiguities(MappedArrays, Base, Core))

using FixedPointNumbers, OffsetArrays, ColorTypes

@testset "ReadonlyMappedArray" begin
    a = [1,4,9,16]
    s = view(a', 1:1, [1,2,4])

    b = @inferred(mappedarray(sqrt, a))
    @test parent(b) === a
    @test eltype(b) == Float64
    @test @inferred(getindex(b, 1)) == 1
    @test b[2] == 2
    @test b[3] == 3
    @test b[4] == 4
    @test_throws ErrorException b[3] = 0
    @test isa(eachindex(b), AbstractUnitRange)
    b = mappedarray(sqrt, a')
    @test isa(eachindex(b), AbstractUnitRange)
    b = mappedarray(sqrt, s)
    @test isa(eachindex(b), CartesianIndices)
    c = Base.unaliascopy(b)
    @test c == b
    @test c !== b
end

@testset "MappedArray" begin
    intsym = Int == Int64 ? :Int64 : :Int32
    a = [1,4,9,16]
    s = view(a', 1:1, [1,2,4])
    c = @inferred(mappedarray(sqrt, x->x*x, a))
    @test parent(c) === a
    @test @inferred(getindex(c, 1)) == 1
    @test c[2] == 2
    @test c[3] == 3
    @test c[4] == 4
    c[3] = 2
    @test a[3] == 4
    @test_throws InexactError(intsym, Int, 2.2^2) c[3] = 2.2  # because the backing array is Array{Int}
    @test isa(eachindex(c), AbstractUnitRange)
    b = @inferred(mappedarray(sqrt, a'))
    @test isa(eachindex(b), AbstractUnitRange)
    c = @inferred(mappedarray(sqrt, x->x*x, s))
    @test isa(eachindex(c), CartesianIndices)

    d = Base.unaliascopy(c)
    @test c == d
    @test c !== d

    sb = similar(b)
    @test isa(sb, Array{Float64})
    @test size(sb) == size(b)

    a = [0x01 0x03; 0x02 0x04]
    b = @inferred(mappedarray(y->N0f8(y,0), x->x.i, a))
    for i = 1:4
        @test b[i] == N0f8(i/255)
    end
    b[2,1] = 10/255
    @test a[2,1] == 0x0a
end

@testset "of_eltype" begin
    a = [0.1 0.3; 0.2 0.4]
    b = @inferred(of_eltype(N0f8, a))
    @test b[1,1] === N0f8(0.1)
    b = @inferred(of_eltype(zero(N0f8), a))
    @test b[1,1] === N0f8(0.1)
    b[2,1] = N0f8(0.5)
    @test a[2,1] == N0f8(0.5)
    @test !(b === a)
    b = @inferred(of_eltype(Float64, a))
    @test b === a
    b = @inferred(of_eltype(0.0, a))
    @test b === a
end

@testset "OffsetArrays" begin
    a = OffsetArray(randn(5), -2:2)
    aabs = mappedarray(abs, a)
    @test axes(aabs) == (-2:2,)
    for i = -2:2
        @test aabs[i] == abs(a[i])
    end
end

@testset "No zero(::T)" begin
    astr = @inferred(mappedarray(length, ["abc", "onetwothree"]))
    @test eltype(astr) == Int
    @test astr == [3, 11]
    a = @inferred(mappedarray(x->x+0.5, Int[]))
    @test eltype(a) == Float64

    # typestable string
    astr = @inferred(mappedarray(uppercase, ["abc", "def"]))
    @test eltype(astr) == String
    @test astr == ["ABC","DEF"]
end

@testset "ReadOnlyMultiMappedArray" begin
    a = reshape(1:6, 2, 3)
#    @test @inferred(axes(a)) == (Base.OneTo(2), Base.OneTo(3))
    b = fill(10.0f0, 2, 3)
    M = @inferred(mappedarray(+, a, b))
    @test @inferred(eltype(M)) == Float32
    @test @inferred(IndexStyle(M)) == IndexLinear()
    @test @inferred(IndexStyle(typeof(M))) == IndexLinear()
    @test @inferred(size(M)) === size(a)
    @test @inferred(axes(M)) === axes(a)
    @test M == a + b
    @test @inferred(M[1]) === 11.0f0
    @test @inferred(M[CartesianIndex(1, 1)]) === 11.0f0

    d = Base.unaliascopy(b)
    @test d == b
    @test d !== b

    c = view(reshape(1:9, 3, 3), 1:2, :)
    M = @inferred(mappedarray(+, c, b))
    @test @inferred(eltype(M)) == Float32
    @test @inferred(IndexStyle(M)) == IndexCartesian()
    @test @inferred(IndexStyle(typeof(M))) == IndexCartesian()
    @test @inferred(axes(M)) === axes(c)
    @test M == c + b
    @test @inferred(M[1]) === 11.0f0
    @test @inferred(M[CartesianIndex(1, 1)]) === 11.0f0
end

@testset "MultiMappedArray" begin
    intsym = Int == Int64 ? :Int64 : :Int32
    a = [0.1 0.2; 0.3 0.4]
    b = N0f8[0.6 0.5; 0.4 0.3]
    c = [0 1; 0 1]
    f = RGB{N0f8}
    finv = c->(red(c), green(c), blue(c))
    M = @inferred(mappedarray(f, finv, a, b, c))
    @test @inferred(eltype(M)) == RGB{N0f8}
    @test @inferred(IndexStyle(M)) == IndexLinear()
    @test @inferred(IndexStyle(typeof(M))) == IndexLinear()
    @test @inferred(size(M)) === size(a)
    @test @inferred(axes(M)) === axes(a)
    @test M[1,1] === RGB{N0f8}(0.1, 0.6, 0)
    @test M[2,1] === RGB{N0f8}(0.3, 0.4, 0)
    @test M[1,2] === RGB{N0f8}(0.2, 0.5, 1)
    @test M[2,2] === RGB{N0f8}(0.4, 0.3, 1)
    M[1,2] = RGB(0.25, 0.35, 0)
    @test M[1,2] === RGB{N0f8}(0.25, 0.35, 0)
    @test a[1,2] == N0f8(0.25)
    @test b[1,2] == N0f8(0.35)
    @test c[1,2] == 0
    @test_throws InexactError(intsym, Int, N0f8(0.45)) M[1,2] = RGB(0.25, 0.35, 0.45)
    R = reinterpret(N0f8, M)
    @test R == N0f8[0.1 0.25; 0.6 0.35; 0 0; 0.3 0.4; 0.4 0.3; 0 1]
    R[2,1] = 0.8
    @test b[1,1] === N0f8(0.8)

    a = view(reshape(0.1:0.1:0.6, 3, 2), 1:2, 1:2)
    M = @inferred(mappedarray(f, finv, a, b, c))
    @test @inferred(eltype(M)) == RGB{N0f8}
    @test @inferred(IndexStyle(M)) == IndexCartesian()
    @test @inferred(IndexStyle(typeof(M))) == IndexCartesian()
    @test @inferred(axes(M)) === axes(a)
    @test M[1,1] === RGB{N0f8}(0.1, 0.8, 0)
    @test_throws ErrorException("indexed assignment fails for a reshaped range; consider calling collect") M[1,2] = RGB(0.25, 0.35, 0)

    d = Base.unaliascopy(M)
    @test d == M
    @test d !== M

    a = reshape(0.1:0.1:0.6, 3, 2)
    @test_throws DimensionMismatch mappedarray(f, finv, a, b, c)
end

@testset "Display" begin
    a = [1,2,3,4]
    b = mappedarray(sqrt, a)
    @test summary(b) == "4-element mappedarray(sqrt, ::Array{Int64,1}) with eltype Float64"
    c = mappedarray(sqrt, x->x*x, a)
    @test summary(c) == "4-element mappedarray(sqrt, x->x * x, ::Array{Int64,1}) with eltype Float64"
    # issue #26
    M = @inferred mappedarray((x1,x2)->x1+x2, a, a)
    io = IOBuffer()
    show(io, MIME("text/plain"), M)
    str = String(take!(io))
    @test occursin("x1 + x2", str)
end


text:
# Gurobi model attributes

############################################
#
#   Low level attribute getter/setters
#
############################################

function get_intattr(model::Model, name::String)
    @assert isascii(name)
    a = Array(Cint, 1)
    ret = @grb_ccall(getintattr, Cint,
        (Ptr{Void}, Ptr{UInt8}, Ptr{Cint}),
        model, name, a);
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    convert(Int, a[1])
end

function get_dblattr(model::Model, name::String)
    @assert isascii(name)
    a = Array(Float64, 1)
    ret = @grb_ccall(getdblattr, Cint,
        (Ptr{Void}, Ptr{UInt8}, Ptr{Float64}),
        model, name, a);
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    a[1]::Float64
end

function get_strattr(model::Model, name::String)
    @assert isascii(name)
    a = Array(Ptr{UInt8}, 1)
    ret = @grb_ccall(getstrattr, Cint,
        (Ptr{Void}, Ptr{UInt8}, Ptr{Ptr{UInt8}}),
        model, name, a)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    unsafe_string(a[1])
end

# array element

function get_intattrelement(model::Gurobi.Model, name::String, element::Int)
    @assert isascii(name)
    a = Array(Cint, 1)
    ret = @grb_ccall(getintattrelement, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Cint}),
        model, name, element - 1, a
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    convert(Int, a[1])
end

function get_dblattrelement(model::Gurobi.Model, name::String, element::Int)
    @assert isascii(name)
    a = Array(Float64, 1)
    ret = @grb_ccall(getdblattrelement, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Float64}),
        model, name, element - 1, a
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    a[1]::Float64
end

function get_charattrelement(model::Gurobi.Model, name::String, element::Int)
    @assert isascii(name)
    a = Array(Cchar, 1)
    ret = @grb_ccall(getcharattrelement, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Cchar}),
        model, name, element - 1, a
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    convert(Char, a[1])
end

# Note: in attrarray API, the start argument is one-based (following Julia convention)

function get_intattrarray!(r::Array{Cint}, model::Model, name::String, start::Integer)
    @assert isascii(name)
    ret = @grb_ccall(getintattrarray, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cint, Ptr{Cint}),
        model, name, start - 1, length(r), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    r
end

function get_intattrarray(model::Model, name::String, start::Integer, len::Integer)
    @assert isascii(name)
    get_intattrarray!(Array(Cint, len), model, name, start)
end

function get_dblattrarray!(r::Array{Float64}, model::Model, name::String, start::Integer)
    @assert isascii(name)
    ret = @grb_ccall(getdblattrarray, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cint, Ptr{Float64}),
        model, name, start - 1, length(r), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    r
end

function get_dblattrarray(model::Model, name::String, start::Integer, len::Integer)
    @assert isascii(name)
    get_dblattrarray!(Array(Float64, len), model, name, start)
end

function get_charattrarray!(r::Array{Cchar}, model::Model, name::String, start::Integer)
    @assert isascii(name)
    ret = @grb_ccall(getcharattrarray, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cint, Ptr{Cchar}),
        model, name, start - 1, length(r), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    map(Char,r)
end

function get_charattrarray(model::Model, name::String, start::Integer, len::Integer)
    @assert isascii(name)
    get_charattrarray!(Array(Cchar, len), model, name, start)
end

# setters

function set_intattr!(model::Model, name::String, v::Integer)
    @assert isascii(name)
    ret = @grb_ccall(setintattr, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint), model, name, v)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_dblattr!(model::Model, name::String, v::Real)
    @assert isascii(name)
    ret = @grb_ccall(setdblattr, Cint,
        (Ptr{Void}, Ptr{UInt8}, Float64), model, name, v)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_strattr!(model::Model, name::String, v::String)
    @assert isascii(name)
    @assert isascii(v)
    ret = @grb_ccall(setstrattr, Cint,
        (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}), model, name, v)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

# array element

function set_intattrelement!(model::Gurobi.Model, name::String, element::Int, v::Integer)
    @assert isascii(name)
    ret = @grb_ccall(setintattrelement, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cint),
        model, name, element - 1, v
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_dblattrelement!(model::Gurobi.Model, name::String, element::Int, v::Real)
    @assert isascii(name)
    ret = @grb_ccall(setdblattrelement, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Float64),
        model, name, element - 1, v
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_charattrelement!(model::Gurobi.Model, name::String, element::Int, v::Char)
    @assert isascii(name)
    ret = @grb_ccall(setcharattrelement, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cchar),
        model, name, element - 1, v
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

# reminder: start is one-based

function set_intattrarray!(model::Model, name::String, start::Integer, len::Integer, values::Vector)
    @assert isascii(name)
    values = convert(Vector{Cint},values)
    ret = @grb_ccall(setintattrarray, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cint, Ptr{Cint}), model, name, start-1, len, ivec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_dblattrarray!(model::Model, name::String, start::Integer, len::Integer, values::Vector)
    @assert isascii(name)
    values = convert(Vector{Float64},values)
    ret = @grb_ccall(setdblattrarray, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cint, Ptr{Float64}), model, name, start-1, len, fvec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_charattrarray!(model::Model, name::String, start::Integer, len::Integer, values::Vector)
    @assert isascii(name)
    values = convert(Vector{Cchar},values)
    ret = @grb_ccall(setcharattrarray, Cint,
        (Ptr{Void}, Ptr{UInt8}, Cint, Cint, Ptr{Cchar}), model, name, start-1, len, cvec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

############################################
#
#   Macros for array definition
#
############################################

macro grb_int_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_intattr(model, $(attrname))
end

macro grb_dbl_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_dblattr(model, $(attrname))
end

macro grb_str_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_strattr(model, $(attrname))
end

############################################
#
#   Model attributes
#
############################################

# basic attributes

@grb_str_attr model_name  "ModelName"

@grb_int_attr num_vars     "NumVars"
@grb_int_attr num_constrs  "NumConstrs"
@grb_int_attr num_sos      "NumSOS"
@grb_int_attr num_qconstrs "NumQConstrs"
@grb_int_attr num_cnzs     "NumNZs"
@grb_int_attr num_qnzs     "NumQNZs"
@grb_int_attr num_qcnzs    "NumQCNZs"

@grb_int_attr num_intvars  "NumIntVars"
@grb_int_attr num_binvars  "NumBinVars"

# derived attribute functions

model_sense(model::Model) = get_intattr(model, "ModelSense") > 0 ? (:minimize) : (:maximize)

is_mip(model::Model) = get_intattr(model, "IsMIP") != 0
is_qp(model::Model)  = get_intattr(model, "IsQP") != 0
is_qcp(model::Model) = get_intattr(model, "IsQCP") != 0

function model_type(model::Model)
    is_qp(model)  ? (:QP)  :
    is_qcp(model) ? (:QCP) : (:LP)
end

function set_sense!(model::Model, sense::Symbol)
    v = sense == :maximize ? -1 :
        sense == :minimize ? 1 :
        throw(ArgumentError("Invalid model sense."))

    set_intattr!(model, "ModelSense", v)
end

# variable related attributes

lowerbounds(model::Model) = get_dblattrarray(model, "LB", 1, num_vars(model))
upperbounds(model::Model) = get_dblattrarray(model, "UB", 1, num_vars(model))
objcoeffs(model::Model) = get_dblattrarray(model, "Obj", 1, num_vars(model))

# note: this takes effect only after update_model! is called:
function set_objcoeffs!(model::Model, c::Vector)
    n = num_vars(model)
    length(c) == n || error("Inconsistent argument dimensions.")
    set_dblattrarray!(model, "Obj", 1, n, c)
end

############################################
#
#   The show method for model
#
#   - Based on attributes
#
############################################

function show(io::IO, model::Model)
    if model.ptr_model != C_NULL
        println(io, "Gurobi Model: $(model_name(model))")
        if is_mip(model)
            println(io, "    type   : $(model_type(model)) (MIP)")
        else
            println(io, "    type   : $(model_type(model))")
        end
        println(io, "    sense  : $(model_sense(model))")
        println(io, "    number of variables             = $(num_vars(model))")
        println(io, "    number of linear constraints    = $(num_constrs(model))")
        println(io, "    number of quadratic constraints = $(num_qconstrs(model))")
        println(io, "    number of sos constraints       = $(num_sos(model))")
        println(io, "    number of non-zero coeffs       = $(num_cnzs(model))")
        println(io, "    number of non-zero qp objective terms  = $(num_qnzs(model))")
        println(io, "    number of non-zero qp constraint terms = $(num_qcnzs(model))")
    else
        println(io, "Gurobi Model: NULL")
    end
end