text:
################################################################################
#
#  Morphism types
#
################################################################################

morphism_type(::Type{T}, ::Type{S}) where {R, T <: AbsAlgAss{R}, S <: AbsAlgAss{R}} = AbsAlgAssMor{T, S, Generic.Mat{R}}

morphism_type(::Type{T}, ::Type{S}) where {R, T <: AbsAlgAss{fmpq}, S <: AbsAlgAss{fmpq}} = AbsAlgAssMor{T, S, fmpq_mat}

morphism_type(::Type{T}, ::Type{S}) where {R, T <: AbsAlgAss{fq}, S <: AbsAlgAss{fq}} = AbsAlgAssMor{T, S, fq_mat}

morphism_type(::Type{T}, ::Type{S}) where {R, T <: AbsAlgAss{fq_nmod}, S <: AbsAlgAss{fq_nmod}} = AbsAlgAssMor{T, S, fq_nmod_mat}

morphism_type(::Type{T}, ::Type{S}) where {R, T <: AbsAlgAss{nmod}, S <: AbsAlgAss{nmod}} = AbsAlgAssMor{T, S, nmod_mat}

morphism_type(::Type{T}, ::Type{S}) where {R, T <: AbsAlgAss{gfp_elem}, S <: AbsAlgAss{gfp_elem}} = AbsAlgAssMor{T, S, gfp_mat}

morphism_type(A::Type{T}) where {T <: AbsAlgAss} = morphism_type(A, A)

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

function basis(A::AbsAlgAss)
  if isdefined(A, :basis)
    return A.basis::Vector{elem_type(A)}
  end
  B = Vector{elem_type(A)}(undef, dim(A))
  for i in 1:dim(A)
    z = Vector{elem_type(base_ring(A))}(undef, dim(A))
    for j in 1:dim(A)
      z[j] = zero(base_ring(A))
    end
    z[i] = one(base_ring(A))
    B[i] = A(z)
  end
  A.basis = B
  return B
end

################################################################################
#
#  Predicates
#
################################################################################

issimple_known(A::AbsAlgAss) = A.issimple != 0

################################################################################
#
#  Associativity, Distributivity test
#
################################################################################

function check_associativity(A::AbsAlgAss)
  for i = 1:dim(A)
    for j = 1:dim(A)
      el = A[i] * A[j]
      for k = 1:dim(A)
        if el * A[k] != A[i] * (A[j] * A[k])
          return false
        end
      end
    end
  end
  return true
end

function check_distributivity(A::AbsAlgAss)
  for i = 1:dim(A)
    for j = 1:dim(A)
      el = A[i]*A[j]
      for k = 1:dim(A)
        if A[i] * (A[j] + A[k]) != el + A[i] * A[k]
          return false
        end
      end
    end
  end
  return true
end

################################################################################
#
#  Dimension of center
#
################################################################################

function dimension_of_center(A::AbsAlgAss)
  C, _ = center(A)
  return dim(C)
end

################################################################################
#
#  Subalgebras
#
################################################################################

# Constructs the algebra e*A
# This is the generic fallback which constructs an associative algebra
function subalgebra(A::AbsAlgAss{T}, e::AbsAlgAssElem{T}, idempotent::Bool = false) where {T}
  @assert parent(e) == A
  B, mB = AlgAss(A)
  C, mC = subalgebra(B, mB\e, idempotent)
  mD = compose_and_squash(mB, mC)
  @assert domain(mD) == C
  return C, mD
end

function subalgebra(A::AbsAlgAss{T}, basis::Array{S}) where {T, S}
  B, mB = AlgAss(A)
  basis_pre = elem_type(B)[mB\(basis[i]) for i in 1:length(basis)]
  C, mC = subalgebra(B, basis_pre)
  mD = compose_and_squash(mB, mC)
  @assert domain(mD) == C
  return C, mD
end

################################################################################
#
#  Decomposition
#
################################################################################

# Assume that A is a commutative algebra over a finite field of cardinality q.
# This functions computes a basis for ker(x -> x^q).
function kernel_of_frobenius(A::AbsAlgAss)
  F = base_ring(A)
  q = order(F)

  b = A()
  B = zero_matrix(F, dim(A), dim(A))
  for i = 1:dim(A)
    b.coeffs[i] = one(F)
    if i > 1
      b.coeffs[i - 1] = zero(F)
    end
    c = b^q - b
    for j = 1:dim(A)
      B[j, i] = c.coeffs[j]
    end
  end

  V = right_kernel_basis(B)
  return [ A(v) for v in V ]
end

@doc Markdown.doc"""
    decompose(A::AbsAlgAss{T}) -> AlgAss{T}

Given a semisimple algebra over a field, this function
returns a decomposition of A as a direct sum of simple algebras.
"""
function decompose(A::AbsAlgAss{T}) where {T}
  if isdefined(A, :decomposition)
    return A.decomposition::Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}
  end

  if issimple_known(A) && A.issimple == 1
    B, mB = AlgAss(A)
    res = Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}[(B, mB)]
  end

  if A isa AlgAss
    D = _decompose(A)
    return D
  end

  if A isa AlgGrp || A isa AlgMat
    B, mB = AlgAss(A)
    D = _decompose(B)
    res = Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}[]
    for (S, mS) in D
      mD = compose_and_squash(mB, mS)
      push!(res, (S, mD))
    end
    A.decomposition = res
    return res
  end
end

function _decompose(A::AbsAlgAss{T}) where {T}
  if iscommutative(A)
    res = _dec_com(A)
  else
    res = _dec_via_center(A)
  end

  A.decomposition = res
  return res
end

function _dec_via_center(A::S) where {T, S <: AbsAlgAss{T}}
  ZA, mZA = center(A)
  Algs = _dec_com(ZA)
  ZA.decomposition = Algs
  res = Tuple{AlgAss{T}, morphism_type(AlgAss{T}, S)}[ subalgebra(A, mZA(BtoZA(one(B))), true) for (B, BtoZA) in Algs]
  for i in 1:length(res)
    res[i][1].issimple = 1
    B, BtoZA = Algs[i] # B is the centre of res[i][1]
    # Build a map from B to res[i][1] via B -> ZA -> A -> res[i][1]
    M = zero_matrix(base_ring(A), dim(B), dim(res[i][1]))
    for j = 1:dim(B)
      t = mZA(BtoZA(B[j]))
      s = res[i][2]\t
      elem_to_mat_row!(M, j, s)
    end
    if dim(res[i][1]) != dim(B)
      res[i][1].center = (B, hom(B, res[i][1], M))
    else
      # res[i][1] is commutative, so we do not cache the centre
      iM = inv(M)
      BtoA = hom(B, A, M*res[i][2].mat, res[i][2].imat*iM)
      res[i] = (B, BtoA)
    end
  end
  A.decomposition = res
  return res
end

function _dec_com(A::AbsAlgAss)
  if characteristic(base_ring(A)) > 0
    return _dec_com_finite(A)
  else
    return _dec_com_gen(A)
  end
end

function _dec_com_gen(A::AbsAlgAss{T}) where {T <: FieldElem}
  if dim(A) == 1
    A.issimple = 1
    B, mB = AlgAss(A)
    return Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}[(B, mB)]
  end

  F = base_ring(A)

  k = dim(A)

  V = elem_type(A)[A[i] for i in 1:k]

  while true
    c = elem_type(F)[ rand(F, -10:10) for i = 1:k ]
    a = dot(c, V)
    f = minpoly(a)

    if degree(f) < 2
      continue
    end
    if isirreducible(f)
      if degree(f) == dim(A)
        A.issimple = 1
        B, mB = AlgAss(A)
        return Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}[(B, mB)]
      end
      continue
    end

    @assert issquarefree(f)

    fac = factor(f)
    R = parent(f)
    factors = Vector{elem_type(R)}()
    for ff in keys(fac.fac)
      push!(factors, ff)
    end
    sols = Vector{elem_type(R)}()
    right_side = elem_type(R)[ R() for i = 1:length(factors) ]
    max_deg = 0
    for i = 1:length(factors)
      right_side[i] = R(1)
      if i != 1
        right_side[i - 1] = R(0)
      end
      s = crt(right_side, factors)
      push!(sols, s)
      max_deg = max(max_deg, degree(s))
    end
    x = one(A)
    powers = Vector{elem_type(A)}()
    for i = 1:max_deg + 1
      push!(powers, x)
      x *= a
    end
    idems = Vector{elem_type(A)}()
    for s in sols
      idem = A()
      for i = 0:degree(s)
        idem += coeff(s, i)*powers[i + 1]
      end
      push!(idems, idem)
    end

    A.issimple = 2

    res = Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}()
    for idem in idems
      S, StoA = subalgebra(A, idem, true)
      decS = _dec_com_gen(S)
      for (B, BtoS) in decS
        BtoA = compose_and_squash(StoA, BtoS)
        push!(res, (B, BtoA))
      end
    end
    return res
  end
end

function _dec_com_finite(A::AbsAlgAss{T}) where T
  if dim(A) == 1
    A.issimple = 1
    B, mB = AlgAss(A)
    return Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}[(B, mB)]
  end

  F = base_ring(A)
  @assert !iszero(characteristic(F))
  V = kernel_of_frobenius(A)
  k = length(V)

  if k == 1
    A.issimple = 1
    B, mB = AlgAss(A)
    return Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}[(B, mB)]
  end

  A.issimple = 2
  c = elem_type(F)[ rand(F) for i = 1:k ]
  M = zero_matrix(F, dim(A), dim(A))
  a = dot(c, V)
  representation_matrix!(a, M)
  f = minpoly(M)
  while degree(f) < 2
    for i = 1:length(c)
      c[i] = rand(F)
    end
    a = dot(c, V)
    zero!(M)
    representation_matrix!(a, M)
    f = minpoly(M)
  end

  #@assert issquarefree(f)
  fac = factor(f)
  R = parent(f)
  factorss = collect(keys(fac.fac))
  sols = Vector{typeof(f)}(undef, length(factorss))
  right_side = typeof(f)[ zero(R) for i = 1:length(factorss) ]
  max_deg = 0
  for i = 1:length(factorss)
    right_side[i] = one(R)
    if 1 != i
      right_side[i - 1] = zero(R)
    end
    sols[i] = crt(right_side, factorss)
    max_deg = max(max_deg, degree(sols[i]))
  end
  powers = Vector{elem_type(A)}(undef, max_deg+1)
  powers[1] = one(A)
  powers[2] = a
  x = a
  for i = 3:max_deg + 1
    x *= a
    powers[i] = x
  end

  idems = Vector{elem_type(A)}()
  for s in sols
    idem = A()
    for i = 0:degree(s)
      idem += coeff(s, i)*powers[i + 1]
    end
    push!(idems, idem)
  end

  res = Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}()
  for idem in idems
    S, StoA = subalgebra(A, idem, true)
    decS = _dec_com_finite(S)
    for (B, BtoS) in decS
      BtoA = compose_and_squash(StoA, BtoS)
      push!(res, (B, BtoA))
    end
  end
  return res

end

################################################################################
#
#  Decomposition as number fields
#
################################################################################

@doc Markdown.doc"""
    as_number_fields(A::AbsAlgAss{fmpq})

Given a commutative algebra over QQ, this function returns a decomposition
of A as direct sum of number fields.
"""
function as_number_fields(A::AbsAlgAss{fmpq})
  if isdefined(A, :maps_to_numberfields)
    return A.maps_to_numberfields::Vector{Tuple{AnticNumberField, AbsAlgAssToNfAbsMor{typeof(A), elem_type(A)}}}
  end

  d = dim(A)

  Adec = decompose(A)

  fields_not_cached = false
  for i = 1:length(Adec)
    if !isdefined(Adec[i][1], :maps_to_numberfields)
      fields_not_cached = true
    end
  end

  if fields_not_cached
    # Compute a LLL reduced basis of the maximal order of A to find "small"
    # polynomials for the number fields.
    OA = maximal_order(A)
    L = lll(basis_mat(OA, copy = false).num)
    n = basis_mat(OA, copy = false).den
    basis_lll = [ elem_from_mat_row(A, L, i, n) for i = 1:d ]
  end

  M = zero_matrix(FlintQQ, 0, d)
  matrices = Vector{fmpq_mat}()
  fields = Vector{AnticNumberField}()
  for i = 1:length(Adec)
    # For each small algebra construct a number field and the isomorphism
    B, BtoA = Adec[i]
    dB = dim(B)
    if !isdefined(B, :maps_to_numberfields)
      local K, BtoK
      found_field = false # Only for debugging
      for j = 1:d
        t = BtoA\basis_lll[j]
        mint = minpoly(t)
        if degree(mint) == dB
          found_field = true
          K, BtoK = _as_field_with_isomorphism(B, t, mint)
          B.maps_to_numberfields = Tuple{AnticNumberField, AbsAlgAssToNfAbsMor{typeof(B), elem_type(B)}}[(K, BtoK)]
          push!(fields, K)
          break
        end
      end
      @assert found_field "This should not happen..."
    else
      K, BtoK = B.maps_to_numberfields[1]
      push!(fields, K)
    end

    if length(Adec) == 1
      A.maps_to_numberfields = Tuple{AnticNumberField, AbsAlgAssToNfAbsMor{typeof(A), elem_type(A)}}[(K, BtoK)]
      return A.maps_to_numberfields
    end

    # Construct the map from K to A
    N = zero_matrix(FlintQQ, degree(K), d)
    for j = 1:degree(K)
      t = BtoA(BtoK\basis(K)[j])
      elem_to_mat_row!(N, j, t)
    end
    push!(matrices, N)
    M = vcat(M, N)
  end
  @assert nrows(M) == d

  invM = inv(M)
  matrices2 = Vector{fmpq_mat}(undef, length(matrices))
  offset = 1
  for i = 1:length(matrices)
    r = nrows(matrices[i])
    N = sub(invM, 1:d, offset:(offset + r - 1))
    matrices2[i] = N
    offset += r
  end

  result = Vector{Tuple{AnticNumberField, AbsAlgAssToNfAbsMor{typeof(A), elem_type(A)}}}()
  for i = 1:length(fields)
    push!(result, (fields[i], AbsAlgAssToNfAbsMor(A, fields[i], matrices2[i], matrices[i])))
  end
  A.maps_to_numberfields = result
  return result
end

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

function rand(A::AbsAlgAss{T}) where T
  c = T[rand(base_ring(A)) for i = 1:dim(A)]
  return A(c)
end

function rand(A::AbsAlgAss{T}, rng::UnitRange{Int}) where T
  c = T[rand(base_ring(A), rng) for i = 1:dim(A)]
  return A(c)
end

function rand(A::AlgAss{fmpq}, rng::UnitRange{Int} = -20:20)
  c = [fmpq(rand(FlintZZ, rng)) for i = 1:dim(A)]
  return A(c)
end

################################################################################
#
#  Generators
#
################################################################################

function gens(A::AbsAlgAss, return_full_basis::Type{Val{T}} = Val{false}) where T
  K = base_ring(A)
  d = dim(A)

  rfb = return_full_basis == Val{true}

  # Sort the basis by the degree of the minpolys (hopefully those with higher
  # degree generate a "bigger" subalgebra)
  minpoly_degrees = [ (i, degree(minpoly(A[i]))) for i = 1:d ]
  sort!(minpoly_degrees, by = x -> x[2], rev = true)

  generators = Vector{elem_type(A)}()
  full_basis = elem_type(A)[ one(A) ] # Contains products of generators which form a full basis
  rfb ? full_basis_indices = Vector{Tuple{Int, Int}}[ Tuple{Int, Int}[] ] : nothing
  B = zero_matrix(K, d, d)

  cur_dim = 0
  for i = 1:d
    if cur_dim == d
      break
    end

    n = length(full_basis)
    b = A[minpoly_degrees[i][1]]
    power = 1
    while cur_dim < d
      for k = 1:d
        B[d, k] = coeffs(b, copy = false)[k]
      end
      new_dim = rref!(B)
      if cur_dim == new_dim
        break
      end
      if power == 1
        push!(generators, b)
      end
      push!(full_basis, b)
      if rfb
        ind = Tuple{Int, Int}[ (length(generators), power) ]
        push!(full_basis_indices, ind)
      end
      cur_dim = new_dim
      cur_dim == d ? break : nothing

      for r = 1:n
        bb = b*full_basis[r]
        for l = 1:n
          t = full_basis[l]*bb
          for k = 1:d
            B[d, k] = coeffs(t, copy = false)[k]
          end
          new_dim = rref!(B)
          if cur_dim == new_dim
            continue
          end
          push!(full_basis, t)
          cur_dim = new_dim
          if rfb
            ind2 = deepcopy(ind)
            prepend!(ind2, full_basis_indices[l])
            append!(ind2, full_basis_indices[r])
            push!(full_basis_indices, ind2)
          end
          cur_dim == d ? break : nothing
        end
        cur_dim == d ? break : nothing
      end
      b *= A[minpoly_degrees[i][1]]
      power += 1
    end
  end

  # Remove the one
  popfirst!(full_basis)
  rfb ? popfirst!(full_basis_indices) : nothing

  if rfb
    return generators, full_basis, full_basis_indices
  else
    return generators
  end
end

################################################################################
#
#  Primitive elements
#
################################################################################

function primitive_element(A::AbsAlgAss)
  a, _ = _primitive_element(A)
  return a
end

function _primitive_element(A::AbsAlgAss)
  error("Not implemented yet")
  return nothing
end

function _primitive_element(A::AbsAlgAss{T}) where T <: Union{nmod, fq, fq_nmod, Generic.Res{fmpz}, fmpq, Generic.ResF{fmpz}, gfp_elem}
  d = dim(A)
  a = rand(A)
  f = minpoly(a)
  while degree(f) < d
    a = rand(A)
    f = minpoly(a)
  end
  return a, f
end

function _as_field(A::AbsAlgAss{T}) where T
  d = dim(A)
  a, mina = _primitive_element(A)
  b = one(A)
  M = zero_matrix(base_ring(A), d, d)
  elem_to_mat_row!(M, 1, b)
  for i in 1:(d - 1)
    b = mul!(b, b, a)
    elem_to_mat_row!(M, i + 1, b)
  end
  B = inv(M)
  N = zero_matrix(base_ring(A), 1, d)
  local f
  let N = N, B = B
    f = function(x)
      for i in 1:d
        N[1, i] = x.coeffs[i]
      end
      return N * B
    end
  end
  return a, mina, f
end

function _as_field_with_isomorphism(A::AbsAlgAss{S}) where { S <: Union{fmpq, gfp_elem, Generic.ResF{fmpz}, fq_nmod, fq} }
  return _as_field_with_isomorphism(A, _primitive_element(A)...)
end

# Assuming a is a primitive element of A and mina its minimal polynomial, this
# functions constructs the field base_ring(A)/mina and the isomorphism between
# A and this field.
function _as_field_with_isomorphism(A::AbsAlgAss{S}, a::AbsAlgAssElem{S}, mina::T) where { S <: Union{fmpq, gfp_elem, Generic.ResF{fmpz}, fq_nmod, fq}, T <: Union{fmpq_poly, gfp_poly, gfp_fmpz_poly, fq_nmod_poly, fq_poly} }
  s = one(A)
  M = zero_matrix(base_ring(A), dim(A), dim(A))
  elem_to_mat_row!(M, 1, s)
  for i = 2:dim(A)
    s = mul!(s, s, a)
    elem_to_mat_row!(M, i, s)
  end

  if base_ring(A) == FlintQQ
    K = number_field(mina, cached = false)[1]
    return K, AbsAlgAssToNfAbsMor(A, K, inv(M), M)
  elseif base_ring(A) isa GaloisField
    Fq = FqNmodFiniteField(mina, Symbol("a"), false)
    return Fq, AbsAlgAssToFqMor(A, Fq, inv(M), M, parent(mina))
  elseif base_ring(A) isa Generic.ResField{fmpz}
    Fq = FqFiniteField(mina, Symbol("a"), false)
    return Fq, AbsAlgAssToFqMor(A, Fq, inv(M), M, parent(mina))
  elseif base_ring(A) isa FqNmodFiniteField || base_ring(A) isa FqFiniteField
    Fr, RtoFr = field_extension(mina)
    return Fr, AbsAlgAssToFqMor(A, Fr, inv(M), M, parent(mina), RtoFr)
  else
    error("Not implemented")
  end
end

################################################################################
#
#  Regular matrix algebra
#
################################################################################

@doc Markdown.doc"""
    regular_matrix_algebra(A::Union{ AlgAss, AlgGrp })

Returns the matrix algebra B generated by the right representation matrices of
the basis elements of A and a map from B to A.
"""
function regular_matrix_algebra(A::Union{ AlgAss, AlgGrp })
  K = base_ring(A)
  B = AlgMat(K, [ representation_matrix(A[i], :right) for i = 1:dim(A) ], isbasis = true, check = false)
  return B, hom(B, A, identity_matrix(K, dim(A)), identity_matrix(K, dim(A)))
end

###############################################################################
#
#  Construction of a crossed product algebra
#
###############################################################################

function find_elem(G::Array{T,1}, el::T) where T
  i=1
  while true
    if el.prim_img==G[i].prim_img
      return i
    end
    i+=1
  end
end

#K/Q is a Galois extension.
function CrossedProductAlgebra(K::AnticNumberField, G::Array{T,1}, cocval::Array{nf_elem, 2}) where T

  n=degree(K)
  m=length(G)
  #=
  Multiplication table
  I order the basis in this way:
  First, I put the basis of the Galois Group, then the product of them with the first
  element of basis of the order and so on...
  =#

  M=Array{fmpq,3}(undef, n*m, n*m, n*m)
  for i=1:n*m
    for j=1:n*m
      for s=1:n*m
        M[i,j,s]=fmpq(0)
      end
    end
  end
  B=basis(K)
  for i=1:n
    for j=1:m
      #I have the element B[i]*G[j]
      for k=1:n
        for h=1:m
          # I take the element B[k]*G[h]
          # and I take the product
          # B[i]*G[j]* B[k]*G[h]=B[i]*G[j](B[k])*c[j,h]*(G[j]*G[h])
          ind=find_elem(G,G[h] * G[j])
          x=B[i]*G[j](B[k])*cocval[j,h]
          #@show i, j, k,h,  ind,B[i],G[j](B[k]),cocval[j,h],  x
          for s=0:n-1
            M[j+(i-1)*n, h+(k-1)*n, ind+s*n]=coeff(x,s)
          end
          #@show M
        end
      end
    end
  end
  return AlgAss(FlintQQ, M)

end

function CrossedProductAlgebra(O::NfOrd, G::Array{T,1}, cocval::Array{nf_elem, 2}) where T

  n=degree(O)
  m=length(G)
  K=nf(O)
  #=
  Multiplication table
  I order the basis in this way:
  First, I put the basis of the Galois Group, then the product of them with the first
  element of basis of the order and so on...
  =#

  M=Array{fmpq,3}(undef, n*m, n*m, n*m)
  for i=1:n*m
    for j=1:n*m
      for s=1:n*m
        M[i,j,s]=fmpq(0)
      end
    end
  end
  B = basis(O, copy = false)
  el = O(0)
  for j=1:m
    for k=1:n
      l =O(G[j](K(B[k])), false)
      for h=1:m
        ind = find_elem(G, G[h] * G[j])
        t = O(cocval[j,h], false)
        for i=1:n
          #I have the element B[i]*G[j]
          # I take the element B[k]*G[h]
          # and I take the product
          # B[i]*G[j]* B[k]*G[h]=B[i]*G[j](B[k])*c[j,h]*(G[j]*G[h])
          mul!(el, B[i], l)
          mul!(el, el, t)
          y = coordinates(el)
          for s=0:n-1
            M[j+(i-1)*m, h+(k-1)*m, ind+s*m] = y[s+1]
          end
        end
      end
    end
  end
  j1 = find_identity(G, *)
  j = find_elem(G, j1)
  O1 = fmpq[0 for i=1:n*m]
  O1[j] = fmpq(1)
  A = AlgAss(FlintQQ, M, O1)
  A.issimple = 1
  return A

end

################################################################################
#
#  Quaternion algebras
#
################################################################################

function quaternion_algebra(a::Int, b::Int)

  M = Array{fmpq,3}(undef, 4,4,4)
  for i = 1:4
    for j = 1:4
      for k = 1:4
        M[i,j,k] = 0
      end
    end
  end
  M[1,1,1] = 1 # 1*1=1
  M[1,2,2] = 1 # 1*i=i
  M[1,3,3] = 1 # 1*j=j
  M[1,4,4] = 1 # 1*ij=1

  M[2,1,2] = 1
  M[2,2,1] = a
  M[2,3,4] = 1
  M[2,4,3] = a

  M[3,1,3] = 1
  M[3,2,4] = -1
  M[3,3,1] = b
  M[3,4,2] = -b

  M[4,1,4] = 1
  M[4,2,3] = -a
  M[4,3,2] = b
  M[4,4,1] = -a*b
  O = fmpq[1, 0, 0, 0]
  return AlgAss(FlintQQ, M, O)

end


text:
using Test

using ABIF
using FormatSpecimens

@testset "ABIF Reader" begin

    function check_abif_parse(filename)
        stream = open(ABIF.Reader, filename)

        for record in stream end
        for (a,b) in collect(stream[1]) end
        for (a,b) in getindex(stream, ABIF.get_tags(stream)) end
        @test typeof(stream) == ABIF.Reader{IOStream}
        @test ABIF.get_tags(stream)[1].name == "AEPt"
        @test ABIF.get_tags(stream, "DATA")[1].name == "DATA"
        @test length(stream["DATA"]) == 12
        @test length(stream[1]) == 1

        @test typeof(getindex(stream, ABIF.get_tags(stream))) == Dict{String,Any}
        @test ABIF.tagelements(stream, "DATA") == 12
    end

    println("Good ab1:")
    for file in list_valid_specimens("ABI")
        filepath = joinpath(path_of_format("ABI"), filename(file))
        println("\t",filepath)
        check_abif_parse(filepath)
    end
    println("Bad ab1:")
    for file in list_invalid_specimens("ABI")
        filepath = joinpath(path_of_format("ABI"), filename(file))
        println("\t",filepath)
        @test_throws Exception check_abif_parse(filepath)
    end

end

ext: jl


text:
"""
    loadfunctionfromfile(filename, functionname[, injection; modulename])

Loads the specified function from the given file.
Encapsulates the loaded file in a module named `modulename` (which exists in
the `Main` module namespace.)

**Note:** This function is potentially very dangerous; the file is loaded in
entirety, including all side effects.
This means that the loaded file may try to do all kinds of naughty things
(like call `exit()`, overwrite files, etc.).
It should **never** be used outside of a sandbox.
"""
function loadfunctionfromfile(
    filename,
    functionname,
    injection::Union{Expr, Nothing}=nothing;
    modulename=Symbol("_Submission"),
)
    functionname = Symbol(functionname)
    modulename = Symbol(modulename)
    # Note: this function should only be run once during a given Julia process.
    # We can't have specific _Submission namespaces since that would result in
    # things not lining up between common code in the submission and solution.
    # e.g. if we have the Fab type, we need to have _Submission.Fab rather than
    # _Submission4561.Fab so that it has the same name everywhere.
    moduleexpr = :(
        module $modulename
        $(injection)

        include($filename)
        end
    )
    try
        Main.eval(moduleexpr)
        # Note that we need to escape the symbol with `QuoteNode` to prevent
        # the symbol from being directly (literally) interpolated.
        # We want `getfield(_Submission, :functionname)` instead of
        # `getfield(_Submission, functionname)` which will raise an error about
        # functionname not being defined.
        func = @eval Main getfield($modulename, $(QuoteNode(functionname)))
        return func
    catch e
        @debug(
            "Caught exception when trying to load function from file.",
            e, filename, functionname,
        )
        error("Couldn't find function $functionname in file $filename.")
    end
end


text:
module TensorCast

using MacroTools

const V = false # verbose debugging

include("parse.jl")

include("macro.jl") # @cast and @reduce

include("matmul.jl") # @mul

include("icheck.jl")

include("string.jl")

include("pretty.jl")

if VERSION < v"1.1.0"
    include("eachslice.jl")
end

include("slice.jl") # slice, glue, orient, etc

include("recursive.jl") # RecursiveArrayTools

include("order.jl") # Reverse{D} etc

include("lazy.jl") # LazyCast

using Requires

function __init__()

    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" begin
        include("static.jl")
    end

    # @require Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" begin
    #     include("strided.jl")
    # end

    @require JuliennedArrays = "5cadff95-7770-533d-a838-a1bf817ee6e0" begin
        include("julienne.jl")
    end

    # @require LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" begin
    # end

    # @require Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" begin
    #   include("flux.jl")
    # end

    # @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" begin
    #     include("zygote.jl")
    # end

    # @require NamedDims = "356022a1-0364-5f58-8944-0da4b18d706f" begin
    #     include("nameddims.jl")
    # end

    @require NamedArrays = "86f7a689-2022-50b4-a561-43c23ac3c673" begin
        include("namedarrays.jl")
    end

    # @require AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" begin
    #     include("axis.jl")
    # end

end

end # module


text:
# Add fallbacks for predict_* which require mean, mode, median.

const BadMeanTypes   = Union{AbstractArray{<:Finite},Table(Finite)}
const BadMedianTypes = Union{AbstractArray{<:Finite},Table(Finite)}

const err_wrong_target_scitype(actual_scitype) = ArgumentError(
    "Attempting to compute mode of predictions made "*
    "by a model with `$actual_scitype` targets. ")

# mode:
predict_mode(m, fitresult, Xnew) =
    mode.(predict(m, fitresult, Xnew))

# mean:
predict_mean(m, fitresult, Xnew) =
    predict_mean(m, fitresult, Xnew, target_scitype(m))
predict_mean(m, fitresult, Xnew, ::Any) =
    mean.(predict(m, fitresult, Xnew))
predict_mean(m, fitresult, Xnew, ::Type{<:BadMeanTypes}) =
    throw(err_wrong_target_scitype(Finite))

# median:
predict_median(m, fitresult, Xnew) =
    predict_median(m, fitresult, Xnew, target_scitype(m))
predict_median(m, fitresult, Xnew, ::Any) =
    median.(predict(m, fitresult, Xnew))
predict_median(m, fitresult, Xnew, ::Type{<:BadMedianTypes}) =
    throw(err_wrong_target_scitype(Finite))

# not in MLJModelInterface as methodswith requires InteractiveUtils
MLJModelInterface.implemented_methods(::FI, M::Type{<:MLJType}) =
    getfield.(methodswith(M), :name) |> unique

# The following serialization fallbacks should live in
# MLJModelInterface when version 2.0 is released. At that time the
# hack block could also be removed.

#####################
# hack block begins #
#####################
const ERR_SERIALIZATION_FAILURE = ErrorException(
"Serialization failure. You are using a model that implements an outdated "*
    "version of the serialization API. If you are using "*
    "a model from XGBoost.jl, try using MLJXGBoostInterface 2.0 or "*
    "or higher. "
)
const ERR_DESERIALIZATION_FAILURE = ErrorException(
    "Deserialization failure. Your model must be deserialized using "*
    "using MLJBase < 0.20 and MLJSerialization < 2.0. If this is an "*
    "XGBoost.jl model, be sure to use MLJXGBoostInterface < 2.0. "
)
MLJModelInterface.save(filename, model, fitresult; kwargs...) =
    throw(ERR_SERIALIZATION_FAILURE)
MLJModelInterface.restore(filename, model, serializable_fitresult) =
    throw(ERR_DESERIALIZATION_FAILURE)
###################
# hack block ends #
###################

MLJModelInterface.save(model, fitresult; kwargs...) = fitresult
MLJModelInterface.restore(model, serializable_fitresult) =
    serializable_fitresult

# to suppress inclusion of abstract types in the model registry.
for T in (:Supervised, :Unsupervised,
          :Interval, :Static, :Deterministic, :Probabilistic)
    ex = quote
        MLJModelInterface.is_wrapper(::Type{$T}) = true
    end
    eval(ex)
end


text:
function one_plus_one(layers::Array, num_generations::Int, f::Function)
    fits = zeros(num_generations)
    parent = Individual(layers::Array)
    evaluate!(parent, f)

    for i in eachindex(fits)
        child = mutate(parent)
        evaluate!(child, f)

        if child.fitness >= parent.fitness
            parent = child
        end

        fits[i] = parent.fitness
        print(i, "/",  num_generations, " > ", parent.fitness, "\t\t\r")
    end
    fits, parent
end

max_stars_repo_name: TemplierPaul/evolution

max_issues_repo_name: TemplierPaul/evolution

max_forks_repo_name: TemplierPaul/evolution


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON
# because of vae and 2stage
using DataFrames
using CSV
using ValueHistories
using Flux
using ConditionalDists
using GenerativeModels
import GenerativeModels: VAE
using Distributions
using DistributionsAD

s = ArgParseSettings()
@add_arg_table! s begin
	"max_seed"
		default = 1
		arg_type = Int
		help = "max_seed"
	"category"
		default = "wood"
		arg_type = String
		help = "category"
	"tab_name"
		default = "vae_LOSS_images_mvtec"
		arg_type = String
		help = "name of tab -> example: vae_LOSS_images, wae-vamp_AUC_images"
	"mi_only"
		arg_type = Int
		default = -1
		help = "index of model in range 1 to 10 or -1 for all models"
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack category, max_seed, tab_name, mi_only, contamination = parsed_args

#######################################################################################
################ THIS PART IS TO BE PROVIDED FOR EACH MODEL SEPARATELY ################

sp = split(tab_name, "_")
enc = sp[1]
criterion = lowercase(sp[2])

modelname = "$(enc)_ocsvm"

function sample_params()
	par_vec = (round.([10^x for x in -4:0.1:2], digits=5),["poly", "rbf", "sigmoid"],[0.01,0.5,0.99])
	argnames = (:gamma,:kernel,:nu)
	return (;zip(argnames, map(x->sample(x, 1)[1], par_vec))...)
end

function joint_fit(models, data_splits)
	info = []
	for (model, data) in zip(models, data_splits)
		push!(info, fit!(model, data))
	end
	return info
end

function joint_prediction(models, data)
	joint_pred = Array{Float32}(undef, length(models), size(data,2))
	for (i,model) in enumerate(models)
		joint_pred[i,:] = predict(model, data)
	end
	return vec(mean(joint_pred, dims=1))
end

function StatsBase.fit(data, parameters, n_models, aux_info)
	# construct model - constructor should only accept kwargs
	models = [GenerativeAD.Models.OCSVM(;parameters...) for _ = 1:n_models]

	# sumbsample and fit train data
	tr_data = data[1][1]
	M,N = size(tr_data)
	split_size = floor(N/n_models)
	data_splits = [tr_data[:,Int(split_size*i+1):Int(split_size*(i+1))] for i = 0:n_models-1]

	try
		global info, fit_t, _, _, _ = @timed joint_fit(models, data_splits)
 	catch e
		# return an empty array if fit fails so nothing is computed
		return (fit_t = NaN,), []
	end

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = fit_t,
		model = nothing
		)

	# now return the different scoring functions
	training_info, [(x->joint_prediction(models, x),  merge(parameters, aux_info))]
end

####################################################################
################ THIS PART IS COMMON FOR ALL MODELS ################
# set a maximum for parameter sampling retries
if abspath(PROGRAM_FILE) == @__FILE__
	try_counter = 0
	max_tries = 1
	cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
	while try_counter < max_tries
		parameters = sample_params()

		for seed in 1:max_seed
			savepath = datadir("experiments/images_mvtec$cont_string/$(modelname)/$(category)/ac=1/seed=$(seed)")
			mkpath(savepath)

			mi_indexes = (mi_only == -1) ? [1:10...] : [mi_only]

			for mi = mi_indexes

				aux_info = (model_index=mi, criterion=criterion)
				data = GenerativeAD.load_data("MVTec-AD", seed=seed, category=category,
					contamination=contamination, img_size=128)
				data, encoding_name, encoder_params = GenerativeAD.Models.load_encoding(tab_name, data, 1, dataset=category, seed=seed, model_index=mi)

				@info "Trying to fit $modelname on $category with parameters $(parameters)..."
				@info "Train/validation/test splits: $(size(data[1][1])) | $(size(data[2][1])) | $(size(data[3][1]))"
				@info "Number of features: $(size(data[1][1]))"
				# here, check if a model with the same parameters was already tested
				@info "Trying to fit $modelname on $category with parameters $(parameters)..."
				if GenerativeAD.check_params(savepath, merge(parameters, aux_info))
					training_info, results = fit(data, parameters, 10, aux_info)
					if training_info.model !== nothing
						tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)),
							Dict("model"=>training_info.model),
							safe = true)
						training_info = merge(training_info, (model = nothing,))
					end
					save_entries = merge(training_info, (modelname = modelname, seed = seed,
						category = category,
						contamination=contamination))
					# now loop over all anomaly score funs
					for result in results
						GenerativeAD.experiment(result..., data, savepath; save_entries...)
					end
					global try_counter = max_tries + 1
				else
					@info "Model already present, sampling new hyperparameters..."
					global try_counter += 1
				end
			end
		end
	end
	(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing
end


text:
using NumCME
using BenchmarkTools
import DifferentialEquations as DE
using Sundials: CVODE_BDF
using LinearAlgebra: BLAS

##  Toggle-switch model
𝕊 = [[1, 0] [-1, 0] [0, 1] [0, -1]]

α₁ = propensity() do x, p
    S₁, S₂ = x
    b₁, b₂, k₁, k₂, a₂₁, a₁₂, n₂₁, n₁₂, γ₁, γ₂, UV, Δtᵤᵥ = p
    b₁ + k₁ / (1.0 + a₂₁ * S₂^n₂₁)
end
α₂ = propensity() do x, p
    S₁, S₂ = x
    b₁, b₂, k₁, k₂, a₂₁, a₁₂, n₂₁, n₁₂, γ₁, γ₂, UV, Δtᵤᵥ = p
    γ₁ * S₁
end
α₃ = propensity() do x, p
    S₁, S₂ = x
    b₁, b₂, k₁, k₂, a₂₁, a₁₂, n₂₁, n₁₂, γ₁, γ₂, UV, Δtᵤᵥ = p
    b₂ + k₂ / (1.0 + a₁₂ * S₁^n₁₂)
end
α₄ = propensity() do t, x, p
    S₁, S₂ = x
    b₁, b₂, k₁, k₂, a₂₁, a₁₂, n₂₁, n₁₂, γ₁, γ₂, UV, Δtᵤᵥ = p
    (γ₂ + (t ≤ Δtᵤᵥ) * 0.002 * UV^2 / (1260 + UV^3)) * S₂
end

# This propensity formulation is mathematically equivalent to α₄ but leads to more computationally efficient CME solves because β₄ is factored into a time-only and a state-only functions
function degradation_rate(t, p)
    b₁, b₂, k₁, k₂, a₂₁, a₁₂, n₂₁, n₁₂, γ₁, γ₂, UV, Δtᵤᵥ = p
    γ₂ + (t ≤ Δtᵤᵥ) * 0.002 * UV^2 / (1260 + UV^3)
end
β₄ = propensity(degradation_rate) do x, p
    x[2]
end

b₁ = 2.2E-3
b₂ = 6.8E-5
k₁ = 1.7E-2
k₂ = 1.6E-2
a₂₁ = 2.6E-3
a₁₂ = 6.1E-3
n₂₁ = 3
n₁₂ = 2.1
γ₁ = 3.8E-4
γ₂ = 3.8E-4
UV = 10.0
Δtᵤᵥ = 3600

θ = [b₁, b₂, k₁, k₂, a₂₁, a₁₂, n₂₁, n₁₂, γ₁, γ₂, UV, Δtᵤᵥ]

propensities_joint = [α₁, α₂, α₃, α₄]
propensities_separable = [α₁, α₂, α₃, β₄]
model_joint = CmeModel(𝕊, propensities_joint, θ)
model_separable = CmeModel(𝕊, propensities_separable, θ)

x₀ = [0, 0]
𝔛₀ = StateSpaceSparse(𝕊, x₀)
p0 = FspVectorSparse(𝔛₀, [x₀ => 1.0])

tspan = (0.0, 8.0 * 3600)
saveat = 0.0:60.0:8*3600.0

fixedrstepfsp = AdaptiveFspSparse(
    ode_method = CVODE_BDF(linear_solver = :GMRES),
    space_adapter = RStepAdapter(20, 5, true)
)

adaptiverstepfsp = AdaptiveFspSparse(
    ode_method = CVODE_BDF(linear_solver = :GMRES),
    space_adapter = SelectiveRStepAdapter(20, 5, true)
)

for threadcount in [1, 4, 8]
    BLAS.set_num_threads(threadcount)
    println("Solving with full R-step expansion and separable propensity format")
    @btime fspsol1 = solve(model_separable, p0, tspan, fixedrstepfsp, saveat = saveat, odertol = 1.0E-4, odeatol = 1.0E-14)
    println("Solving with selective R-step expansion and separable propensity format")
    @btime fspsol2 = solve(model_separable, p0, tspan, adaptiverstepfsp, saveat = saveat, odertol = 1.0E-4, odeatol = 1.0E-14)
    println("Solving with full R-step expansion and non-separable propensity format")
    @btime fspsol3 = solve(model_joint, p0, tspan, fixedrstepfsp, saveat = saveat, odertol = 1.0E-4, odeatol = 1.0E-14)
    println("Solving with selective R-step expansion and non-separable propensity format")
    @btime fspsol4 = solve(model_joint, p0, tspan, adaptiverstepfsp, saveat = saveat, odertol = 1.0E-4, odeatol = 1.0E-14)
end


text:
using ConductanceBasedModels
using Test

@testset "ConductanceBasedModels.jl" begin
    # Write your tests here.
end


text:
# Our colleague's project depends on a useless package
# (https://github.com/crstnbr/Useless.jl)
using Useless

life_universe_everything(question) = 42

println(life_universe_everything("What is the meaning of life?"))

# Step 1: instantiate the environment
# Step 2: run the code: julia --project=. code.jl


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# # Example 2
#
# Two-dimensional fields are read from the netcdf file generated in `Example1`, and then re-written to a new netcdf file.

# +
using NCTiles,NCDatasets

inputs = "input/"
outputs = "output/"

file_in=outputs*"ex1/ex1_NCDatasets.nc"
~isfile(file_in) ? error("Running Example1 first is needed to run Example2") : nothing

file_out=outputs*"ex2/ex2_NCDatasets.nc"
if ~ispath(outputs*"ex2/"); mkpath(outputs*"ex2/"); end

README = readlines(joinpath(inputs,"README"))

# Get all the metadata from the file and set up NCvars
ncvars,ncdims,fileatts = readncfile(file_in)

# Rewrite to a file
write(ncvars,file_out,README=README,globalattribs=fileatts)
# -


text:
using Tar, Downloads

PNG_SUITE_DIR = "PngSuite"
PNG_SUITE_PATH = joinpath(PNG_TEST_PATH, PNG_SUITE_DIR)
PNG_SUITE_FILE = joinpath(PNG_TEST_PATH, "PngSuite.tgz")

# See https://github.com/JuliaIO/PNGFiles.jl/pull/5 for discussion of these errors
const imdiff_tolerance = Dict(
    "basi2c16.png" => 0.010046128738460332,
    "basi6a08.png" => 0.021200980392156808,
    "basn2c16.png" => 0.010046128738460332,
    "basn6a08.png" => 0.021200980392156808,
    "bgan6a08.png" => 0.021200980392156808,
    "bgwn6a08.png" => 0.021200980392156808,
    "g03n2c08.png" => 0.05052849264705868,
    "g05n2c08.png" => 0.041383272058823424,
    "g07n2c08.png" => 0.09764476102941164,
    "g25n2c08.png" => 0.203565410539216,
    "oi1n2c16.png" => 0.010046128738460332,
    "oi2n2c16.png" => 0.010046128738460332,
    "oi4n2c16.png" => 0.010046128738460332,
    "oi9n2c16.png" => 0.010046128738460332,
    "pp0n2c16.png" => 0.010046128738460332,
    "ps1n2c16.png" => 0.010046128738460332,
    "ps2n2c16.png" => 0.010046128738460332,
    "tp0n2c08.png" => 0.012859987745098002,
    "tp0n3p08.png" => 0.012331495098039179,
)

if !isdir(PNG_SUITE_PATH)
    mkdir(PNG_SUITE_PATH)
    Downloads.download("https://github.com/JuliaIO/PNG.jl/releases/download/PngSuite-2017jul19/PngSuite-2017jul19.tar.gz", PNG_SUITE_FILE)
    Tar.extract(PNG_SUITE_FILE, PNG_SUITE_PATH)
    rm(PNG_SUITE_FILE)
end

struct _Palleted; end
const pngsuite_colormap = Dict("0g" => Gray, "2c" => RGB, "3p" => _Palleted, "4a" => GrayA, "6a" => RGBA)
# based on values in bgbn4a08.png, bggn4a16.png, bgwn6a08.png, bgyn6a16.png

function parse_pngsuite(x::AbstractString)
    code = splitext(x)[1]
    return (
        case=code[1:end-5],
        is_interlaced=code[end-4]=='i',
        color_type=pngsuite_colormap[code[end-3:end-2]],
        bit_depth=parse(Int, code[end-1:end]),
    )
end
parse_pngsuite(x::Symbol) = parse_pngsuite(String(x))

@testset "PngSuite" begin
    for fpath in glob(joinpath("./**/$(PNG_SUITE_DIR)", "[!x]*[!_new][!_new_im].png"))
        case = last(splitpath(fpath))
        case_info = parse_pngsuite(case)
        C = case_info.color_type
        b = case_info.bit_depth
        is_gray = C in (Gray, GrayA)

        @testset "$(case)" begin
            global read_in_pngf = PNGFiles.load(fpath, gamma = is_gray ? 1.0 : nothing)
            @test typeof(read_in_pngf) <: AbstractMatrix

            path, ext = splitext(fpath)
            newpath = path * "_new" * ext

            open(io->PNGFiles.save(io, read_in_pngf), newpath, "w") #test IO method
            @test PNGFiles.save(newpath, read_in_pngf) === nothing
            global read_in_immag = _standardize_grayness(ImageMagick.load(fpath))

            @testset "$(case): PngSuite/ImageMagick read type equality" begin
                if case_info.case in ("tbb", "tbg", "tbr", "tbw")  # transaprency adds an alpha layer
                    if C === RGB
                        C = RGBA
                    elseif C === Gray
                        C = GrayA
                    end
                end
                if C === _Palleted
                    # _Palleted images could have an alpha channel, but its not evident from "case"
                    @test eltype(read_in_pngf) == eltype(read_in_immag)
                else
                    basictype = b > 8 ? Normed{UInt16,16} : Normed{UInt8,8}
                    @test eltype(read_in_pngf) == C{basictype}
                end
            end
            if b >= 8 # ImageMagick.jl does not read in sub 8 bit images correctly
                @testset "$(case): ImageMagick read values equality" begin
                    imdiff_val = imdiff(collect(read_in_pngf), read_in_immag)
                    onfail(@test imdiff_val <= get(imdiff_tolerance, case, 0.01)) do
                        PNGFiles._inspect_png_read(fpath)
                        _add_debugging_entry(fpath, case, imdiff_val)
                    end
                end
            end
            @testset "$(case): IO is idempotent" begin
                @test all(read_in_pngf .≈ PNGFiles.load(newpath))
                @test all(read_in_pngf .≈ open(io->PNGFiles.load(io), newpath))
            end
        end
    end

    ## TODO: Malformed pngs that should error. This throws `signal (6): Aborted` since we
    ## don't work with `png_jmpbuf` properly.
    # for fpath in glob(joinpath("./**/$(PNG_SUITE_DIR)", "[x]*.png"))
    #     case = splitpath(fpath)[end]
    #     @info case
    #     @testset "$(case)" begin
    #         @test_throws ErrorException PNGFiles.load(fpath)
    #     end
    # end

    @testset "Background chunk" begin
        for fpath in glob(joinpath("./**/$(PNG_SUITE_DIR)", "bg*[!_new][!_new_im].png"))
            case = last(splitpath(fpath))
            bg_color_type = case[3]
            @testset "$(case)" begin
                if bg_color_type == 'a'
                    @test all(PNGFiles.load(fpath, background=true) .≈ PNGFiles.load(fpath, background=false))
                    continue
                end

                path, ext = splitext(fpath)
                if case == "bgbn4a08.png"
                    bg=Gray{N0f8}(0.0)
                elseif case == "bggn4a16.png"
                    bg=Gray{N0f16}(0xab84/typemax(UInt16))
                elseif case == "bgwn6a08.png"
                    bg=RGB{N0f8}(1.0)
                elseif case == "bgyn6a16.png"
                    bg=RGB{N0f16}(1.0, 1.0, 0.0)
                end

                img_bg = PNGFiles.load(fpath, background=true, gamma=1.0)
                img_userbg = PNGFiles.load(fpath, background=bg, gamma=1.0)
                img_nobg = PNGFiles.load(fpath, background=false, gamma=1.0)

                @testset "Background chunk is not ignored" begin
                    @test !(img_bg ≈ img_nobg)
                end

                @testset "File <-> user provided backgrounds" begin
                    @test img_bg ≈ img_userbg
                end

                @testset "Background chunk roundtripping" begin
                    new_path1 = path * "_bg_new1" * ext
                    new_path2 = path * "_bg_new2" * ext
                    PNGFiles.save(new_path1, img_nobg, background=bg, file_gamma=1.0)
                    new1_img_bg = PNGFiles.load(new_path1, background=true, gamma=1.0)
                    new1_img_nobg = PNGFiles.load(new_path1, background=false, gamma=1.0)
                    PNGFiles.save(new_path2, new1_img_nobg, background=bg, file_gamma=1.0)
                    new2_img_bg = PNGFiles.load(new_path2, background=true, gamma=1.0)
                    @test img_bg ≈ new1_img_bg
                    @test new1_img_bg ≈ new2_img_bg
                    Base.rm(new_path1)
                    Base.rm(new_path2)
                end
            end
        end
    end
end


text:
# Factorize two different matrices (one real, one complex) and solve for multiple rhs and free memory
@everywhere using MUMPS
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end
using Compat
include("getDivGrad.jl");

Ar = getDivGrad(30,31,33);
n  = size(Ar,1);
Ac  = Ar + im*spdiagm(rand(n),0);

nrhs = 10;
rhsr = randn(n,nrhs) + im*randn(n,nrhs);
rhsc = randn(n,nrhs);

# println("Factorize real matrix on main instance and destroy")
# tic();
# fact = 0;
# for i=1:100;
# 	fact = factorMUMPS(Ar,1)
# 	destroyMUMPS(fact)
# end
# toc();

println("Factorize real matrix on remote worker and destroy")
tic();
fact = 0;
for i=1:100;
	fact = remotecall_fetch(factorMUMPS,workers()[1], Ar,1)
	remotecall_fetch(destroyMUMPS,workers()[1],fact)
end
toc();

println("DONE!")


text:
using Distributed

# parallelized for loop
@everywhere function spec_loop(spec::SpecParams, disk::DiskParams, Nloop::T; top::Float64=NaN) where T<:Integer
    rms0 = zeros(Nloop)
    avg0 = zeros(Nloop)
    for j in 1:Nloop
        # synthesize spectra
        lambdas, outspec = synthesize_spectra(spec, disk, seed_rng=false, top=top)

        # extract velocities
        v_grid, ccf = calc_ccf(lambdas, outspec, spec, normalize=true)
        rvs, sigs = calc_rvs_from_ccf(v_grid, ccf)

        # compute rms and mean
        rms0[j] = calc_rms(rvs)
        avg0[j] = mean(rvs)
    end
    return mean(avg0), std(avg0), mean(rms0), std(rms0)
end

# small function to check and create file structure
function check_plot_dirs()
    # directories
    dirs = [homedir() * "/grass_output/",
            homedir() * "/grass_output/plots/",
            homedir() * "/grass_output/data/"]

    # create dirs if they dont exist, and return dir names
    for dir in dirs
        if !isdir(dir)
            mkdir(dir)
        end
    end
    return dirs
end

# small function to parse command line arguments
function parse_args(args)
    # chain of if/else to set booleans
    if isempty(args)
        run = true
        plot = true
    else
        if "run" in args
            run = true
        else
            run = false
        end

        if "plot" in args
            plot = true
        else
            plot = false
        end
    end

    # print statements for clarity
    if (run & plot)
        println(">>> Running and plotting " * string(PROGRAM_FILE) * "...")
    elseif run
        println(">>> Running " * string(PROGRAM_FILE) * "...")
    elseif plot
        println(">>> Plotting " * string(PROGRAM_FILE) * "...")
    else
        println(">>> Only defining functions in " * string(PROGRAM_FILE) * "...")
    end
    return run, plot
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
function addearliestharvest!(df)
    df[!,:days_in_garden_until_first_harvest] =  Array{Union{Missing, Day},1}(missing, nrow(df))
    df[!,:earliest_harvest] = Array{Union{Missing, Date},1}(missing, nrow(df))

    for row in eachrow(df)
        if !ismissing(row[:days_in_garden_min]) & !ismissing(row[:period_of_harvest_days])
            row[:days_in_garden_until_first_harvest] = row[:days_in_garden_min] - row[:period_of_harvest_days]
            if !ismissing(row[:date_in])
                row[:earliest_harvest] = Date(row[:date_in] + row[:days_in_garden_until_first_harvest])
            end
        end
    end
end


text:
using DataDeps

register(DataDep(
  "wine",
  "http://archive.ics.uci.edu/ml/datasets/Wine",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data",
  "6be6b1203f3d51df0b553a70e57b8a723cd405683958204f96d23d7cd6aea659",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:14,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
    struct Cone <: Shape

Cone shape model.

"""
struct Cone <: Shape
    id::Int64
    h::Float64
    obj_to_world::Transformation
    world_to_obj::Transformation
end

Cone() = Cone(1, 1.0, Transformation(), Transformation())
Cone(id::Integer) = Cone(id, 1.0, Transformation(), Transformation())
Cone(T::Transformation) = Cone(1, 1.0, T, inv(T))
Cone(id::Integer, T::Transformation) = Cone(id, 1.0, T, inv(T))
Cone(id::Integer, h::Real, T::Transformation) = Cone(id, h, T, inv(T))

function cone_between_points(id::Integer, P1::Point3, r1::Real, P2::Point3, r2::Real)
    if r2 ≈ r1
        error("Cone degenerate into cylinder (TODO: implement)")
    end
    reversed = r2 > r1
    if reversed
        r1, r2 = r2, r1
    end

    s = norm(P2 - P1)
    h = 1.0 - r2/r1
    zscale = s * r1 / (r1 - r2)

    flip = reversed ? translation(Z_AXIS * h) * rotation(X_AXIS, π) : Transformation()
    scale = scaling(r1, r1, zscale)
    rot = rotate_z_to(P2 - P1)
    trans = translation(P1...)

    C = Cone(id, h, trans*rot*scale*flip)

    return C
end

"""
    can_intersect(s::Cone) = true

Is a paraboloid intersectable?

TODO: should this function be replaced by an abstract type like `IntersectableShape`?

"""
can_intersect(s::Cone) = true

"""
    obj_bounds(P::Cone)

The bounding box or a paraboloid in object space.

"""
function obj_bounds(cone::Cone)
    v0 = Point3(-1.0, -1.0, 0.0)
    v1 = Point3(1.0, 1.0, cone.h)
    return BoundingBox(v0, v1)
end

"""
    world_bounds(P::Cone)

The bounding box or a paraboloid transformed to world space.

"""
function world_bounds(cone::Cone)
    cone.obj_to_world(obj_bounds(cone))
end

"""
    (T::Transformation)(P::Cone)

Apply a `Transformation` to a `Cone` object, return a new object.

"""
function (T::Transformation)(P::Cone)
    T2 = T * P.obj_to_world
    Cone(P.id, P.h, T2, inv(T2))
end

function cone_coefs(ray::Ray)
    A = ray.direction.x^2 + ray.direction.y^2  - ray.direction.z^2
    B = 2 * (ray.origin.x * ray.direction.x + ray.origin.y * ray.direction.y - ray.origin.z * ray.direction.z + ray.direction.z)
    C = (ray.origin.x^2 + ray.origin.y^2 - ray.origin.z^2 + 2*ray.origin.z - 1)
    return A, B, C
end

function shape_intersect(r::Ray, cone::Cone)
    ray = cone.world_to_obj(r)

    hit, tnear, tfar = quadratic(cone_coefs(ray)...)

    # Not hit or both hits out of bounds
    if !hit || tnear > r.tmax || tfar < r.tmin || (tnear < r.tmin && tfar > r.tmax)
        return nothing, NaN, NaN
    end

    # Compute hit point, check if it's on our cone
    P1 = ray(tnear)
    P2 = ray(tfar)

    if (0.0 <= P1.z <= cone.h && tnear >= r.tmin)
        P = P1
        t = tnear
    elseif (0.0 <= P2.z <= cone.h && tfar <= r.tmax)
        P = P2
        t = tfar
    else
        return nothing, NaN, NaN
    end

    phi = atan(P.y, P.x)

    n = Normal3(cos(phi)*sqrt(2)/2, sin(phi)*sqrt(2)/2, sqrt(2)/2)
    s = normalize(Vector3(1.0-P.x, 1.0-P.y, 1.0-P.z))

    n = flip_normal(n, ray.direction)

    DG = DifferentialGeometry(cone.obj_to_world(P), cone.obj_to_world(n), cone.obj_to_world(s))

    return DG, t, 5e-4 * t

end

"""
    intersectP(r::Ray, cone::Cone)

"Quick" true/false intersect test between given ray and paraboloid.

"""
function intersectP(r::Ray, cone::Cone)
    ray = cone.world_to_obj(r)

    hit, tnear, tfar = quadratic(cone_coefs(ray)...)

    # Not hit or both hits out of bounds
    if !hit || tnear > r.tmax || tfar < r.tmin || (tnear < r.tmin && tfar > r.tmax)
        return false
    end

    # Compute hit points (on dual cone)
    P1 = ray(tnear)
    P2 = ray(tfar)

    # Return true if either point is on the "actual" cone surface
    return (0.0 <= P1.z <= cone.h && tnear > r.tmin) || (0.0 <= P2.z <= cone.h && tfar < r.tmax)
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
"""
```
t = find_threshold(Balanced(), histogram, edges)
```

Balanced histogram thresholding weighs a histogram and compares the overall weight
of each side of the histogram. Each iteration weight is removed from the heavier
side and the start/middle/end points are recalculated. The agorithm continues
untill the start and end points have converged to meet the middle point.

# Output

Returns a real number `t` that specifies the threshold.

# Details

If after the start and end points have converged the final position is at the
initial start or end point then the histogram must have a single peak and has
failed to find a threshold. In this case the algorithm will fall back to using
the `UnimodalRosin` method to select a threshold.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "cameraman" image in the `TestImages` package.

```julia

using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img, 256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(Balanced(), counts[1:end], edges)
```

# Reference

1. "BI-LEVEL IMAGE THRESHOLDING - A Fast Method", Proceedings of the First International Conference on Bio-inspired Systems and Signal Processing, 2008. Available: [10.5220/0001064300700076](https://doi.org/10.5220/0001064300700076)
"""
function find_threshold(algorithm::Balanced, histogram::AbstractArray, edges::AbstractRange)
    # set initial start/middle/end points and weigths
    Iₛ = 1
    Iₑ = length(histogram)
    Iₘ = round(Int, (Iₛ + Iₑ) / 2)
    Wₗ = 0
    Wᵣ = 0
    for i = 1:Iₘ
        Wₗ += histogram[i]
    end
    for i = Iₘ:Iₑ
        Wᵣ += histogram[i]
    end
    while Iₛ < Iₑ
        if Wₗ < Wᵣ
            Wᵣ -= histogram[Iₑ]
            Iₑ -= 1
            if (Iₛ + Iₑ) / 2 < Iₘ
                Wₗ -= histogram[Iₘ]
                Wᵣ += histogram[Iₘ]
                Iₘ -= 1
            end
        else
            Wₗ -= histogram[Iₛ]
            Iₛ += 1
            if (Iₛ + Iₑ) / 2 > Iₘ
                Wₗ += histogram[Iₘ + 1]
                Wᵣ -= histogram[Iₘ + 1]
                Iₘ += 1
            end
        end
    end
    if Iₘ == 1 || Iₘ == length(histogram)
        @warn "Failed to threshold. Falling back to `UnimodalRosin` method."
        return find_threshold(UnimodalRosin(), histogram, edges)
    else
        return edges[Iₘ]
    end
end


text:
mutable struct Multi
    handle :: Ptr{Cvoid}
    timer  :: Ptr{Cvoid}
end

function Multi()
    timer = jl_malloc(Base._sizeof_uv_timer)
    uv_timer_init(timer)
    multi = Multi(curl_multi_init(), timer)
    finalizer(multi) do multi
        uv_close(multi.timer, cglobal(:jl_free))
        curl_multi_cleanup(multi.handle)
    end
    add_callbacks(multi)
    return multi
end

# libuv callbacks

struct CURLMsg
   msg  :: CURLMSG
   easy :: Ptr{Cvoid}
   code :: CURLcode
end

function check_multi_info(multi::Multi)
    while true
        p = curl_multi_info_read(multi.handle, Ref{Cint}())
        p == C_NULL && return
        message = unsafe_load(convert(Ptr{CURLMsg}, p))
        if message.msg == CURLMSG_DONE
            easy_handle = message.easy
            easy_p_ref = Ref{Ptr{Cvoid}}()
            @check curl_easy_getinfo(easy_handle, CURLINFO_PRIVATE, easy_p_ref)
            easy = unsafe_pointer_to_objref(easy_p_ref[])::Easy
            @assert easy_handle == easy.handle
            easy.code = message.code
            close(easy.progress)
            close(easy.buffers)
        else
            @async @info("unknown CURL message type", msg = message.msg)
        end
    end
end

function event_callback(
    uv_poll_p :: Ptr{Cvoid},
    status    :: Cint,
    events    :: Cint,
)::Cvoid
    ## TODO: use a member access API
    multi_p = unsafe_load(convert(Ptr{Ptr{Cvoid}}, uv_poll_p))
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    sock_p = uv_poll_p + Base._sizeof_uv_poll
    sock = unsafe_load(convert(Ptr{curl_socket_t}, sock_p))
    flags = 0
    events & UV_READABLE != 0 && (flags |= CURL_CSELECT_IN)
    events & UV_WRITABLE != 0 && (flags |= CURL_CSELECT_OUT)
    @check curl_multi_socket_action(multi.handle, sock, flags)
    check_multi_info(multi)
end

function timeout_callback(uv_timer_p::Ptr{Cvoid})::Cvoid
    ## TODO: use a member access API
    multi_p = unsafe_load(convert(Ptr{Ptr{Cvoid}}, uv_timer_p))
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    @check curl_multi_socket_action(multi.handle, CURL_SOCKET_TIMEOUT, 0)
    check_multi_info(multi)
end

# curl callbacks

function timer_callback(
    handle_p   :: Ptr{Cvoid},
    timeout_ms :: Clong,
    multi_p    :: Ptr{Cvoid},
)::Cint
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    @assert handle_p == multi.handle
    if timeout_ms == 0
        @check curl_multi_socket_action(multi.handle, CURL_SOCKET_TIMEOUT, 0)
        check_multi_info(multi)
    elseif timeout_ms > 0
        timeout_cb = @cfunction(timeout_callback, Cvoid, (Ptr{Cvoid},))
        uv_timer_start(multi.timer, timeout_cb, timeout_ms, 0)
    else
        uv_timer_stop(multi.timer)
    end
    return 0
end

function socket_callback(
    easy      :: Ptr{Cvoid},
    sock      :: curl_socket_t,
    action    :: Cint,
    multi_p    :: Ptr{Cvoid},
    uv_poll_p :: Ptr{Cvoid},
)::Cint
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    if action in (CURL_POLL_IN, CURL_POLL_OUT, CURL_POLL_INOUT)
        if uv_poll_p == C_NULL
            uv_poll_p = uv_poll_alloc()
            uv_poll_init(uv_poll_p, sock)
            ## TODO: use a member access API
            unsafe_store!(convert(Ptr{Ptr{Cvoid}}, uv_poll_p), multi_p)
            sock_p = uv_poll_p + Base._sizeof_uv_poll
            unsafe_store!(convert(Ptr{curl_socket_t}, sock_p), sock)
            @check curl_multi_assign(multi.handle, sock, uv_poll_p)
        end
        events = 0
        action != CURL_POLL_IN  && (events |= UV_WRITABLE)
        action != CURL_POLL_OUT && (events |= UV_READABLE)
        event_cb = @cfunction(event_callback, Cvoid, (Ptr{Cvoid}, Cint, Cint))
        uv_poll_start(uv_poll_p, events, event_cb)
    elseif action == CURL_POLL_REMOVE
        if uv_poll_p != C_NULL
            uv_poll_stop(uv_poll_p)
            uv_close(uv_poll_p, cglobal(:jl_free))
            @check curl_multi_assign(multi.handle, sock, C_NULL)
        end
    else
        @async @error("socket_callback: unexpected action", action)
    end
    return 0
end

function add_callbacks(multi::Multi)
    # stash multi handle pointer in timer
    multi_p = pointer_from_objref(multi)
    ## TODO: use a member access API
    unsafe_store!(convert(Ptr{Ptr{Cvoid}}, multi.timer), multi_p)

    # set timer callback
    timer_cb = @cfunction(timer_callback, Cint, (Ptr{Cvoid}, Clong, Ptr{Cvoid}))
    @check curl_multi_setopt(multi.handle, CURLMOPT_TIMERFUNCTION, timer_cb)
    @check curl_multi_setopt(multi.handle, CURLMOPT_TIMERDATA, multi_p)

    # set socket callback
    socket_cb = @cfunction(socket_callback,
        Cint, (Ptr{Cvoid}, curl_socket_t, Cint, Ptr{Cvoid}, Ptr{Cvoid}))
    @check curl_multi_setopt(multi.handle, CURLMOPT_SOCKETFUNCTION, socket_cb)
    @check curl_multi_setopt(multi.handle, CURLMOPT_SOCKETDATA, multi_p)
end


text:
###############
# H-OBCA: Hierarchical Optimization-based Collision Avoidance - a path planner for autonomous parking
# Copyright (C) 2018
# Alexander LINIGER [liniger@control.ee.ethz.ch; Automatic Control Lab, ETH Zurich]
# Xiaojing ZHANG [xiaojing.zhang@berkeley.edu; MPC Lab, UC Berkeley]
# Atsushi SAKAI [atsushisakai@global.komatsu; Komatsu Ltd / MPC Lab]
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
###############
# The paper describing the theory can be found here:
# 	X. Zhang, A. Liniger and F. Borrelli; "Optimization-Based Collision Avoidance"; Technical Report, 2017, [https://arxiv.org/abs/1711.03449]
#   X. Zhang, A. Liniger, A. Sakai and F. Borrelli; "Autonomous  Parking  using  Optimization-Based  Collision  Avoidance"; Technical Report, 2018 [add URL]
###############

###############
# veloSmooth: a velocity smoother
###############

function veloSmooth(v,amax,Ts)
	v_ex  = zeros(length(v)+40,1)
	v_bar = zeros(4,length(v)+40)
	v_bar2 = zeros(4,length(v)+40)
	v_barMM = zeros(1,length(v))

	for i = 1:length(v)
		for j = 1:4
			v_bar[j,i+19] = v[i];
			v_ex[i+19] = v[i];
		end
	end

	v_cut1 = 0.25*abs(v[1])
	v_cut2 = 0.25*abs(v[1])+abs(v[1])

	accPhase = Int(round(abs(v[1])/amax/Ts))

	index1 = find((diff(v_ex).>v_cut1) .& (diff(v_ex).<v_cut2))
	index2 = find(diff(v_ex).>v_cut2)

	index3 = find((diff(v_ex).<-v_cut1) .& (diff(v_ex).>-v_cut2))
	index4 = find(diff(v_ex).<-v_cut2)

	if length(index1) >=1 && index1[1]==19
		index1[1] = index1[1]+1
	end
	if length(index3) >=1 && index3[1]==19
		index3[1] = index3[1]+1
	end

	for j = 1:length(index1)
		if v_ex[index1[j]] > v_cut1 || v_ex[index1[j]+1] > v_cut1
			v_bar[1,index1[j]:index1[j]+accPhase] = linspace(0,abs(v[1]),accPhase+1)''
		elseif v_ex[index1[j]] < -v_cut1 || v_ex[index1[j]+1] < -v_cut1
			v_bar[1,index1[j]-accPhase+1:index1[j]+1] = linspace(-abs(v[1]),0,accPhase+1)''
		end
	end

	for j = 1:length(index3)
		if v_ex[index3[j]] > v_cut1 || v_ex[index3[j]+1] > v_cut1
			v_bar[2,index3[j]-accPhase+1:index3[j]+1] = linspace(abs(v[1]),0,accPhase+1)''
		elseif v_ex[index3[j]] < -v_cut1 || v_ex[index3[j]+1] < -v_cut1
			v_bar[2,index3[j]:index3[j]+accPhase] = linspace(0,-abs(v[1]),accPhase+1)''
		end
	end

	for j = 1:length(index2)
		v_bar[3,index2[j]-accPhase:index2[j]+accPhase] = linspace(-abs(v[1]),abs(v[1]),2*accPhase+1)''
	end

	for j = 1:length(index4)
		v_bar[4,index4[j]-accPhase:index4[j]+accPhase] = linspace(abs(v[1]),-abs(v[1]),2*accPhase+1)''
	end

	for i = 20:length(v)+19
		for j = 1:4
			if v_bar[j,i] == 0
				v_bar2[j,i] = v_bar[j,i]
			elseif sign(v_ex[i]) != sign(v_bar[j,i])
				v_bar2[j,i] = v_ex[i]
			else
				v_bar2[j,i] = v_bar[j,i]
			end
		end
	end

	for i = 20:length(v)+19
		if v_ex[i] > 0
			v_barMM[i-19] = minimum(v_bar2[:,i])
		else
			v_barMM[i-19] = maximum(v_bar2[:,i])
		end
	end

	a = diff(v_barMM')./Ts

	return v_barMM', a

end

function getSmoothProfile(rx,ry,ryaw,Ts,sampleN,motionStep)
	# extract (smooth) velocity profile from Hybrid A* solution
	rv = zeros(length(rx),1)
	for i=1:length(rx)
		if i < length(rx)
			rv[i] = (rx[i+1] - rx[i])/(Ts/sampleN)*cos(ryaw[i]) + (ry[i+1]-ry[i])/(Ts/sampleN)*sin(ryaw[i])
		else
			rv[i] = 0
		end
	end

	# smoothen velocity 0.3 m/s^2 max acceleration
	v,a = veloSmooth(rv,0.3,Ts/sampleN)
	# compute steering angle
	delta = atan.(diff(ryaw)*L/motionStep.*sign.(v[1:end-1]));

	# down-sample
	rx_sampled = rx[1:sampleN:end]
	ry_sampled = ry[1:sampleN:end]
	ryaw_sampled = ryaw[1:sampleN:end]
	v_sampled = v[1:sampleN:end]

	a_sampled = a[1:sampleN:end]
	delta_sampled = delta[1:sampleN:end]

	# construct warm-start solution
	xWS = [rx_sampled ry_sampled ryaw_sampled v_sampled]
	uWS = [delta_sampled a_sampled]
	N = length(rx_sampled)-1

	return xWS, uWS, N

end

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
include("src/App.jl")

function main()
  App.start()
end

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
module Downhill

using LinearAlgebra
using Logging
using Logging: Info, default_logcolor, default_metafmt

export FixedRateDescent, MomentumDescent, NesterovMomentum
export SteepestDescent, HyperGradDescent, CGDescent
export BFGS, CholBFGS

export optimize!, optimize

include("utils.jl")
include("abstract_types.jl")
include("wrappers.jl")
include("linesearch.jl")
include("grad_descent.jl")
include("rate_descent.jl")
include("hypergradient.jl")
include("quasinewton.jl")
include("chol_bfgs.jl")
include("optimize.jl")
include("conjgrad.jl")

end # module


text: struct ELBO <: VariationalObjective end


text:
using PyPlot

include("simulate_exp_gp2.jl")
nt = 10
t = collect(linspace(0,10,nt))
beta = complex(0.05,1.0)
#beta = 0.3+im*0.0
alpha = 1.0+im*0.0
#println(alpha)
#println(beta)
#beta = abs(randn())+im*randn()
#alpha = randn()+im*randn()
c1 = real(alpha)*real(beta)+imag(alpha)*imag(beta)
c2 = real(alpha)*real(beta)-imag(alpha)*imag(beta)
while c1 < 0 || c2 < 0
  beta = 0.01*abs(randn())+im*randn()
  alpha = randn()+im*randn()
  c1 = real(alpha)*real(beta)+imag(alpha)*imag(beta)
  c2 = real(alpha)*real(beta)-imag(alpha)*imag(beta)
end
println(c1," ",c2)
clf()
#for i=1:1
ndev = randn(nt)
data = simulate_exp_gp2(t,alpha,beta,ndev)
#datac = simulate_exp_gp2(t,conj(alpha),conj(beta),ndev)
#datac = simulate_exp_gp2(t,alpha,conj(beta),ndev)
# Now build matrix & Cholesky:
mat = zeros(nt,nt)
for i=1:nt
  for j=1:nt
    dt = abs(t[i]-t[j])
    mat[i,j]=exp(-real(beta)*dt)*(real(alpha)*cos(imag(beta)*dt)+imag(alpha)*sin(imag(beta)*dt))
  end
end
chmat=chol(mat)
sim=*(transpose(chmat),ndev)
# Now compute cholesky decomposition analytically:
chol_anal = zeros(Complex{Float64},nt,nt)
# Loop over columns (see 10/26/16 notes):
gam = 0.0+im*0.0
for i=1:nt
  if i > 1
    gam = sqrt(1.-exp(-2.0*beta*(t[i]-t[i-1])))
  else
    gam = 1.0
  end
# Loop over rows:
  for j=i:nt
    chol_anal[i,j]=sqrt(alpha)*gam*exp(-beta*(t[j]-t[i]))
  end
end
plot(t,0.5*real(data),color="red")
plot(t,sim,color="blue")
plot(t,sim-0.5*real(data),color="green")
#println(abs(alpha)," ",std(real(data)))
#end


text:
using Gadfly
using CSV
using DataFrames
using Random
using Dates
import Cairo, Fontconfig

filepath = ARGS[1]
data = CSV.read(filepath, DataFrame)

nodes = unique(data[!, :node])

inds = randperm(length(nodes))[1:5]

@show somenodes = nodes[inds]

somenodes = [6433, 6126, 7758, 6736, 6726]

filter!(row -> row[:node] in somenodes, data)

p = plot(
	data,
	x=:time,
	y=:flow,
	color=:node,
	Geom.line,
	Guide.xticks(ticks=[DateTime(2021), DateTime(2026), DateTime(2031)]),
	Guide.xlabel("Year"),
	Guide.ylabel("Power flow (MWh)")
)
img = PNG("node_inj.png", 6inch, 4inch)
draw(img, p)


text:
using ExaPF
using LinearAlgebra
using LineSearches
using Printf
using KernelAbstractions
using UnicodePlots
using Statistics

include("algo.jl")

# Augmented Lagrangian method
function auglag(datafile; bfgs=false, iter_max=200, itout_max=1,
                       feasible_start=false)

    # Load problem.
    pf = ExaPF.PowerSystem.PowerNetwork(datafile, 1)
    polar = PolarForm(pf, CPU())

    x0 = ExaPF.initial(polar, State())
    p = ExaPF.initial(polar, Parameters())
    if feasible_start
        prob = run_reduced_ipopt(datafile; hessian=false, feasible=true)
        uk = prob.x
    else
        uk = ExaPF.initial(polar, Control())
    end
    u0 = copy(uk)
    u_start = copy(u0)
    wk = copy(uk)
    u_prev = copy(uk)

    buffer = ExaPF.get(polar, ExaPF.PhysicalState())
    constraints = Function[ExaPF.state_constraint, ExaPF.power_constraints]
    nlp = ExaPF.ReducedSpaceEvaluator(polar, x0, uk, p; constraints=constraints,
                                      ε_tol=1e-10)
    # Init a penalty evaluator with initial penalty c₀
    c0 = 0.1

    aug = ExaPF.AugLagEvaluator(nlp, u0; c₀=c0, scale=true)
    ωtol = 1e-5 #1 / c0

    # initialize arrays
    grad = similar(uk)
    # ut for line search
    ut = similar(uk)
    fill!(grad, 0)
    grad_prev = copy(grad)
    obj_prev = Inf
    norm_grad = Inf

    outer_costs = Float64[]
    cost_history = Float64[]
    grad_history = Float64[]
    feas_history = Float64[]

    ls_algo = BackTracking()
    if bfgs
        H = InverseLBFGSOperator(Float64, length(uk), 50, scaling=true)
        α0 = 1.0
    else
        H = I
        α0 = 1e-4
        α0 = 1.0
    end
    αi = α0
    u♭ = nlp.u_min
    u♯ = nlp.u_max
    ηk = 1.0 / (c0^0.1)

    for i_out in 1:itout_max
        iter = 1
        uk .= u_start
        converged = false
        # Inner iteration: projected gradient algorithm
        # uk, norm_grad, n_iter = projected_gradient(aug, uk; α0=α0)
        uk, norm_grad, n_iter, h_costs, h_grads, h_feas = ngpa(aug, uk; α_bb=α0, α♯=α0)

        # Evaluate current position in the original space
        cons = zeros(ExaPF.n_constraints(nlp))
        ExaPF.constraint!(nlp, cons, uk)
        obj = ExaPF.objective(nlp, uk)
        inf_pr = ExaPF.primal_infeasibility(nlp, cons)
        @printf("#Outer %d %-4d %.3e %.3e \n", i_out, n_iter, obj, inf_pr)

        # History
        push!(outer_costs, obj)
        append!(cost_history, h_costs)
        append!(grad_history, h_grads)
        append!(feas_history, h_feas)

        if (norm_grad < 1e-6) && (inf_pr < 1e-8)
            break
        end

        ρ = 0e-6
        u_start .= ρ * u0 .+ (1 - ρ) .* uk

        # Update the parameters (see Nocedal & Wright, page 521)
        if norm(abs.(aug.infeasibility), Inf) <= ηk
            ExaPF.update_multipliers!(aug)
            ηk = ηk / (aug.ρ^0.9)
        else
            ExaPF.update_penalty!(aug; η=10.0)
            ηk = 1.0 / (aug.ρ^0.1)
        end
    end
    # uncomment to plot cost evolution
    # plt = lineplot(log10.(cost_history), title = "Cost history", width=80);
    # println(plt)

    cons = zeros(ExaPF.n_constraints(nlp))
    ExaPF.constraint!(nlp, cons, uk)
    ExaPF.sanity_check(nlp, uk, cons)
    println("  #objective ", aug.counter.objective)
    println("  #gradient  ", aug.counter.gradient)

    return uk, cost_history, feas_history, grad_history
end

datafile = joinpath(dirname(@__FILE__), "..", "test", "data", "case57.m")

u_opt, ch, fh, gh = @time auglag(datafile; bfgs=false, itout_max=10, feasible_start=false,
                   iter_max=1000)


text:
# Macro for defining object-based backend methods.
#

# Generate single call.
macro blis_ccall_object(funcname, ret, named_param_pairs...)
    # C and Julia function names.
    cfunc  = string("bli_", funcname)
    jlfunc = Symbol("bli_", funcname, "!")
    ccfunc = Expr(:string, cfunc)

    # Separate variable names and their types.
    names  = named_param_pairs[1:2:end]
    params = named_param_pairs[2:2:end]

    # For types, replace xObj with appropriate C and Jl call types.
    # ccparams = ( iType -> eval(:( ( xObj -> $iType )(Ptr{BliObjBase}) )) ).(params)
    # jlparams = ( iType -> eval(:( ( xObj -> $iType )(    BliObj     ) )) ).(params)
    # None-Base types cannot be "verified" with the above `eval`
    #  (will yield "BLIS.ObjeckBackend.BliObj" as variable name hence not found).
    # Do this string (symbol) replacing.
    ccparams = (iSym -> if iSym == :xObj
                :( Ptr{BliObjBase} ) else iSym end).(params)
    jlparams = (iSym -> if iSym == :xObj
                :BliObj else iSym end).(params)

    # Convert C typenames into a static symbol of tuple.
    ccparam_syms = Expr(:tuple, ccparams...)

    # Convert Jl typenames into typed-asserted parameter definition.
    jldefs = ( (iArg, iType) -> :( $iArg::$iType ) ).(names, jlparams)

    # Actual arguments with BliObj type should be
    #  converted to BliObjBase when passed to ccall.
    jlargs = ( (iArg, iType) -> begin
                   if iType == :BliObj
                       :( Ref($iArg.obj) )
                   else
                       iArg
                   end
               end ).(names, jlparams)

    # Docstring.
    funcdoc = Expr(:string, """
        Autogenerated BLIS Object-API function.
        For usage refer the
        [BLIS docs for `$cfunc`.](https://github.com/flame/blis/blob/0.8.0/docs/BLISObjectAPI.md#$funcname)
        """)

    return quote
        $(esc(jlfunc))( $(jldefs...) ) = begin
            ccall(dlsym(libblis, $(esc(ccfunc))), $ret, $ccparam_syms, $(jlargs...))
        end

        @doc $funcdoc $jlfunc
    end
end


text:
#!/usr/bin/env julia

function main()
    f = split.(readlines("input"))

    depth = 0
    position = 0
    for i in range(1, length(f))
        direction = f[i][1]
        amount = parse(Int, f[i][2])
        if direction == "forward"
            position += amount
        elseif direction == "down"
            depth += amount
        elseif direction == "up"
            depth -= amount
        end
    end

    print(depth * position)
end


text:
using HASlib

a = 1
k = 1

hg = HASlib.HexGrid(a)
chan = HASlib.get_channels(hg.l, k, 30, 1)
corr = HASlib.CorrugatedMorse(1,1/7,0.001,hg)
w = HASlib.get_w(0,chan,k,corr)


text:
export Grid

struct Grid
    cells::Array{Cell}
    starts::Array{Int64}
    exits::Array{Int64}
end

function get_empty(cells::Array{Cell}, x::Int64, y::Int64)
    xs = x > 0 ? [x] : shuffle!(collect(1:size(cells, 2)))
    ys = y > 0 ? [y] : shuffle!(collect(1:size(cells, 1)))
    for xi in xs
        for yi in ys
            if object(cells[yi, xi]) == "empty"
                return xi, yi
            end
        end
    end
    error("No empty cell found for $x, $y")
end

function Grid(m::Dict)
    cells = Array{Cell}(undef, m["height"], m["width"])
    for i in eachindex(cells)
        cells[i] = Cell()
    end
    for o in get(m, "objects", [])
        color = get(o, "color", "white")
        obj = get(o, "type", "empty")
        x = get(o, "x", 0)
        y = get(o, "y", 0)
        x, y = get_empty(cells, x, y)
        cells[y, x] = Cell(color, obj)
    end
    inds = LinearIndices(cells)
    starts = Array{Int64}(undef, 0)
    if "starts" in keys(m)
        for s in m["starts"]
            x = get(s, "x", 0)
            y = get(s, "y", 0)
            if (x == 0 || y == 0)
                x, y = get_empty(cells, x, y)
            end
            x, y = get_empty(cells, x, y)
            push!(starts, inds[y, x])
        end
    else
        for i in eachindex(cells)
            if object(cells[i]) != "wall"
                push!(starts, i)
            end
        end
    end
    exits = Array{Int64}(undef, 0)
    if "exits" in keys(m)
        for e in m["exits"]
            x = get(e, "x", 0)
            y = get(e, "y", 0)
            if (x == 0 || y == 0)
                x, y = get_empty(cells, x, y)
            end
            push!(exits, inds[y, x])
        end
    end
    Grid(cells, starts, exits)
end

function Grid(map::String="maps/default.yaml")
    Grid(YAML.load_file(map))
end


text:
# Periodic boundary conditions

module PBC

using LinearAlgebra
using ..Dorothy.Geometry
using ..Dorothy.Properties
using ..Dorothy.Utils
using StaticArrays

export
		TriclinicPBC, RhombododecahedralPBC, OrthorhombicPBC, Kspace, kspace,

		pbcgeometry, pbcmatrix, pbcvolume, rhombododecahedral,
		truncatedoctahedral, hexagonal, orthorhombic, tetragonal, cubic, ispbc,

		TriclinicCell, RhombododecahedralCell, OrthorhombicCell, pbccell,

		wrappos, wrappos!, pbcpos, pbcpos!, eachimage, nearestpos, nearestpos!,
		mindist, sqmindist, UnwrapStrategy, UnwrapByExtent, UnwrapByGap,
		unwrap!

const RealTriple = Tuple{Real,Real,Real}
const PBCGeometry = Tuple{RealTriple,RealTriple}

abstract type TriclinicPBC end

abstract type RhombododecahedralPBC <: TriclinicPBC end

abstract type OrthorhombicPBC <: TriclinicPBC end

struct Kspace end

const kspace = Kspace()

pbcgeometry(cell::TriclinicPBC) = pbcgeometry(pbcmatrix(cell))

pbcgeometry(M::AbstractMatrix{<:Real}) = pbcgeometry(Basis3D(M))

function pbcgeometry(M::Basis3D)
	A, B, C = [M[:,i] for i = 1:3]
	a, b, c = [norm(V) for V in (A, B, C)]
	α, β, γ = acos(B⋅C / (b*c)), acos(A⋅C / (a*c)), acos(A⋅B / (a*b))
	(a,b,c), (α,β,γ)
end

function pbcmatrix((a,b,c)::RealTriple, (α,β,γ)::RealTriple)
	Ax, Ay, Az = a, 0.0, 0.0 # A lies along the positive X axis
	Bx, By, Bz = b*cos(γ), b*sin(γ), 0.0 # B lies in the X-Y plane
	Cx = c*cos(β)
	Cy = (b*c*cos(α) - Bx*Cx) / By
	Cz = √(c^2 - Cx^2 - Cy^2)
	Basis3D([Ax Bx Cx;
	         Ay By Cy;
			 Az Bz Cz])
end

Base.:(==)(cell1::TriclinicPBC, cell2::TriclinicPBC) =
		pbcmatrix(cell1) ≈ pbcmatrix(cell2)

Base.:(*)(cell::TriclinicPBC, R::AbstractVector{<:Real}) = cell.T * R

Base.broadcastable(cell::TriclinicPBC) = Ref(cell)

Base.inv(cell::TriclinicPBC) = cell.inv

Geometry.dims(cell::TriclinicPBC) = Vector3D(sum(pbcmatrix(cell); dims=2))

Geometry.center(cell::TriclinicPBC) = dims(cell) ./ 2

Geometry.extent(cell::TriclinicPBC) =
		BoundingBox(Vector3D(0.0, 0.0, 0.0), dims(cell))

Geometry.isinside(R::AbstractVector{<:Real}, cell::TriclinicPBC) =
		isinside(Vector3D(R), cell)

Geometry.isinside(R::Vector3D, cell::TriclinicPBC) =
		isinside(inv(cell) * R, kspace)

Geometry.isinside(R::Vector3D, cell::OrthorhombicPBC) =
		(0.0 <= R.x < dims(cell).x) && (0.0 <= R.y < dims(cell).y) &&
		(0.0 <= R.z < dims(cell).z)

Geometry.isinside(K::Vector3D, ::Kspace) =
		(0.0 <= K.x < 1.0) && (0.0 <= K.y < 1.0) && (0.0 <= K.z < 1.0)

function pbcvolume((a,b,c)::RealTriple, (α,β,γ)::RealTriple)
	cosα, cosβ, cosγ = cos(α), cos(β), cos(γ)
	a*b*c * √(1 - cosα^2 - cosβ^2 - cosγ^2 + 2*cosα*cosβ*cosγ)
end

pbcvolume(cell::TriclinicPBC) = pbcvolume(pbcmatrix(cell))

pbcvolume(M::AbstractMatrix{<:Real}) = pbcvolume(Basis3D(M))

function pbcvolume(M::Basis3D)
	A, B, C = [M[1:3,i] for i = 1:3]
	A⋅(B×C)
end

Geometry.volume(cell::TriclinicPBC) = pbcvolume(cell)

rhombododecahedral(a::Real) = (a, a, a), (τ/6, τ/6, τ/4)

rhombododecahedral((a,b,c)::RealTriple, (α,β,γ)::RealTriple) =
		rhombododecahedral(a)

function truncatedoctahedral(a::Real)
	θ = arcos(1/3)
	(a, a, a), (θ, τ/2 - θ, θ)
end

truncatedoctahedral((a,b,c)::RealTriple,(α,β,γ)::RealTriple) =
		truncatedoctahedral(a)

hexagonal((a,b,c)::RealTriple, (α,β,γ)::RealTriple) =
		(a, a, c), (τ/4, τ/4, τ/3)

orthorhombic(a::Real, b::Real, c::Real) = (a, b, c), (τ/4, τ/4, τ/4)

orthorhombic((a,b,c)::RealTriple, (α,β,γ)::RealTriple) =
		orthorhombic(a, b, c)

tetragonal(a::Real, c::Real) = (a, a, c), (τ/4, τ/4, τ/4)

tetragonal((a,b,c)::RealTriple,(α,β,γ)::RealTriple) = tetragonal(a, c)

cubic(a::Real) = (a, a, a), (τ/4, τ/4, τ/4)

cubic((a,b,c)::RealTriple, (α,β,γ)::RealTriple) = cubic(a)

pbcrmsd(((a1,b1,c1),(α1,β1,γ1))::PBCGeometry,
		((a2,b2,c2),(α2,β2,γ2))::PBCGeometry) =
		sqrt((abs(a2-a1)^2 + abs(b2-b1)^2 + abs(c2-c1)^2) / 3),
		sqrt((abs(α2-α1)^2 + abs(β2-β1)^2 + abs(γ2-γ1)^2) / 3)

const defpbctol = 0.01, 0.01

function ispbc(((a1,b1,c1),(α1,β1,γ1))::PBCGeometry, geometry;
		tol::Tuple{Real,Real} = defpbctol)
	(a2,b2,c2), (α2,β2,γ2) = geometry((a1,b1,c1), (α1,β1,γ1))
	rmsdedges, rmsdangles =
			pbcrmsd(((a1,b1,c1),(α1,β1,γ1)), ((a2,b2,c2),(α2,β2,γ2)))
	edgetol, angletol = tol
	(rmsdedges <= edgetol) && (rmsdangles <= angletol)
end

ispbc(M::AbstractMatrix{<:Real}, geometry; tol::Tuple{Real,Real} = defpbctol) =
		ispbc(pbcgeometry(M), geometry; tol = tol)

ispbc(cell::TriclinicPBC, geometry; tol::Tuple{Real,Real} = defpbctol) =
		ispbc(pbcgeometry(cell), geometry; tol = tol)

struct TriclinicCell <: TriclinicPBC
	T::LinearTransformation
	inv::LinearTransformation

	TriclinicCell(M::AbstractMatrix{<:Real}) =
			new(LinearTransformation(M), LinearTransformation(inv(M)))
end

TriclinicCell((a,b,c)::RealTriple,(α,β,γ)::RealTriple) =
		TriclinicCell(pbcmatrix((a,b,c), (α,β,γ)))

TriclinicCell(cell::TriclinicPBC) = TriclinicCell(pbcmatrix(cell))

(cell::TriclinicCell)(x) = transform(cell.T, x)

pbcmatrix(cell::TriclinicCell) = cell.T.M

struct RhombododecahedralCell <: RhombododecahedralPBC
	T::LinearTransformation
	inv::LinearTransformation
	r2::Float64

	function RhombododecahedralCell(a::Real)
		M = pbcmatrix((a,a,a), (τ/6, τ/6, τ/4))
		new(LinearTransformation(M), LinearTransformation(inv(M)), (a/2)^2)
	end
end

pbcmatrix(cell::RhombododecahedralCell) = cell.T.M

function pbcgeometry(cell::RhombododecahedralCell)
	a = 2*sqrt(cell.r2)
	(a,a,a), (τ/6,τ/6,τ/4)
end

(cell::RhombododecahedralCell)(x) = transform(cell.T, x)

struct OrthorhombicCell <: OrthorhombicPBC
	T::Scaling
	inv::Scaling

	OrthorhombicCell(V::AbstractVector{<:Real}) =
			new(Scaling(V), Scaling(1.0 ./ V))
end

OrthorhombicCell(x::Real, y::Real, z::Real) = OrthorhombicCell([x, y, z])

OrthorhombicCell(xy::Real, z::Real) = OrthorhombicCell([xy, xy, z])

OrthorhombicCell(xyz::Real) = OrthorhombicCell([xyz, xyz, xyz])

pbcmatrix(cell::OrthorhombicCell) = Diagonal(cell.T.V)

pbcgeometry(cell::OrthorhombicCell) = (dims(cell)...,), (τ/4,τ/4,τ/4)

Geometry.dims(cell::OrthorhombicCell) = cell.T.V

Geometry.volume(cell::OrthorhombicCell) = prod(dims(cell))

(cell::OrthorhombicCell)(x) = transform(cell.T, x)

pbccell(::Nothing; tol::Tuple{Real,Real} = defpbctol) = nothing

pbccell(cell::TriclinicPBC; tol::Tuple{Real,Real} = defpbctol) = cell

pbccell(cell::TriclinicCell; tol::Tuple{Real,Real} = defpbctol) =
		pbccell(pbcgeometry(cell)...; tol = tol)

pbccell(M::AbstractMatrix{<:Real}; tol::Tuple{Real,Real} = defedgetol) =
		pbccell(pbcgeometry(M)...; tol = tol)

function pbccell((a,b,c)::RealTriple, (α,β,γ)::RealTriple;
		tol::Tuple{Real,Real} = defpbctol)
	if ispbc(((a,b,c), (α,β,γ)), rhombododecahedral; tol = tol)
		RhombododecahedralCell(a)
	elseif ispbc(((a,b,c), (α,β,γ)), orthorhombic; tol = tol)
		OrthorhombicCell(a, b, c)
	else
		TriclinicCell((a,b,c), (α,β,γ))
	end
end

wrappos(R::AbstractVector{<:Real}, cell::TriclinicPBC) =
		wrappos(Vector3D(R), cell)

wrappos(R::Vector3D, cell::TriclinicPBC) = cell * wrappos(inv(cell) * R, kspace)

wrappos(K::Vector3D, ::Kspace) = (K - floor.(K))

wrappos!(R::AbstractVector{Vector3D}, cell::TriclinicPBC) =
		wrappos!(cog, R, cell)

function wrappos!(f::Function, R::AbstractVector{Vector3D}, cell::TriclinicPBC)
	@boundscheck length(R) > 0 || error("cannot wrap zero positions")
	C = f(R)
	T = wrappos(C, cell) - C
	if ! iszero(T)
		translation(T)(R)
	end
	R
end

pbcpos(R::AbstractVector{<:Real}, cell::TriclinicPBC) =
		pbcpos(Vector3D(R), cell)

function pbcpos(R::Vector3D, cell::TriclinicPBC)
	Kw = wrappos(inv(cell) * R, kspace)
	Rw = cell * Kw
	Rw, Kw
end

pbcpos(R::Vector3D, ::Nothing) = R, R

pbcpos(R::AbstractVector{Vector3D}, cell::Union{TriclinicPBC,Nothing}) =
		pbcpos!(similar(R), similar(R), R, cell)

function pbcpos!(Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		R::AbstractVector{Vector3D}, cell::Union{TriclinicPBC,Nothing})
	@boundscheck length(Rw) == length(Kw) == length(R) ||
			error("size mismatch between position arrays")
	for i in eachindex(R)
		Rw[i], Kw[i] = pbcpos(R[i], cell)
	end
	Rw, Kw
end

struct PBCImageIterator{T<:TriclinicPBC}
	Kw::Vector3D
	cell::T

	PBCImageIterator{T}(Kw::Vector3D, cell::T) where {T<:TriclinicPBC} =
			new(Kw, cell)
end

PBCImageIterator(Kw::Vector3D, cell::T) where {T<:TriclinicPBC} =
		PBCImageIterator{T}(Kw, cell)

Base.eltype(::Type{PBCImageIterator}) = Vector3D

Base.length(::PBCImageIterator) = 26

const pbcfirstshell = @SVector [
	Vector3D(-1.0,-1.0,-1.0), Vector3D(-1.0,-1.0,0.0), Vector3D(-1.0,-1.0,1.0),
	Vector3D(-1.0,0.0,-1.0),  Vector3D(-1.0,0.0,0.0),  Vector3D(-1.0,0.0,1.0),
	Vector3D(-1.0,1.0,-1.0),  Vector3D(-1.0,1.0,0.0),  Vector3D(-1.0,1.0,1.0),
	Vector3D(0.0,-1.0,-1.0),  Vector3D(0.0,-1.0,0.0),  Vector3D(0.0,-1.0,1.0),
	Vector3D(0.0,0.0,-1.0),                            Vector3D(0.0,0.0,1.0),
	Vector3D(0.0,1.0,-1.0),   Vector3D(0.0,1.0,0.0),   Vector3D(0.0,1.0,1.0),
	Vector3D(1.0,-1.0,-1.0),  Vector3D(1.0,-1.0,0.0),  Vector3D(1.0,-1.0,1.0),
	Vector3D(1.0,0.0,-1.0),   Vector3D(1.0,0.0,0.0),   Vector3D(1.0,0.0,1.0),
	Vector3D(1.0,1.0,-1.0),   Vector3D(1.0,1.0,0.0),   Vector3D(1.0,1.0,1.0)]

function Base.iterate(iter::PBCImageIterator, n = 1)
	if n > 26
		nothing
	else
		iter.cell * (iter.Kw + pbcfirstshell[n]), n + 1
	end
end

eachimage(R::Vector3D, cell::TriclinicPBC) =
		PBCImageIterator(wrappos(inv(cell) * R, kspace), cell)

eachimage(Rw::Vector3D, Kw::Vector3D, cell::TriclinicPBC) =
		PBCImageIterator(Kw, cell)

nearestpos(R::AbstractVector{<:Real}, O::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}) =
		nearestpos(Vector3D(R), Vector3D(O), cell)

nearestpos(R::Vector3D, O::Vector3D, cell::Union{TriclinicPBC,Nothing}) =
		nearestpos(pbcpos(R, cell)..., O, cell)

function nearestpos(Rw::Vector3D, Kw::Vector3D, O::Vector3D, cell::TriclinicPBC)
	Rwmin = Rw
	d2min = sqdist(Rw, O)
	for Rwi in eachimage(Rw, Kw, cell)
		d2 = sqdist(Rwi, O)
		if d2 < d2min
			d2min = d2
			Rwmin = Rwi
		end
	end
	Rwmin
end

function nearestpos(Rw::Vector3D, Kw::Vector3D, O::Vector3D,
		cell::RhombododecahedralPBC)
	Rwmin = Rw
	d2min = sqdist(Rw, O)
	if d2min < cell.r2
		return Rwmin
	end
	for Rwi in eachimage(Rw, Kw, cell)
		d2 = sqdist(Rwi, O)
		if d2 < cell.r2
			return Rwi
		elseif d2 < d2min
			d2min = d2
			Rwmin = Rwi
		end
	end
	Rwmin
end

function nearestpos(Rw::Vector3D, Kw::Vector3D, O::Vector3D,
		cell::OrthorhombicPBC)
	dx, dy, dz = O - Rw
	x, y, z = Rw
	if dx > center(cell).x
		x += dims(cell).x
	elseif dx < -center(cell).x
		x -= dims(cell).x
	end
	if dy > center(cell).y
		y += dims(cell).y
	elseif dy < -center(cell).y
		y -= dims(cell).y
	end
	if dz > center(cell).z
		z += dims(cell).z
	elseif dz < -center(cell).z
		z -= dims(cell).z
	end
	Vector3D(x, y, z)
end

nearestpos(Rw::Vector3D, Kw::Vector3D, O::Vector3D, ::Nothing) = Rw

nearestpos!(R::AbstractVector{Vector3D}, O::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}) = nearestpos!(R, Vector3D(O), cell)

nearestpos!(f, R::AbstractVector{Vector3D}, O::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}) =
		nearestpos!(f, R, Vector3D(O), cell)

nearestpos!(R::AbstractVector{Vector3D}, O::Vector3D,
		cell::Union{TriclinicPBC,Nothing}) = nearestpos!(cog, R, O, cell)

function nearestpos!(f, R::AbstractVector{Vector3D}, O::Vector3D,
		cell::Union{TriclinicPBC,Nothing})
	if ! isempty(R)
		C::Vector3D = f(R)
		T = nearestpos(C, O, cell) - C
		if ! iszero(T)
			translation(T)(R)
		end
	end
	R
end

mindist(R1::AbstractVector{<:Real}, R2::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}) = sqrt(sqmindist(R1, R2, cell))

mindist(R1w::AbstractVector{<:Real}, K1w::AbstractVector{<:Real},
		R2w::AbstractVector{<:Real}, K2w::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}) =
		sqrt(sqmindist(R1w, K1w, R2w, K2w, cell))

sqmindist(R1::AbstractVector{<:Real}, R2::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}) =
		sqmindist(Vector3D(R1), Vector3D(R2), cell)

sqmindist(R1w::AbstractVector{<:Real}, K1w::AbstractVector{<:Real},
		R2w::AbstractVector{<:Real}, K2w::AbstractVector{<:Real},
		cell::Union{TriclinicPBC,Nothing}) = sqmindist(Vector3D(R1w),
		Vector3D(K1w), Vector3D(R2w), Vector3D(K2w), cell)

sqmindist(R1::Vector3D, R2::Vector3D, cell::Union{TriclinicPBC,Nothing}) =
		sqmindist(pbcpos(R1)..., pbcpos(R2)..., cell)

function sqmindist(R1w::Vector3D, K1w::Vector3D, R2w::Vector3D, K2w::Vector3D,
		cell::TriclinicPBC)
	d2min = sqdist(R1w, R2w)
	for R1wi in eachimage(R1w, K1w, cell)
		d2 = sqdist(R1wi, R2w)
		if d2 < d2min
			d2min = d2
		end
	end
	d2min
end

function sqmindist(R1w::Vector3D, K1w::Vector3D, R2w::Vector3D, K2w::Vector3D,
		cell::RhombododecahedralPBC)
	d2min = sqdist(R1w, R2w)
	if d2min < cell.r2
		return d2min
	end
	for R1wi in eachimage(R1w, K1w, cell)
		d2 = sqdist(R1wi, R2w)
		if d2 < cell.r2
			return d2
		elseif d2 < d2min
			d2min = d2
		end
	end
	d2min
end

function sqmindist(R1w::Vector3D, K1w::Vector3D, R2w::Vector3D, K2w::Vector3D,
		cell::OrthorhombicPBC)
	x, y, z = abs.(R2w - R1w)
	if x > center(cell).x
		x = dims(cell).x - x
	end
	if y > center(cell).y
		y = dims(cell).y - y
	end
	if z > center(cell).z
		z = dims(cell).z - z
	end
	x^2 + y^2 + z^2
end

sqmindist(R1w::Vector3D, K1w::Vector3D, R2w::Vector3D, K2w::Vector3D,
		::Nothing) = sqdist(R1w, R2w)

abstract type UnwrapStrategy end

struct UnwrapByExtent <: UnwrapStrategy end

function unwrap!(R::AbstractVector{Vector3D}, cell::TriclinicPBC,
		strategy::UnwrapByExtent)
	length(R) > 1 || return R
	inv(cell)(R)
	unwrap!(R, kspace, strategy)
	cell(R)
end

function unwrap!(K::AbstractVector{Vector3D}, ::Kspace,
		strategy::UnwrapByExtent)
	length(K) > 1 || return K
	Kdims = dims(extent(K))
	if Kdims.x > 0.5
		for i in eachindex(K)
			x, y, z = K[i]
			x -= 0.5
			x -= floor(x)
			x += 0.5
			K[i] = Vector3D(x,y,z)
		end
	end
	if Kdims.y > 0.5
		for i in eachindex(K)
			x, y, z = K[i]
			y -= 0.5
			y -= floor(y)
			y += 0.5
			K[i] = Vector3D(x,y,z)
		end
	end
	if Kdims.z > 0.5
		for i in eachindex(K)
			x, y, z = K[i]
			z -= 0.5
			z -= floor(z)
			z += 0.5
			K[i] = Vector3D(x,y,z)
		end
	end
	K
end

struct UnwrapByGap <: UnwrapStrategy
	Dmax::Vector3D
	buffer::Vector{Vector3D}

	function UnwrapByGap(cell::TriclinicPBC;
			dmax::Real = maximum(values(covalent_radii)))
		(a,b,c), angles = pbcgeometry(cell)
		Dmax = Vector3D(dmax,dmax,dmax) ./ Vector3D(a,b,c)
		new(Dmax, Vector3D[])
	end
end

function unwrap!(R::AbstractVector{Vector3D}, cell::TriclinicPBC,
		strategy::UnwrapByGap)
	length(R) > 1 || return R
	inv(cell)(R)
	unwrap!(R, kspace, strategy)
	cell(R)
end

function unwrap!(K::AbstractVector{Vector3D}, ::Kspace, strategy::UnwrapByGap)
	length(K) > 1 || return K
	Ks = resize!(strategy.buffer, length(K))
	Ks .= K
	sort!(Ks, by = R -> R.x)
	for i = 2:length(Ks)
		if Ks[i].x - Ks[i-1].x > strategy.Dmax.x
			d = Ks[i].x
			for j in eachindex(K)
				x, y, z = K[j]
				x -= d
				x -= floor(x)
				x += d
				K[j] = Vector3D(x,y,z)
			end
			break
		end
	end
	sort!(Ks, by = R -> R.y)
	for i = 2:length(Ks)
		if Ks[i].y - Ks[i-1].y > strategy.Dmax.y
			d = Ks[i].y
			for j in eachindex(K)
				x, y, z = K[j]
				y -= d
				y -= floor(y)
				y += d
				K[j] = Vector3D(x,y,z)
			end
			break
		end
	end
	sort!(Ks, by = R -> R.z)
	for i = 2:length(Ks)
		if Ks[i].z - Ks[i-1].z > strategy.Dmax.z
			d = Ks[i].z
			for j in eachindex(K)
				x, y, z = K[j]
				z -= d
				z -= floor(z)
				z += d
				K[j] = Vector3D(x,y,z)
			end
			break
		end
	end
	K
end

end # module


text:
using SDPA_GMP
using Test

@testset "SDPA_GMP.jl" begin

end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractSelectionRule end

function reset!(::AbstractSelectionRule)
    return nothing
end

"""
    SelectUniform(n::Integer)

Incoming cuts are placed into aggregates uniformly, so that each aggregate has at most `n` cuts. Behaves as [`PartialAggregation`](@ref).

"""
mutable struct SelectUniform <: AbstractSelectionRule
    n::Integer

    function SelectUniform(n::Integer)
        return new(max(1,n))
    end
end

function select(rule::SelectUniform, aggregates::Vector{<:AggregatedOptimalityCut}, ::HyperPlane{OptimalityCut})
    if num_subproblems(aggregates[1]) == rule.n - 1
        return 1, true
    end
    return 1, false
end

function str(::SelectUniform)
    return "uniform selection"
end

"""
    SelectDecaying(T₀::Integer, T̲::Integer = 1, γ::T)

Behaves like [`SelectUniform`](@ref), but the uniform aggregate size decays by `γ` each iteration, starting from `T₀`. `T̲` is an optional lower bound on the aggregate size.

"""
mutable struct SelectDecaying <: AbstractSelectionRule
    n::Int
    T̲::Int
    γ::Float64
    idx::Int

    function SelectDecaying(T₀::Integer, γ::Float64)
        return new(max(1,T₀), 1, γ, 1)
    end

    function SelectDecaying(T₀::Integer, T̲::Integer, γ::Float64)
        return new(max(T₀), max(1,T̲), γ, 1)
    end
end

function select(rule::SelectDecaying, aggregates::Vector{<:AggregatedOptimalityCut}, ::HyperPlane{OptimalityCut})
    if num_subproblems(aggregates[1]) == rule.n - 1
        return 1, true
    end
    return 1, false
end

function reset!(rule::SelectDecaying)
    rule.n = max(rule.T̲, round(Int, rule.γ*rule.n))
    rule.idx = 1
    return nothing
end

function str(::SelectDecaying)
    return "decaying aggregation level selection"
end

"""
    SelectRandom(max = Inf)

Incoming cuts are placed into aggregates randomly. An optional maximum number of cuts `max` can be specified.

"""
struct SelectRandom <: AbstractSelectionRule
    max::Float64
end

SelectRandom(; max = Inf) = SelectRandom(max)

function select(rule::SelectRandom, aggregates::Vector{<:AggregatedOptimalityCut}, ::HyperPlane{OptimalityCut})
    idx = rand(eachindex(aggregates))
    while num_subproblems(aggregates[idx]) > rule.max
        idx = rand(eachindex(aggregates))
    end
    return idx, (num_subproblems(aggregates[idx]) == rule.max-1)
end

function str(rule::SelectRandom)
    if rule.max < Inf
        return "random selection of at most $(rule.max)"
    end
    return "random selection"
end

"""
    SelectClosest(τ::AbstractFloat; distance::Function = absolute_distance)

Incoming cuts are placed into the closest aggregate, according the supplied `distance` function. An empty aggregate is chosen if no aggregate is within the tolerance `τ`

The following distance measures are available
- [`absolute_distance`](@ref)
- [`angular_distance`](@ref)
- [`spatioangular_distance`](@ref

"""
struct SelectClosest <: AbstractSelectionRule
    τ::Float64
    distance::Function
end

function SelectClosest(τ::AbstractFloat; distance::Function = absolute_distance)
    return SelectClosest(τ, distance)
end

function select(rule::SelectClosest, aggregates::Vector{<:AggregatedOptimalityCut}, cut::HyperPlane{OptimalityCut})
    separations = map(c -> rule.distance(c, cut), aggregates)
    (dist, idx) = findmin(separations)
    if dist > rule.τ
        zero_idx = findfirst(iszero, aggregates)
         if zero_idx == nothing
            return idx, false
        end
        return zero_idx, false
    end
    return idx, false
end

function str(::SelectClosest)
    return "distance based selection"
end

"""
    SortByReference(τ::AbstractFloat; distance::Function = absolute_distance)

Incoming cuts are placed into an aggregate based on the distance to a reference cut, according the supplied `distance` function. Behaves as [`SelectClosest`](@ref) if not withing the tolerance `τ` to the reference cut.

The following distance measures are available
- [`absolute_distance`](@ref)
- [`angular_distance`](@ref)
- [`spatioangular_distance`](@ref

"""
mutable struct SortByReference <: AbstractSelectionRule
    distance::Function
    reference::AggregatedOptimalityCut{Float64}
    buffer::AggregatedOptimalityCut{Float64}
end

function SortByReference(; distance::Function = absolute_distance)
    return SortByReference(distance, zero(AggregatedOptimalityCut{Float64}), zero(AggregatedOptimalityCut{Float64}))
end

function select(rule::SortByReference, aggregates::Vector{<:AggregatedOptimalityCut}, cut::HyperPlane{OptimalityCut})
    # Fill buffer with total aggregate each iteration
    rule.buffer += cut
    if iszero(rule.reference)
        # Single-cut before reference has been calculated
        return 1, false
    else
        τs = LinRange(0., 1., length(aggregates))
        dist = rule.distance(cut, agg)
        i = something(findfirst(τ -> dist <= τ, τs), 1)
        return i, false
    end
end

function reset!(rule::SortByReference)
    # Swap in buffer and reset when all cuts have been seen
    rule.reference = zero(AggregatedOptimalityCut{Float64})
    rule.reference += rule.buffer
    rule.buffer = zero(AggregatedOptimalityCut{Float64})
    return nothing
end

function str(::SortByReference)
    return "sort by reference based selection"
end