text:
struct UniformScalingMap{T} <: LinearMap{T}
    λ::T
    M::Int
    function UniformScalingMap(λ::Number, M::Int)
        M < 0 && throw(ArgumentError("size of UniformScalingMap must be non-negative, got $M"))
        return new{typeof(λ)}(λ, M)
    end
end
UniformScalingMap(λ::Number, M::Int, N::Int) =
    (M == N ? UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square"))
UniformScalingMap(λ::T, sz::Dims{2}) where {T} =
    (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square"))

MulStyle(::UniformScalingMap) = FiveArg()

# properties
Base.size(A::UniformScalingMap) = (A.M, A.M)
Base.parent(A::UniformScalingMap) = A.λ
Base.isreal(A::UniformScalingMap) = isreal(A.λ)
LinearAlgebra.issymmetric(::UniformScalingMap) = true
LinearAlgebra.ishermitian(A::UniformScalingMap) = isreal(A)
LinearAlgebra.isposdef(A::UniformScalingMap) = isposdef(A.λ)

# comparison of UniformScalingMap objects
Base.:(==)(A::UniformScalingMap, B::UniformScalingMap) = (A.λ == B.λ && A.M == B.M)

# special transposition behavior
LinearAlgebra.transpose(A::UniformScalingMap) = A
LinearAlgebra.adjoint(A::UniformScalingMap)   = UniformScalingMap(conj(A.λ), size(A))

# multiplication with scalar
Base.:(*)(A::UniformScaling, B::LinearMap) = A.λ * B
Base.:(*)(A::LinearMap, B::UniformScaling) = A * B.λ
Base.:(*)(α::Number, J::UniformScalingMap) = UniformScalingMap(α * J.λ, size(J))
Base.:(*)(J::UniformScalingMap, α::Number) = UniformScalingMap(J.λ * α, size(J))
# needed for disambiguation
Base.:(*)(α::RealOrComplex, J::UniformScalingMap) = UniformScalingMap(α * J.λ, size(J))
Base.:(*)(J::UniformScalingMap, α::RealOrComplex) = UniformScalingMap(J.λ * α, size(J))

# multiplication with vector
Base.:(*)(A::UniformScalingMap, x::AbstractVector) =
    length(x) == A.M ? A.λ * x : throw(DimensionMismatch("*"))

# multiplication with vector/matrix
for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
    @eval begin
        function _unsafe_mul!(y::$outtype, J::UniformScalingMap, x::$intype)
            _scaling!(y, J.λ, x, true, false)
            return y
        end
        function _unsafe_mul!(y::$outtype, J::UniformScalingMap, x::$intype,
                    α::Number, β::Number)
            _scaling!(y, J.λ, x, α, β)
            return y
        end
    end
end

function _scaling!(y, λ, x, α, β)
    if (iszero(α) || iszero(λ))
        iszero(β) && (fill!(y, zero(eltype(y))); return y)
        isone(β) && return y
        # β != 0, 1
        rmul!(y, β)
        return y
    elseif isone(α)
        if iszero(β)
            isone(λ) && return copyto!(y, x)
            y .= λ .* x
            return y
        elseif isone(β)
            isone(λ) && return y .+= x
            y .+= λ .* x
            return y
        else # β != 0, 1
            isone(λ) && (axpby!(one(eltype(x)), x, β, y); return y)
            y .= y .* β .+ λ .* x
            return y
        end
    else # α != 0, 1
        if iszero(β)
            isone(λ) && (y .= x .* α; return y)
            y .= λ .* x .* α
            return y
        elseif isone(β)
            isone(λ) && (axpby!(α, x, β, y); return y)
            y .+= λ .* x .* α
            return y
        else # β != 0, 1
            isone(λ) && (y .= y .* β .+ x .* α; return y)
            y .= y .* β .+ λ .* x .* α
            return y
        end # β-cases
    end # α-cases
end # function _scaling!

# combine LinearMap and UniformScaling objects in linear combinations
Base.:(+)(A₁::LinearMap, A₂::UniformScaling) = A₁ + UniformScalingMap(A₂.λ, size(A₁, 1))
Base.:(+)(A₁::UniformScaling, A₂::LinearMap) = UniformScalingMap(A₁.λ, size(A₂, 1)) + A₂
Base.:(-)(A₁::LinearMap, A₂::UniformScaling) = A₁ - UniformScalingMap(A₂.λ, size(A₁, 1))
Base.:(-)(A₁::UniformScaling, A₂::LinearMap) = UniformScalingMap(A₁.λ, size(A₂, 1)) - A₂


text:
function bt04()
#
# Brugnano and Trigiante
#
    p = poly([1;1;1; -1*[1;1;1;1]; (.5+im)*[1;1;1]; (.5-im)*[1;1;1];
                      .5*(1+im)*[1;1]; .5*(1-im)*[1;1]])

    z = [1 3; -1 4; .5+im 3; .5-im 3; .5*(1+im) 2; .5*(1-im) 2]

    p, PolyZeros(z)
end

ext: jl


text:
# the old

"Gradient ∇Y()"
function gradient(y::RandVar, ω::Ω, vals = linearize(ω))
  # y(ω)
  indomainₛ(y, ω)
  function unpackcall(xs)
    -logerr(indomainₛ(y, unlinearize(xs, ω)))
  end
  ForwardDiff.gradient(unpackcall, vals)
end

function gradient(y::RandVar, sω::SimpleΩ{I, V}, vals) where {I, V <: AbstractArray}
  @assert false
  sω = unlinearize(vals, sω)
  sωtracked = SimpleΩ(Dict(i => param(v) for (i, v) in sω.vals))
  # @grab vals
  l = -err(y(sωtracked))
  # @grab sωtracked
  # @grab y
  # @grab l
  # @assert false
  @assert !(isnan(l))
  Flux.back!(l)
  totalgrad = 0.0
  # @grab sωtracked
  for v in values(sωtracked.vals)
    @assert !(any(isnan(v)))

    @assert !(any(isnan(v.grad)))
    totalgrad += mean(v.grad)
  end
  # @show totalgrad
  sω_ = SimpleΩ(Dict(i => v.data for (i, v) in sωtracked.vals))
  linearize(sω_)
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

@userplot ParetoPlot

@recipe function f(p::ParetoPlot)
  # get user input
  data = p.args[1]

  x = sort(data, rev=true)
  n = length(x)
  logp = log.([i/(n+1) for i in 1:n])
  logx = log.(x)

  seriestype --> :scatter
  xlabel --> "-log(i/(n+1))"
  ylabel --> "log(x*)"

  -logp, logx
end


text:
# using RoMEPlotting
using ZMQ
using Caesar, Caesar.ZmqCaesar
using IncrementalInference, RoME
using Test

## See tests in DFG, IIF, and RoME for more complete coverage of features

@testset "ZMQ Interface" begin
    # Unit tests
    include("multilangzmq/callbackCompatibilityTest.jl")
    # Integration tests
    include("multilangzmq/zmqInternal.jl")
    # End to end tests
    # Do this manually...
    # include("multilangzmq/serverTest.jl")
end

include("testPose2AprilTag4Corner.jl")


text:
function latexinline(x; kwargs...)
    latexstr = latexstring( latexraw(x; kwargs...) )
    COPY_TO_CLIPBOARD && clipboard(latexstr)
    return latexstr
end

latexinline(x::Union{AbstractArray, Tuple}; kwargs...) = [ latexinline(i; kwargs...) for i in x]
latexinline(args...; kwargs...) = latexinline(args; kwargs...)


text:
export subalgebra, decompose, radical, is_central, is_simple

_base_ring(A::AbsAlgAss) = base_ring(A)

################################################################################
#
#  Morphism types
#
################################################################################

morphism_type(::Type{T}, ::Type{S}) where {R, T <: AbsAlgAss{R}, S <: AbsAlgAss{R}} = AbsAlgAssMor{T, S, Generic.MatSpaceElem{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(::Type{T}, ::Type{S}) where {R, T <: AbsAlgAss{gfp_fmpz_elem}, S <: AbsAlgAss{gfp_fmpz_elem}} = AbsAlgAssMor{T, S, gfp_fmpz_mat}

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

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

@doc Markdown.doc"""
    basis(A::AbsAlgAss) -> Vector{AbsAlgAssElem}

Returns the basis of $A$.
"""
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::Vector{elem_type(A)}
end

################################################################################
#
#  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

@attr is_central(A::AbsAlgAss) = dimension_of_center(A) == 1

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

# This is the generic fallback which constructs an associative algebra
@doc Markdown.doc"""
    subalgebra(A::AbsAlgAss, e::AbsAlgAssElem, idempotent::Bool = false,
               action::Symbol = :left)
      -> AlgAss, AbsAlgAssMor

Given an algebra $A$ and an element $e$, this function constructs the algebra
$e \cdot A$ (if `action == :left`) respectively $A \cdot e$ (if `action == :right`)
and a map from this algebra to $A$.
If `idempotent` is `true`, it is assumed that $e$ is idempotent in $A$.
"""
function subalgebra(A::AbsAlgAss{T}, e::AbsAlgAssElem{T}, idempotent::Bool = false, action::Symbol = :left) where {T}
  @assert parent(e) == A
  B, mB = AlgAss(A)
  C, mC = subalgebra(B, mB\e, idempotent, action)
  mD = compose_and_squash(mB, mC)
  @assert domain(mD) == C
  return C, mD
end

@doc Markdown.doc"""
    subalgebra(A::AbsAlgAss, basis::Vector{AbsAlgAssElem})
      -> AlgAss, AbsAlgAssMor

Returns the subalgebra $A$ generated by the elements in `basis` and a map
from this algebra to $A$.
"""
function subalgebra(A::AbsAlgAss{T}, basis::Vector{ <: AbsAlgAssElem{T} }) where T
  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) -> Array{Tuple{AlgAss, AbsAlgAssMor}}

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

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

  res = _decompose(A)
  A.decomposition = res
  return res
end

# Generic function for everything besides AlgAss
function decompose(A::AbsAlgAss{T}) where T
  return __decompose(A)
end

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

  B, mB = AlgAss(A)

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

  D = _decompose(B)::Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, AlgAss{T})}}
  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
  Z, ZtoB = center(B)
  if dim(Z) != dim(B)
    if isdefined(A, :center)
      @assert A.center[1] === Z
    end
    A.center = (Z, compose_and_squash(mB, ZtoB))
  end
  return res
end

function _decompose(A::AbsAlgAss{T}) where {T}
  @assert _issemisimple(A) != 2 "Algebra is not semisimple"
  if is_commutative(A)
    return _dec_com(A)::Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}
  else
    return _dec_via_center(A)::Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}
  end
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].is_simple = 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{T}) where {T}
  if characteristic(base_ring(A)) > 0
    return _dec_com_finite(A)::Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}
  else
    return _dec_com_gen(A)::Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}}
  end
end

function _dec_com_gen(A::AbsAlgAss{T}) where {T <: FieldElem}
  if dim(A) == 1
    A.is_simple = 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 is_irreducible(f)
      if degree(f) == dim(A)
        A.is_simple = 1
        B, mB = AlgAss(A)
        return Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}[(B, mB)]
      end
      continue
    end

    @assert is_squarefree(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.is_simple = 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)::Vector{Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(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.is_simple = 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.is_simple = 1
    B, mB = AlgAss(A)
    return Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}[(B, mB)]
  end

  A.is_simple = 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 is_squarefree(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"""
    components(::Type{Field}, A::AbsAlgAss)
      -> Vector{Tuple{Field, Morphism}}

Given an étale algebra $A$, return the simple components of $A$
as fields $K$ together with the projection $A \to K$.
"""
function components(::Type{Field}, A::AbsAlgAss)
  @assert iscommutative(A)
  return as_number_fields(A)
end

@doc Markdown.doc"""
    component(::Type{Field}, A::AbsAlgAss, i::Int)
      -> Vector{Tuple{Field, Morphism}}

Given an étale algebra $A$ and index $i$, return the $i$-th simple components
of $A$ as a field $K$ together with the projection $A \to K$.
"""
function component(::Type{Field}, A::AbsAlgAss, i::Int)
  nf = as_number_fields(A)
  return nf[i]
end

@doc Markdown.doc"""
    as_number_fields(A::AbsAlgAss{fmpq})
      -> Vector{Tuple{AnticNumberField, AbsAlgAssToNfAbsMor}}

Given a commutative algebra $A$ over $\mathbb Q$, this function returns a
decomposition of $A$ as direct sum of number fields and maps from $A$ to
these fields.
"""
function as_number_fields(A::AbsAlgAss{T}) where {T}
  return __as_number_fields(A)
end

function __as_number_fields(A::AbsAlgAss{T}; use_maximal_order::Bool = true) where {T}
  if isdefined(A, :maps_to_numberfields)
    NF = A.maps_to_numberfields::Vector{Tuple{_ext_type(T), _abs_alg_ass_to_nf_abs_mor_type(A)}}
    return NF
  end

  result = _as_number_fields(A, use_maximal_order = use_maximal_order)
  @assert all(domain(AtoK) === A for (_, AtoK) in result)
  A.maps_to_numberfields = result
  return result
end

_ext_type(::Type{fmpq}) = AnticNumberField

_ext_type(::Type{nf_elem}) = NfRel{nf_elem}

function _as_number_fields(A::AbsAlgAss{T}; use_maximal_order::Bool = true) where {T}
  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 && T === fmpq && use_maximal_order
    # 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_matrix(OA, copy = false).num)
    n = basis_matrix(OA, copy = false).den
    basis_lll = elem_type(A)[ elem_from_mat_row(A, L, i, n) for i = 1:d ]
  elseif fields_not_cached
    basis_lll = basis(A)
  end

  KK = base_ring(A)

  M = zero_matrix(KK, 0, d)
  matrices = Vector{dense_matrix_type(T)}()
  fields = Vector{_ext_type(T)}()
  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{_ext_type(T), _abs_alg_ass_to_nf_abs_mor_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]::Tuple{_ext_type(T), _abs_alg_ass_to_nf_abs_mor_type(B)}
      push!(fields, K)
    end

    # TODO: We can do a short cut, but the map must be BtoK \circ inv(BtoA)
    #if length(Adec) == 1
    #  res = Tuple{_ext_type(T), _abs_alg_ass_to_nf_abs_mor_type(A)}[(K, BtoK)]
    #  A.maps_to_numberfields = res
    #  return res
    #end

    # Construct the map from K to A
    N = zero_matrix(KK, 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{dense_matrix_type(T)}(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{_ext_type(T), _abs_alg_ass_to_nf_abs_mor_type(A)}}()
  for i = 1:length(fields)
    push!(result, (fields[i], AbsAlgAssToNfAbsMor(A, fields[i], matrices2[i], matrices[i])))
  end

  return result
end

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

Random.gentype(::Type{T}) where {T<:AbsAlgAss} = elem_type(T)

function rand(rng::AbstractRNG, Asp::Random.SamplerTrivial{<:AbsAlgAss{T}}) where T
  A = Asp[]
  c = rand(rng, base_ring(A), dim(A))
  return A(c)
end

function rand(A::AbsAlgAss{nf_elem}, rng::UnitRange{Int} = -10:10)
  c = nf_elem[rand(base_ring(A), rng) 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
#
################################################################################

# Reduces the vector v w. r. t. M and writes it in the i-th row of M.
# M should look like this:
#     (0 1 * 0 *)
#     (1 0 * 0 *)
# M = (0 0 0 1 *)
#     (0 0 0 0 0)
#     (0 0 0 0 0),
# i. e. "almost" in rref, but the rows do not have to be sorted.
# For a column c of M pivot_rows[c] should be the row with the pivot or 0.
# The function changes M, v and pivot_rows in place!
function _add_row_to_rref!(M::MatElem{T}, v::Vector{T}, pivot_rows::Vector{Int}, i::Int) where { T <: FieldElem }
  @assert ncols(M) == length(v)
  @assert ncols(M) == length(pivot_rows)
  @assert 1 <= i && i <= nrows(M)

  pivot_found = false
  pivot_col = 0
  s = one(base_ring(M))
  for c = 1:ncols(M)
    if iszero(v[c])
      continue
    end
    if pivot_rows[c] == 0
      # We found an entry in a column of v, where no other row of M has the pivot.
      if pivot_found
        # We already found a pivot
        continue
      end

      pivot_found = true
      pivot_col = c
      pivot_rows[c] = i
      continue
    end

    r = pivot_rows[c]
    # Reduce the entries of v by the row r of M
    t = -v[c] # we assume M[r, c] == 1 (it is the pivot)
    for j = c + 1:ncols(M)
      Mrj = M[r, j]
      if iszero(Mrj)
        continue
      end

      s = mul!(s, t, Mrj)
      v[j] = addeq!(v[j], s)
    end
    v[c] = zero!(v[c])
  end
  if !pivot_found
    return false
  end

  # Make the pivot 1
  t = inv(v[pivot_col])
  for j = pivot_col + 1:ncols(M)
    if iszero(v[j])
      continue
    end

    v[j] = mul!(v[j], v[j], t)
  end
  v[pivot_col] = one(base_ring(M))

  # Reduce the rows above the newly inserted one
  for r = 1:i - 1
    Mrp = M[r, pivot_col]
    if iszero(Mrp)
      continue
    end

    t = -Mrp
    for c = pivot_col + 1:ncols(M)
      s = mul!(s, t, v[c])
      M[r, c] = addeq!(M[r, c], s)
    end
    M[r, pivot_col] = zero(base_ring(M))
  end

  for c = 1:ncols(M)
    M[i, c] = deepcopy(v[c])
  end

  return true
end

@doc Markdown.doc"""
    gens(A::AbsAlgAss, return_full_basis::Typel{Val{T}} = Val{false};
         thorough_search::Bool = false) where T
      -> Vector{AbsAlgAssElem}

Returns a subset of `basis(A)`, which generates $A$ as an algebra over
`base_ring(A)`.
If `return_full_basis` is set to `Val{true}`, the function also returns a
`Vector{AbsAlgAssElem}` containing a full basis consisting of monomials in
the generators and a `Vector{Vector{Tuple{Int, Int}}}` containing the
information on how these monomials are built. E. g.: If the function returns
`g`, `full_basis` and `v`, then we have
`full_basis[i] = prod( g[j]^k for (j, k) in v[i] )`.
If `thorough_search` is `true`, the number of returned generators is possibly
smaller. This will in general increase the runtime. It is not guaranteed that
the number of generators is minimal in any case.
"""
function gens(A::AbsAlgAss, return_full_basis::Type{Val{T}} = Val{false}; thorough_search::Bool = false) where T
  d = dim(A)
  if return_full_basis === Val{false}
    if isdefined(A, :gens)
      return A.gens::Vector{elem_type(A)}
    end
  end

  if thorough_search
    # 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)
  else
    is_gen = falses(d)
  end

  generators = Vector{elem_type(A)}()
  full_basis = elem_type(A)[ one(A) ] # Contains products of generators which form a full basis
  elts_in_gens = Vector{Tuple{Int, Int}}[ Tuple{Int, Int}[] ]
  B = zero_matrix(base_ring(A), d, d)
  pivot_rows = zeros(Int, d)
  new_elements = Set{Int}()

  s = one(A)
  t = one(A)

  cur_dim = 0
  cur_basis_elt = 1
  while cur_dim != d
    if isempty(new_elements)
      # We have to add a generator
      new_gen = A()
      new_elt = false
      while true
        if thorough_search
          i = minpoly_degrees[cur_basis_elt][1]
        else
          i = rand(1:dim(A))
          while is_gen[i]
            i = rand(1:dim(A))
          end
          is_gen[i] = true
        end
        new_gen = A[i]
        cur_basis_elt += 1
        new_elt = _add_row_to_rref!(B, coefficients(new_gen), pivot_rows, cur_dim + 1)
        if new_elt
          break
        end
      end
      push!(generators, new_gen)
      b = new_gen
      power = 1
      # Compute the powers of new_gen
      while new_elt
        cur_dim += 1
        push!(full_basis, b)
        if power == 1 && length(generators) != 1
          push!(new_elements, length(full_basis))
        end
        ind = Tuple{Int, Int}[ (length(generators), power) ]
        push!(elts_in_gens, ind)
        cur_dim == d ? break : nothing
        b *= new_gen
        power += 1
        new_elt = _add_row_to_rref!(B, coefficients(b), pivot_rows, cur_dim + 1)
      end
      continue
    else
      i = pop!(new_elements)
      b = full_basis[i]
    end

    # Compute all possible products involving b
    n = length(full_basis)
    for r = 1:n
      s = mul!(s, b, full_basis[r])
      for l = 1:n
        if !is_commutative(A)
          t = mul!(t, full_basis[l], s)
        else
          t = s
        end
        new_elt = _add_row_to_rref!(B, coefficients(t), pivot_rows, cur_dim + 1)
        if !new_elt
          continue
        end
        push!(full_basis, deepcopy(t))
        cur_dim += 1
        coord = _merge_elts_in_gens!(elts_in_gens[l], deepcopy(elts_in_gens[i]), elts_in_gens[r])
        push!(elts_in_gens, coord)
        if thorough_search && coord[1][2] == 1 && coord[end][2] == 1
          push!(new_elements, length(full_basis))
        end
        if is_commutative(A)
          break
        end
        cur_dim == d ? break : nothing
      end
      cur_dim == d ? break : nothing
    end
  end

  # Remove the one
  popfirst!(full_basis)
  popfirst!(elts_in_gens)

  if !isdefined(A, :gens)
    A.gens = generators
  end

  if return_full_basis == Val{true}
    return generators, full_basis, elts_in_gens
  else
    return generators
  end
end

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

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

function primitive_element(A::AbsAlgAss{fmpq})
  if isdefined(A, :maps_to_numberfields)
    return primitive_element_via_number_fields(A)
  end
  a, _ = _primitive_element(A)
  return a
end

# TODO: Fix this with the types once a new version is released
#function _primitive_element(A::AbsAlgAss)
#  error("Not implemented yet")
#  return nothing
#end

# If T == fmpq, we try to find a small primitive element by
# going "via number fields". There a procedure using LLL
# is implemented to find primitive elements with small minimal
# polynomial. Note that this could be improved by calling into
# simplify for number fields. But it is a bit tricky.
function _primitive_element(A::AbsAlgAss{fmpq})
  a = primitive_element_via_number_fields(A)
  return a, minpoly(a)
end

function __primitive_element(A::S) where {T <: FinFieldElem, S <: AbsAlgAss{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 primitive_element_via_number_fields(A::AbsAlgAss{fmpq})
  fields_and_maps = as_number_fields(A)
  a = A()
  for (K, AtoK) in fields_and_maps
    a += AtoK\gen(K)
  end
  return a
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{fmpq}) #<: Union{fmpq, gfp_elem, Generic.ResF{fmpz}, fq_nmod, fq} }
  return _as_field_with_isomorphism(A, _primitive_element(A)...)
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, 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

  return __as_field_with_isomorphism(A, mina, M)
end

function __as_field_with_isomorphism(A::AbsAlgAss{S}, f, M) where {S <: Union{fmpq, nf_elem}}
  K = number_field(f, cached = false)[1]
  return K, AbsAlgAssToNfAbsMor(A, K, inv(M), M)
end

function __as_field_with_isomorphism(A::AbsAlgAss{gfp_elem}, f::gfp_poly, M::gfp_mat)
  Fq = FqNmodFiniteField(f, Symbol("a"), false)
  return Fq, AbsAlgAssToFqMor(A, Fq, inv(M), M, parent(f))
end

function __as_field_with_isomorphism(A, f::gfp_fmpz_poly, M)
#function __as_field_with_isomorphism(A::AbsAlgAss{Generic.ResF{fmpz}}, f::gfp_fmpz_poly, M::Generic.MatSpaceElem{Generic.ResF{fmpz}})
  Fq = FqFiniteField(f, Symbol("a"), false)
  return Fq, AbsAlgAssToFqMor(A, Fq, inv(M), M, parent(f))
end

function __as_field_with_isomorphism(A::AbsAlgAss{S}, f::T, M::U) where { S <: Union{ fq_nmod, fq }, T <: Union{ fq_nmod_poly, fq_poly }, U <: Union{ fq_nmod_mat, fq_mat } }
  Fr, RtoFr = field_extension(f)
  return Fr, AbsAlgAssToFqMor(A, Fr, inv(M), M, parent(f), RtoFr)
end

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

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

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 = matrix_algebra(K, [ representation_matrix(A[i], :right) for i = 1:dim(A) ], isbasis = true)
  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::Vector{T}, el::T) where T
  i=1
  while true
    if image_primitive_element(el) == image_primitive_element(G[i])
      return i
    end
    i+=1
  end
end

#K/Q is a Galois extension.
function CrossedProductAlgebra(K::AnticNumberField, G::Vector{T}, cocval::Matrix{nf_elem}) 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::Vector{T}, cocval::Matrix{nf_elem}) 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.is_simple = 1
  return A

end

###############################################################################
#
#  Trace Matrix
#
###############################################################################

function _assure_trace_basis(A::AbsAlgAss{T}) where T
  if isdefined(A, :trace_basis_elem)
    return nothing
  end

  mt = multiplication_table(A, copy = false)
  A.trace_basis_elem = Vector{T}(undef, dim(A))
  for i = 1:dim(A)
    A.trace_basis_elem[i] = sum(mt[i, j, j] for j = 1:dim(A))
  end
  return nothing
end

@doc Markdown.doc"""
    trace_matrix(A::AbsAlgAss) -> MatElem

Returns a matrix $M$ over the base ring of $A$ such that
$M_{i, j} = \mathrm{tr}(b_i \cdot b_j)$, where $b_1, \dots, b_n$ is the
basis of $A$.
"""
function trace_matrix(A::AbsAlgAss)
  F = base_ring(A)
  n = dim(A)
  M = zero_matrix(F, n, n)
  for i = 1:n
    M[i, i] = tr(A[i]^2)
  end
  for i = 1:n
    for j = i + 1:n
      x = tr(A[i]*A[j])
      M[i, j] = x
      M[j, i] = x
    end
  end
  return M
end

################################################################################
#
#  Change of ring
#
################################################################################

# This is the generic fallback which constructs an associative algebra
@doc Markdown.doc"""
    restrict_scalars(A::AbsAlgAss{nf_elem}, Q::FlintRationalField)
    restrict_scalars(A::AbsAlgAss{fq_nmod}, Fp::GaloisField)
    restrict_scalars(A::AbsAlgAss{fq}, Fp::Generic.ResField{fmpz})
      -> AlgAss, Function, Function

Given an algebra $A$ over a field $L$ and the prime field $K$ of $L$, this
function returns the restriction $B$ of $A$ to $K$ and maps from $A$ to $B$
and from $B$ to $A$.
"""
function restrict_scalars(A::AbsAlgAss{T}, K::Field) where {T}
  #K == base_ring(A) && throw(error("Not yet implemented"))
  B, BtoA = AlgAss(A)::Tuple{AlgAss{T}, morphism_type(AlgAss{T}, typeof(A))}
  C, BtoC, CtoB = restrict_scalars(B, K)

  function CtoA(x::AbsAlgAssElem)
    @assert parent(x) === C
    return BtoA(CtoB(x))
  end

  function AtoC(x::AbsAlgAssElem)
    @assert parent(x) === A
    return BtoC(BtoA\x)
  end

  return C, AtoC, CtoA
end

################################################################################
#
#  Radical
#
################################################################################

@doc Markdown.doc"""
     radical(A::AbsAlgAss) -> AbsAlgAssIdl

Returns the Jacobson-Radical of $A$.
"""
function radical(A::AbsAlgAss{T}) where { T } #<: Union{ gfp_elem, Generic.ResF{fmpz}, fq_nmod, fq, fmpq, nf_elem } }
  return ideal_from_gens(A, _radical(A), :twosided)
end

# Section 2.3.2 in W. Eberly: Computations for Algebras and Group Representations
# TODO: Fix the type
function _radical(A::AbsAlgAss{T}) where { T } #<: Union{ gfp_elem, Generic.ResF{fmpz} } }
  F = base_ring(A)
  p = characteristic(F)
  k = flog(fmpz(dim(A)), p)

  MF = trace_matrix(A)
  d, B = nullspace(MF)
  if d == 0
    return elem_type(A)[]
  end

  C = transpose(B)
  # Now, iterate: we need to find the kernel of tr((xy)^(p^i))/p^i mod p
  # on the subspace generated by C
  # Hard to believe, but this is linear!!!!
  MZ = zero_matrix(FlintZZ, dim(A), dim(A))
  pl = fmpz(1)
  a = A()
  for l = 1:k
    pl = p*pl
    M = zero_matrix(F, dim(A), nrows(C))
    for i = 1:nrows(C)
      c = elem_from_mat_row(A, C, i)
      for j = 1:dim(A)
        a = mul!(a, c, A[j])
        MF = representation_matrix(a)
        for m = 1:nrows(MF)
          for n = 1:ncols(MF)
            MZ[m, n] = lift(MF[m, n])
          end
        end
        t = tr(MZ^Int(pl))
        @assert iszero(mod(t, pl))
        M[j, i] = F(divexact(t, pl))
      end
    end
    d, B = nullspace(M)
    if d == 0
      return elem_type(A)[]
    end
    C = transpose(B)*C
  end

  return elem_type(A)[ elem_from_mat_row(A, C, i) for i = 1:nrows(C) ]
end

function _radical(A::AbsAlgAss{T}) where { T <: Union{ fq_nmod, fq } }
  F = base_ring(A)

  p = characteristic(F)
  if T <: fq_nmod
    Fp = GF(Int(p))
  else
    Fp = GF(p)
  end
  A2, A2toA = restrict_scalars(A, Fp)
  n = degree(F)

  if n == 1
    J = _radical(A2)
    return elem_type(A)[ A2toA(b) for b in J ]
  end

  k = flog(fmpz(dim(A)), p)
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = Qx(push!(fmpq[ -fmpq(coeff(gen(F)^n, i)) for i = 0:(n - 1) ], fmpq(1)))
  K, a = number_field(f, "a")

  MF = trace_matrix(A2)
  d, B = nullspace(MF)
  if d == 0
    return elem_type(A)[]
  end

  C = transpose(B)
  pl = fmpz(1)
  MK = zero_matrix(K, dim(A), dim(A))
  MQx = zero_matrix(Qx, dim(A), dim(A))
  a = A()
  for l = 1:k
    pl = p*pl
    M = zero_matrix(Fp, dim(A2), nrows(C))
    for i = 1:nrows(C)
      c = A2toA(elem_from_mat_row(A2, C, i))
      for j = 1:dim(A)
        a = mul!(a, c, A[j])
        MF = representation_matrix(a)
        MK = _lift_fq_mat!(MF, MK, MQx)
        t = tr(MK^Int(pl))
        @assert all([ iszero(mod(coeff(t, s), pl)) for s = 0:(n - 1) ])
        jn = (j - 1)*n
        for s = 1:n
          M[jn + s, i] = Fp(divexact(numerator(coeff(t, s - 1)), pl))
        end
      end
    end
    d, B = nullspace(M)
    if d == 0
      return elem_type(A)[]
    end
    C = transpose(B)*C
  end

  return elem_type(A)[ A2toA(elem_from_mat_row(A2, C, i)) for i = 1:nrows(C) ]
end

function _lift_fq_mat!(M1::MatElem{T}, M2::MatElem{nf_elem}, M3::MatElem{fmpq_poly}) where { T <: Union{ fq_nmod, fq } }
  @assert ncols(M1) == ncols(M2) && ncols(M1) == ncols(M3)
  @assert nrows(M1) == nrows(M2) && nrows(M1) == nrows(M3)
  n = degree(base_ring(M1))
  K = base_ring(M2)
  R = base_ring(M3)
  for i = 1:nrows(M1)
    for j = 1:ncols(M1)
      # Sadly, there is no setcoeff! for nf_elem...
      for k = 0:(n - 1)
        M3[i, j] = setcoeff!(M3[i, j], k, fmpq(coeff(M1[i, j], k)))
      end
      ccall((:nf_elem_set_fmpq_poly, libantic), Nothing, (Ref{nf_elem}, Ref{fmpq_poly}, Ref{AnticNumberField}), M2[i, j], M3[i, j], K)
    end
  end
  return M2
end

function _radical(A::AbsAlgAss{T}) where { T <: Union{ fmpq, NumFieldElem } }
  M = trace_matrix(A)
  n, N = nullspace(M)
  b = Vector{elem_type(A)}(undef, n)
  t = zeros(base_ring(A), dim(A))
  # the construct A(t) will make a copy (hopefully :))
  for i = 1:n
    for j = 1:dim(A)
      t[j] = N[j, i]
    end
    b[i] = A(t)
  end
  return b
end

################################################################################
#
#  is_simple
#
################################################################################

# TODO: Fix this once we have a new Nemo version

# TODO: Make sure this always returns 1 or 2. So far we only have _radical for
# algebras over Fp, Fq, QQ, and number fields
#function _issemisimple(A::AbsAlgAss)
#  return A.issemisimple
#end

function _issemisimple(A::AbsAlgAss{T}) where { T } #<: Union{ gfp_elem, Generic.ResF{fmpz}, fq, fq_nmod, fmpq, nf_elem } }
  if A.issemisimple == 0
    if isempty(_radical(A))
      A.issemisimple = 1
    else
      A.issemisimple = 2
    end
  end

  return A.issemisimple
end

is_simple_known(A::AbsAlgAss) = A.is_simple != 0

function is_simple(A::AbsAlgAss)
  if A.is_simple != 0
    return A.is_simple == 1
  end

  if _issemisimple(A) == 2
    A.is_simple = 2
    return false
  end
  # Still can't be certain that A is semisimple, since _issemisimple does not
  # always work...

  Adec = decompose(A)
  if length(Adec) == 1
    A.is_simple = 1
    return true
  end
  A.is_simple = 2
  return false
end

################################################################################
#
#  Trace signature
#
################################################################################

function trace_signature(A::AbsAlgAss{nf_elem}, P::InfPlc)
  M = trred_matrix(basis(A))
  Ky, y = PolynomialRing(base_ring(A), "y", cached = false)
  f = charpoly(Ky, M)
  npos = n_positive_roots(f, P; multiplicities = true)
  return (npos, degree(f) - npos)
end

function trace_signature(A::AbsAlgAss{fmpq})
  O = get_attribute(A, :any_order)
  if O === nothing
    M = trred_matrix(basis(A))
  else
    _M = trred_matrix(O::order_type(A))
    M = change_base_ring(QQ, _M)
  end

  Ky, y = PolynomialRing(base_ring(A), "y", cached = false)
  f = charpoly(Ky, M)
  npos = n_positive_roots(f, multiplicities = true)
  return npos, degree(f) - npos
end


text:
"""
    metric(f::Function; name=nothing) -> Function

Converts a julia function `f` to a python metric function to be used by Keras.
Assumes that `f` is of the form `f(::Tensor, ::Tensor) -> Tensor`
"""
@pydef mutable struct metric
    function __init__(self, f::Function, name=nothing)
        self[:f] = f
        self[:__name__] = name == nothing ? typeof(f).name.mt.name : name
    end

    function __call__(self, x::PyObject, y::PyObject)
        self[:f](Tensor(x), Tensor(y)).o
    end
end


text:
import Base: download
import Dates: unix2datetime, now

@static if Sys.iswindows()
    function _download(url::AbstractString, filename::AbstractString,
        verbose::Bool)
        res = ccall((:URLDownloadToFileW, :urlmon), stdcall, Cuint,
                    (Ptr{Cvoid}, Cwstring, Cwstring, Cuint, Ptr{Cvoid}),
                    C_NULL, url, filename, 0, C_NULL)
        if res != 0
            error("Download failed.")
        end
    end
else
    function _download(url::AbstractString, filename::AbstractString,
        verbose::Bool)
        if verbose
            run(`curl -o $filename -L $url`)
        else
            run(`curl -s -o $filename -L $url`)
        end
    end
end

"""
    download(rf::RemoteFile; quiet::Bool=false, verbose::Bool=false,
        force::Bool=false)

Download `rf`.

- `quiet`: Do not print messages.
- `verbose`: Print all messages.
- `force`: Force download and overwrite existing files.
"""
function download(rf::RemoteFile; verbose::Bool=false, quiet::Bool=false,
    force::Bool=false)
    file = joinpath(rf.dir, rf.file)

    if isfile(file) && !force
        if !isoutdated(rf)
            verbose && @info "File '$(rf.file)' is up-to-date."
            return
        end
    end

    tempfile = tempname()

    if !isdir(rf.dir)
        mkpath(rf.dir)
    end

    !quiet && @info "Downloading file '$(rf.file)' from '$(rf.uri)'."
    tries = 0
    success = false
    while tries < rf.retries
        tries += 1
        try
            _download(string(rf.uri), tempfile, verbose)
            success = true
        catch err
            if isa(err, ErrorException)
                if !quiet
                    @warn "Download failed with error: $(err.msg)"
                    @info "Retrying in $(rf.wait) seconds."
                end
                sleep(rf.wait)
            else
                rethrow(err)
            end
        end
    end

    if success
        update = true
        if isfile(file) && !force
            if samecontent(tempfile, file)
                update = false
                verbose && @info "File '$(rf.file)' has not changed. Update skipped."
            else
                verbose && @info "File '$(rf.file)' was successfully updated."
            end
        else
            verbose && @info "File '$(rf.file)' was successfully downloaded."
        end
        update && mv(tempfile, file, force=true)
    else
        if rf.failed == :warn && isfile(file)
            @warn "Download of '$(rf.uri)' failed after $(rf.retries) retries. Local file was not updated."
        elseif rf.failed == :error || (rf.failed == :warn && !isfile(file))
            error("Download of '$(rf.uri)' failed after $(rf.retries) retries.")
        end
    end
    rm(tempfile, force=true)
end

function samecontent(file1, file2)
    h1 = hash(open(read, file1, "r"))
    h2 = hash(open(read, file2, "r"))
    h1 == h2
end


text:
"""
Individual Types
"""

"""
Any individual is a subtype of AbstractIndividual.
The individual should be a mutable struct with fields that store
data that need to be passed through steps of an experiment
"""
abstract type AbstractIndividual end

"""
A BasicIndividual is a minimal concrete type of the abstract type
AbstractIndividual, and it includes only the fields that are necessary
to complete an experiment.
"""
mutable struct BasicIndividual <: AbstractIndividual
    grn::Array{AbstractFloat,2}
end


text:
import Pkg

function _force_remove_path(path)::Nothing
    try
        rm(path; force = true, recursive = true)
    catch
    end
    return nothing
end

original_directory = pwd()
project_root = joinpath(splitpath(@__DIR__)...)
cd(project_root)

this_registry_toml_path = joinpath(project_root, "Registry.toml")
this_registry_configuration = Pkg.TOML.parsefile(this_registry_toml_path)
if !haskey(this_registry_configuration, "packages")
    this_registry_configuration["packages"] = Dict{String,Any}()
end

this_registry_name_to_path = Dict{String, String}()
this_registry_name_to_uuid = Dict{String, String}()
this_registry_all_packages = String[]
for pair in this_registry_configuration["packages"]
    uuid = pair[1]
    name = pair[2]["name"]
    path = pair[2]["path"]
    push!(this_registry_all_packages, name,)
    this_registry_name_to_path[name] = path
    this_registry_name_to_uuid[name] = uuid
end

unique!(this_registry_all_packages)
sort!(this_registry_all_packages)
n = length(this_registry_all_packages)
@debug(
    "this_registry_all_packages",
    length(this_registry_all_packages),
    this_registry_all_packages,
    repr(this_registry_all_packages),
    )

packages_to_force_delete = String[
    strip(s) for s in ARGS
    ]

unique!(packages_to_force_delete)
sort!(packages_to_force_delete)
@debug(
    "packages_to_force_delete",
    length(packages_to_force_delete),
    packages_to_force_delete,
    repr(packages_to_force_delete),
    )

info_messages = Vector{String}(undef, 0)
n = length(packages_to_force_delete)
for i = 1:n
    name = packages_to_force_delete[i]
    @debug("Deleting \"$(name)\" (package $(i) of $(n))")
    if name in this_registry_all_packages
        our_uuid = this_registry_name_to_uuid[name]
        our_path = this_registry_name_to_path[name]
        our_full_path = joinpath(project_root, our_path)
        rm(our_full_path; force = true, recursive = true)
        delete!(this_registry_configuration["packages"], our_uuid)
    end
    our_hypothetical_path_uppercase = joinpath(
        "packages",
        "$(uppercase(name[1:1]))",
        "$(name)",
        )
    our_hypothetical_path_lowercase = joinpath(
        "packages",
        "$(lowercase(name[1:1]))",
        "$(name)",
        )
    our_full_hypothetical_path_uppercase = joinpath(
        project_root,
        our_hypothetical_path_uppercase,
        )
    our_full_hypothetical_path_lowercase = joinpath(
        project_root,
        our_hypothetical_path_lowercase,
        )
    rm(our_full_hypothetical_path_uppercase; force = true, recursive = true)
    rm(our_full_hypothetical_path_lowercase; force = true, recursive = true)
    push!(info_messages, "Deleted $(name) from our registry")
end

rm(this_registry_toml_path; force = true, recursive = true)
open(this_registry_toml_path, "w") do f
    Pkg.TOML.print(f, this_registry_configuration)
end

temp_registry_toml = read(this_registry_toml_path, String)
temp_registry_toml_parsed = Pkg.TOML.parse(temp_registry_toml)
if !haskey(temp_registry_toml_parsed, "packages")
    temp_output = string(temp_registry_toml,
                         "\n\n\n\n",
                         "[packages]",
                         "\n\n\n\n")
    rm(this_registry_toml_path; force = true, recursive = true)
    write(this_registry_toml_path, temp_output)
end
push!(info_messages, "Wrote Registry.toml file to $(this_registry_toml_path)")

for message in info_messages
    @info(message)
end

cd(original_directory)


text:
module Examples

  using ..AlphaZero

  include("../games/tictactoe/main.jl")
  export Tictactoe

  include("../games/connect-four/main.jl")
  export ConnectFour

  include("../games/grid-world/main.jl")
  export GridWorld

  include("../games/mancala/main.jl")
  export Mancala

  include("../games/jsp/main.jl")
  export JSP

  include("../games/fjsp/main.jl")
  export FJSP

  include("../games/fjspt/main.jl")
  export FJSPT

  const games = Dict(
    "grid-world" => GridWorld.Training.experiment.gspec,
    "tictactoe" => Tictactoe.Training.experiment.gspec,
    "connect-four" => ConnectFour.Training.experiment.gspec,
    "mancala" => Mancala.Training.experiment.gspec,
    "jsp" => JSP.Training.experiment.gspec,
    "fjsp" => FJSP.Training.experiment.gspec,
    "fjspt" => FJSPT.Training.experiment.gspec)
    # "ospiel_ttt" => OSpielTictactoe.GameSpec()
  # ospiel_ttt is added from openspiel_example.jl when OpenSpiel.jl is imported

  const experiments = Dict(
    "grid-world" => GridWorld.Training.experiment,
    "tictactoe" => Tictactoe.Training.experiment,
    "connect-four" => ConnectFour.Training.experiment,
    "mancala" => Mancala.Training.experiment,
    "jsp" => JSP.Training.experiment,
    "fjsp" => FJSP.Training.experiment,
    "fjspt" => FJSPT.Training.experiment)
    # "ospiel_ttt" => OSpielTictactoe.Training.experiment

end


text:
module bearingDistributions
using Distributions
using LinearAlgebra
using Plots

include("App.jl")
using .App

export Bearing, intersectGrid, bisectionRoots, plotConfidenceInterval, App

"""
Bearing object.

Defined by a mean direction θ, a standard deviation of direction σ, and a position (x,y).
"""
struct Bearing
    θ
    σ
    x
    y
end

"""
findIntersect(b1::Bearing, b2::Bearing)

This solves the simultaneous equations of the two paramterized lines as below:
    x₁ + cos(θ₁)t₁ = x₂ + cos(θ₂)t₂
    y₁ + sin(θ₁)t₁ = y₂ + sin(θ₂)t₂

    x₁ - x₂ = -cos(θ₁)t₁ + cos(θ₂)t₂
    y₁ - y₂ = -sin(θ₁)t₁ + sin(θ₂)t₂

    This can be represented in matrix form as:
    b = Ax̄

    where:
    x̄ = [t₁, t₂]
    A = |-cos(θ₁) cos(θ₂)|
        |-sin(θ₁) sin(θ₂)|

    b = [x₁ - x₂, y₁ - y₂]

    Solved by:
    x̄ = A⁻¹ ⋅ b

    This matrix is singular when θ₁ and θ₂ are parallel, i.e.  mod(θ₁,π) == mod(θ₂,π)
    In this case, there is a heuristic that sets the "intersection" as the midpoint between the two points.

    Furthermore, if tᵢ<0 then there is no defined intersection point (lines pointed away from each other)
"""
function findIntersect(b1::Bearing, b2::Bearing)
    if mod(b1.θ-b2.θ,π)≈0
        point = ([b1.x;b1.y] + [b2.x;b2.y])/2
        return point
    end

    A = [-cos(b1.θ) cos(b2.θ);
         -sin(b1.θ) sin(b2.θ)]

    b = [b1.x - b2.x, b1.y - b2.y];

    t = inv(A) * b;
    if any(t .< 0)
        point = ([b1.x;b1.y] + [b2.x;b2.y])/2
        return point
    end
    point = [b1.x, b1.y] + [cos(b1.θ), sin(b1.θ)]*t[1];
    return point
end

"""
boundingBox(b1::Bearing, b2::Bearing)
    Return the box [xₗ xᵤ;yₗ yᵤ] around the intersection of two bearings.

    Size the box such that the origin of the bearings will always be included in the box.
"""
function boundingBox(b1::Bearing, b2::Bearing)
    pᵢ = findIntersect(b1, b2);
    p₁ = [b1.x;b1.y];
    p₂ = [b2.x;b2.y];
    d₁ = pᵢ - p₁;
    d₂ = pᵢ - p₂;
    D₁ = norm(d₁);
    D₂ = norm(d₂);
    D = 1.5*max(D₁, D₂);
    box = pᵢ .+ 1*[-D D;-D D];
    box
end

"""
rangeFromBox(box::Array{Number,2}; length = 51)

    Wrapper to splat out the box array.
"""
function rangeFromBox(box::Array{T,2} where T<:Number ; length = 51)
    x,y = rangeFromBox(box...; length)
    return x,y
end

"""
rangeFromBox(x1,y1,x2,y2; length)

    convert `x1`:`x2` and `y1`:`y2` into ranges of length `length`
"""
function rangeFromBox(x1,y1,x2,y2; length)
    x = range(x1, x2; length);
    y = range(y1, y2; length);
    return x,y
end

"""
probabilityGrid(b::Bearing, x::AbstractArray, y::AbstractArray)

    Compute the probability grid over ranges `x` and `y` of bearing `b`
"""
function probabilityGrid(b::Bearing, x::AbstractArray, y::AbstractArray)
    θ = mod(b.θ + π,2π) - π
    D = Normal(0, b.σ);
    Θ = atan.(y .- b.y, x' .- b.x) .- θ;
    Θ = mod.(Θ .+ π, 2π) .-π
    P = pdf.(D, Θ);
    ΔA = convert(Float64,x.step * y.step);
    P = P./sum(P[:])/ΔA;
    return x, y, P
end

"""
intersectGrid(b1::Bearing, b2::Bearing; length = 51)

    Compute the probability grid of the intersection of `b1` and `b2` over ranges `x` and `y`

    Return the probability grids for `b1`, `b2`, and the intersection.

"""
function intersectGrid(b1::Bearing, b2::Bearing; length = 51)
    box = boundingBox(b1, b2);
    x,y = rangeFromBox(box...; length)
    _,_,P₁ = probabilityGrid(b1, x, y);
    _,_,P₂ = probabilityGrid(b2, x, y);
    Pᵢ = P₁ .* P₂
    ΔA = convert(Float64,x.step * y.step);
    Pᵢ = Pᵢ ./ sum(Pᵢ[:])
    Pᵢ = Pᵢ ./ ΔA;
    return x, y, Pᵢ, P₁, P₂
end

"""
plotProbabilityGrids(x::AbstractArray, y::AbstractArray, P...)

    Plot of tuple of probability grids `P` on ranges x and y with the default contour levels
"""
function plotProbabilityGrids(x::AbstractArray, y::AbstractArray, P...)
    p = contour(;aspect_ratio=:equal);
    for P in P
        contour!(x,y,P)
    end
    display(p)
end

"""
plotConfidenceInterval(α::Number, x::AbstractArray, y::AbstractArray, P...)

    Plot the tuple of probability grids `P` on ranges x and y at the defined confidence interval `α`
"""
function plotConfidenceInterval(α::Number, x::AbstractArray, y::AbstractArray, P...)
    handle = contour(;aspect_ratio=:equal);
    ΔA = convert(Float64,x.step*y.step);
    zeros = [bisectionRoots((p) -> ΔA*sum(Pⱼ[Pⱼ .> p]) - α, minimum(Pⱼ),maximum(Pⱼ)) for Pⱼ in P]
    for (i,P) in enumerate(P)
        contour!(x,y,P; levels=[zeros[i]])
    end
    return handle
end

"""
bisectionRoots(f::Function,xₗ::Number,xᵤ::Number; max_iteration = 100, xtol = 1e-3)

    Computes the root of function `f(x)` on the bounded range `[xₗ xᵤ]` with the Bisection algorithm.

    Exit condition is defined by xᵤ-xₗ < xtol
"""
function bisectionRoots(f::Function,xₗ::Number,xᵤ::Number; max_iteration = 100, xtol = 1e-3)
    fₗ  = f(xₗ);
    fᵤ = f(xᵤ);

    if fₗ == 0
        return xₗ
    elseif fᵤ == 0
        return xᵤ
    end

    if (fₗ<0 && fᵤ<0) || (fₗ>0 && fᵤ>0)
        return NaN
    end

    xᵢ = (xₗ + xᵤ)/2;
    iter = 1;

    while abs(xᵤ-xₗ) > xtol
        iter > max_iteration ? break : iter+=1;
        if sign(f(xᵢ))==sign(fₗ)
            xₗ = xᵢ
            xᵢ = (xᵢ + xᵤ)/2
        else
            xᵤ = xᵢ
            xᵢ = (xₗ + xᵢ)/2
        end
    end
    return xᵢ
end

end


text:
function Vector1D(arg0::jdouble)
    return Vector1D((jdouble,), arg0)
end

function Vector1D(arg0::jdouble, arg1::Vector1D)
    return Vector1D((jdouble, Vector1D), arg0, arg1)
end

function Vector1D(arg0::jdouble, arg1::Vector1D, arg2::jdouble, arg3::Vector1D)
    return Vector1D((jdouble, Vector1D, jdouble, Vector1D), arg0, arg1, arg2, arg3)
end

function Vector1D(arg0::jdouble, arg1::Vector1D, arg2::jdouble, arg3::Vector1D, arg4::jdouble, arg5::Vector1D)
    return Vector1D((jdouble, Vector1D, jdouble, Vector1D, jdouble, Vector1D), arg0, arg1, arg2, arg3, arg4, arg5)
end

function Vector1D(arg0::jdouble, arg1::Vector1D, arg2::jdouble, arg3::Vector1D, arg4::jdouble, arg5::Vector1D, arg6::jdouble, arg7::Vector1D)
    return Vector1D((jdouble, Vector1D, jdouble, Vector1D, jdouble, Vector1D, jdouble, Vector1D), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function add(obj::Vector1D, arg0::JVector)
    return jcall(obj, "add", Vector1D, (JVector,), arg0)
end

function add(obj::Vector1D, arg0::jdouble, arg1::JVector)
    return jcall(obj, "add", Vector1D, (jdouble, JVector), arg0, arg1)
end

function distance(::Type{Vector1D}, arg0::Vector1D, arg1::Vector1D)
    return jcall(Vector1D, "distance", jdouble, (Vector1D, Vector1D), arg0, arg1)
end

function distance(obj::Vector1D, arg0::Point)
    return jcall(obj, "distance", jdouble, (Point,), arg0)
end

function distance1(obj::Vector1D, arg0::JVector)
    return jcall(obj, "distance1", jdouble, (JVector,), arg0)
end

function distance_inf(::Type{Vector1D}, arg0::Vector1D, arg1::Vector1D)
    return jcall(Vector1D, "distanceInf", jdouble, (Vector1D, Vector1D), arg0, arg1)
end

function distance_inf(obj::Vector1D, arg0::JVector)
    return jcall(obj, "distanceInf", jdouble, (JVector,), arg0)
end

function distance_sq(::Type{Vector1D}, arg0::Vector1D, arg1::Vector1D)
    return jcall(Vector1D, "distanceSq", jdouble, (Vector1D, Vector1D), arg0, arg1)
end

function distance_sq(obj::Vector1D, arg0::JVector)
    return jcall(obj, "distanceSq", jdouble, (JVector,), arg0)
end

function dot_product(obj::Vector1D, arg0::JVector)
    return jcall(obj, "dotProduct", jdouble, (JVector,), arg0)
end

function equals(obj::Vector1D, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_norm(obj::Vector1D)
    return jcall(obj, "getNorm", jdouble, ())
end

function get_norm1(obj::Vector1D)
    return jcall(obj, "getNorm1", jdouble, ())
end

function get_norm_inf(obj::Vector1D)
    return jcall(obj, "getNormInf", jdouble, ())
end

function get_norm_sq(obj::Vector1D)
    return jcall(obj, "getNormSq", jdouble, ())
end

function get_space(obj::Vector1D)
    return jcall(obj, "getSpace", Space, ())
end

function get_x(obj::Vector1D)
    return jcall(obj, "getX", jdouble, ())
end

function get_zero(obj::Vector1D)
    return jcall(obj, "getZero", Vector1D, ())
end

function hash_code(obj::Vector1D)
    return jcall(obj, "hashCode", jint, ())
end

function is_infinite(obj::Vector1D)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::Vector1D)
    return jcall(obj, "isNaN", jboolean, ())
end

function negate(obj::Vector1D)
    return jcall(obj, "negate", Vector1D, ())
end

function normalize(obj::Vector1D)
    return jcall(obj, "normalize", Vector1D, ())
end

function scalar_multiply(obj::Vector1D, arg0::jdouble)
    return jcall(obj, "scalarMultiply", Vector1D, (jdouble,), arg0)
end

function subtract(obj::Vector1D, arg0::JVector)
    return jcall(obj, "subtract", Vector1D, (JVector,), arg0)
end

function subtract(obj::Vector1D, arg0::jdouble, arg1::JVector)
    return jcall(obj, "subtract", Vector1D, (jdouble, JVector), arg0, arg1)
end

function to_string(obj::Vector1D)
    return jcall(obj, "toString", JString, ())
end

function to_string(obj::Vector1D, arg0::NumberFormat)
    return jcall(obj, "toString", JString, (NumberFormat,), arg0)
end


text:
using CUDArt
@everywhere CUDArt.device((myid()-1) % CUDArt.devcount())
using KUparser, KUnet, HDF5, JLD, ArgParse, Compat
VERSION < v"0.4-" && eval(Expr(:using,:Dates))

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "datafiles"
        help = "JLD files for input: first one will be used for training"
        nargs = '+'
        required = true
        "--out"
        help = "JLD file for saving last network"
        #default = ""
        "--in"
        help = "JLD file for loading init network"
        #default = ""
        "--best"
        help = "JLD file for best (dev) network"
        #default = ""
        "--nogpu"
        help = "Do not use gpu"
        action = :store_true
        "--ispunct"
        help = "Do not score punct as determined by this function"
        #default = "zn11punct"
        "--parser"
        help = "Parsing algorithm to use: g(reedy)parser, b(eam)parser, or o(racle)parser for static training"
        default = "bparser"
        "--arctype"
        help = "Move set to use: ArcEager{R1,13}, ArcHybrid{R1,13}"
        default = "ArcEager13"
        "--feats"
        help = "Features to use from KUparser.Flist"
        default = "zn11cpv"
        "--hidden"
        help = "One or more hidden layer sizes"
        nargs = '+'
        arg_type = Int
        #default = [16384]
        "--minepochs"
        help = "Minimum number of epochs to train"
        arg_type = Int
        #default = 100
        "--maxepochs"
        help = "Maximum number of epochs to train"
        arg_type = Int
        #default = 100
        "--epochs"
        help = "Exact number of epochs to train"
        arg_type = Int
        #default = 100
        # Just start Julia with desired workers instead
        # "--ncpu"
        # help = "Number of workers for multithreaded parsing"
        # arg_type = Int
        #default = 16
        "--usepmap"
        help = "Use multiple processes"
        action = :store_true
        "--strain"
        help = "Number of sentences to parse for a training batch"
        arg_type = Int
        default = 2500
        "--stest"
        help = "Number of sentences to parse for a testing batch"
        arg_type = Int
        default = 2500
        "--pbatch"
        help = "Minibatch argument for parse"
        arg_type = Int
        default = 25
        "--tbatch"
        help = "Minibatch argument for train"
        arg_type = Int
        default = 128
        "--nbeam"
        help = "Beam size for bparser"
        arg_type = Int
        default = 64
        "--adagrad"
        help = "Use adagrad"
        action = :store_true
        # arg_type = Float32
        # nargs = '+'
        # default = [1f-8]
        "--dropout"
        help = "Dropout probability"
        arg_type = Float32
        nargs = '+'
        # default = [0.1f0, 0.5f0]
        "--l1reg"
        help = "L1 regularization parameter"
        arg_type = Float32
        nargs = '+'
        "--l2reg"
        help = "L2 regularization parameter"
        arg_type = Float32
        nargs = '+'
        "--lr"
        help = "Learning rate"
        arg_type = Float32
        nargs = '+'
        default = [0.01f0]
        "--maxnorm"
        help = "If nonzero upper limit on weight matrix row norms"
        arg_type = Float32
        nargs = '+'
        "--momentum"
        help = "Momentum"
        arg_type = Float32
        nargs = '+'
        "--nesterov"
        help = "Apply nesterov's accelerated gradient"
        arg_type = Float32
        nargs = '+'
        "--seed"
        help = "Random number seed (-1 leaves the default seed)"
        arg_type = Int
        default = -1
        "--shuffle"
        help = "Whether to shuffle the data every epoch"
        action = :store_true
    end
    args = parse_args(s)
    for (arg,val) in args
        print("$arg:$val ")
    end
    println("")
    return args
end

function getscore(parses, corpus, ispunct)
    if ispunct == nothing
        @show e = evalparse(parses, corpus)
        return e[1]
    else
        @show e = evalparse(parses, corpus; ispunct=ispunct)
        return e[2]
    end
end

function loaddata(files)
    data = Corpus[]
    deprel = postag = nothing
    ndeps = 0
    for f in files
        @show f
        @date d=load(f)
        push!(data, d["corpus"])
        if deprel == nothing
            deprel = d["deprel"]
            postag = d["postag"]
            ndeps = length(deprel)
        else
            @assert deprel == d["deprel"]
            @assert postag == d["postag"]
        end
    end
    return (data, ndeps)
end

function initnet(args, pt, ndeps)
    args["in"] != nothing && !isfile(args["in"]) && (warn("Cannot find net file, creating blank net"); args["in"]=nothing)
    if args["in"] == nothing
        net = Layer[]
        length(args["dropout"]) > 1 && push!(net, Drop(args["dropout"][1]))
        for h in args["hidden"]
            append!(net, [Mmul(h), Bias(), Relu()])
            !isempty(args["dropout"]) && push!(net, Drop(args["dropout"][end]))
        end
        yrows = pt(1,ndeps).nmove
        append!(net, [Mmul(yrows), Bias()])
        push!(net, args["parser"]=="bparser" ? ScalLoss() : XentLoss())
    else
        @date net = loadnet(args["in"])
    end
    for k in fieldnames(KUparam)
        haskey(args, string(k)) || continue
        v = args[string(k)]
        if isempty(v)
            continue
        elseif length(v)==1
            @eval setparam!($net; $k=$v[1])
        else
            @assert length(v)==length(net) "$k should have 1 or $(length(net)) elements"
            for i=1:length(v)
                @eval setparam!($net[$i]; $k=$v[$i])
            end
        end
    end
    return net
end

function main()
    args = parse_commandline()
    @show nworkers()
    KUnet.gpu(!args["nogpu"])
    args["nogpu"] && blas_set_num_threads(20)
    args["seed"] >= 0 && setseed(args["seed"])
    (data, ndeps) = loaddata(args["datafiles"])
    @show (map(size, data), ndeps)
    @show feats = eval(parse("Flist."*args["feats"]))
    @show pt = eval(parse(args["arctype"]))
    net = initnet(args, pt, ndeps)
    display(net); println()
    ispunct = (args["ispunct"]==nothing ? nothing : eval(parse("KUparser."*args["ispunct"])))
    @show ispunct
    accuracy = Array(Float32, length(data))
    (bestscore,bestepoch,epoch)=(0,0,0)

    while true
        println("\n$(now()) epoch => $(epoch += 1)")
        icorpus=1
        corpus=data[icorpus]          # use data[1] as the training corpus
        parses=pt[]
        for c1=1:args["strain"]:length(corpus)
            c2=min(length(corpus), c1+args["strain"]-1)
            sentences=sub(corpus, c1:c2)
            @show (epoch, icorpus, c1, c2)
            if (args["parser"] == "oparser") # || ((epoch==1) && (args["in"]==nothing)))
                @date (p,x,y) = oparse(pt, sentences, ndeps, feats; usepmap=args["usepmap"])
            elseif (args["parser"] == "gparser")
                @date (p,x,y) = gparse(pt, sentences, ndeps, feats, net; usepmap=args["usepmap"], nbatch=args["pbatch"], returnxy=true)
            elseif (args["parser"] == "bparser")
                @date (p,x,y) = bparse(pt, sentences, ndeps, feats, net, args["nbeam"]; usepmap=args["usepmap"], nbatch=args["pbatch"], returnxy=true)
            else
                error("Unknown parser")
            end
            @show map(size, (p,x,y))
            @assert isa(p, Vector{pt})
            @assert isa(x, Vector{Matrix{wtype(sentences)}})
            @assert isa(y, Vector{Matrix{wtype(sentences)}})
            @assert length(x)==length(y)
            @everywhere gc(); @everywhere gc()
            append!(parses, p)
            @date for i=1:length(x)
                train(net, x[i], y[i]; batch=args["tbatch"], shuffle=args["shuffle"])
            end
        end
        accuracy[icorpus] = getscore(parses, corpus, ispunct)
        (args["out"] != nothing) && (@date savenet(args["out"], net))

        for icorpus=2:length(data)
            corpus=data[icorpus]
            parses=pt[]
            for c1=1:args["stest"]:length(corpus)
                c2=min(length(corpus), c1+args["stest"]-1)
                @show (:corpus, icorpus, c1, c2)
                sentences=sub(corpus, c1:c2)
                if in(args["parser"], ("oparser", "gparser"))
                    @date p = gparse(pt, sentences, ndeps, feats, net; usepmap=args["usepmap"], nbatch=args["pbatch"], returnxy=false)
                elseif (args["parser"] == "bparser")
                    @date p = bparse(pt, sentences, ndeps, feats, net, args["nbeam"]; usepmap=args["usepmap"], nbatch=args["pbatch"], returnxy=false)
                else
                    error("Unknown parser")
                end
                @show size(p)
                @assert isa(p, Vector{pt})
                @everywhere gc(); @everywhere gc()
                append!(parses, p)
            end
            accuracy[icorpus] = getscore(parses, corpus, ispunct)
        end

        println("DATA:\t$epoch\t"*join(accuracy, '\t')); flush(STDOUT)

        # We look at accuracy[2] (dev score) by default,
        # accuracy[1] (trn score) if there is no dev.

        myscore = (length(accuracy) > 1 ? accuracy[2] : accuracy[1])

        # Update best score and save best net

        if (myscore > bestscore)
            @show (bestscore,bestepoch)=(myscore,epoch)
            (args["best"] != nothing) && (@date savenet(args["best"], net))
        end

        # If epochs specified, do exactly epoch.
        # Otherwise quit if we have not made progress in the last half.
        # If minepochs specified, do not quit before minepoch.
        # If maxepochs specified, quit before maxepoch.

        if (args["epochs"] != nothing)
            (epoch >= args["epochs"]) && break
        elseif (epoch >= 2*bestepoch)
            (args["minepochs"]==nothing || epoch >= args["minepochs"]) && break
        elseif (args["maxepochs"]!=nothing)
            (epoch >= args["maxepochs"]) && break
        end

    end # while true
end # main

main()

# DEAD CODE
    # Initialize training set using oparse on first corpus
    # @date (p,x,y) = oparse(pt, data[1], ndeps, args["ncpu"], feats)
    # @date (p,x,y) = oparse(pt, data[1], ndeps, feats)
    # @show evalparse(p, data[1]); p=nothing
    # for i=2:length(data)
    #     #@date p = oparse(pt, data[i], ndeps, args["ncpu"])
    #     @date p = oparse(pt, data[i], ndeps)
    #     @show evalparse(p, data[i]); p=nothing
    # end

    # else # if DBG
    #     @date net = loadnet("foo11nnet3.jld")
    #     # @date d = load("zn11cpv.jld")
    #     # (x,y) = (d["xtrn"],d["ytrn"])
    #     @date x = Array(Float32, 3105,1820392)
    #     @date y = Array(Float32, 24,1820392)
    # end  # if DBG
        # @meminfo
        # if !DBG
        #     if isa(eltype(x), Array)
        #         @assert length(x)==length(y)
        #         for i=1:length(x)
        #             @show map(size, (x[i], y[i]))
        #             @date train(net, x[i], y[i]; batch=args["tbatch"], shuffle=true)
        #         end
        #     else
        #         @show map(size, (x, y))
        #         @date train(net, x, y; batch=args["tbatch"], shuffle=true)
        #     end
        # end # if !DBG

        # @meminfo
        # if (args["parser"] == "bparser")
        #     # The first corpus gives us the new training set
        #     p = x = y = nothing; gc()
        #     @date (@everywhere gc())
        #     @date (@everywhere require("KUparser"))
        #     @date (for w in workers(); @fetchfrom w require("errlog.jl"); end)
        #     @meminfo
        #     @date pxy = bparse(pt, data[1], ndeps, feats, net, args["nbeam"], args["pbatch"], args["ncpu"]; xy=true)
        #     @meminfo
        #     p = map(z->z[1], pxy); @show map(size, p)
        #     x = map(z->z[2], pxy); @show map(size, x)
        #     y = map(z->z[3], pxy); @show map(size, y)
        #     @meminfo
        #     #@show map(getbytes, (p,x,y))
        #     @show e = evalparse(mycat(p), data[1]); p=nothing
        #     accuracy[1] = e[2]
        #     for i=2:length(data)
        #         @date p = bparse(pt, data[i], ndeps, feats, net, args["nbeam"], args["pbatch"], args["ncpu"])
        #         @show e = evalparse(mycat(p), data[i]); p=nothing
        #         accuracy[i] = e[2]
        #     end
        #     #@date Main.restartmachines()
        #     #@date Main.restartcuda()
        # elseif (args["parser"] == "oparser")
        #     # We never change the training set with oparser, just report accuracy
        #     for i=1:length(data)
        #         @date p = gparse(pt, data[i], ndeps, feats, net, args["pbatch"], args["ncpu"])
        #         @show e = evalparse(p, data[i]); p=nothing
        #         accuracy[i] = e[2]  # e[2] is UAS excluding punct
        #     end
        # elseif (args["parser"] == "gparser")
        #     # The first corpus gives us the new training set
        #     p = x = y = nothing;
        #     @date (@everywhere gc())
        #     @date (@everywhere require("KUparser"))
        #     @date (@everywhere require("errlog.jl"))
        #     @date (p,x,y) = gparse(pt, data[1], ndeps, feats, net, args["pbatch"], args["ncpu"]; xy=true)
        #     @show e = evalparse(p, data[1]); p=nothing
        #     accuracy[1] = e[2]
        #     for i=2:length(data)
        #         @date p = gparse(pt, data[i], ndeps, feats, net, args["pbatch"], args["ncpu"])
        #         @show e = evalparse(p, data[i]); p=nothing
        #         accuracy[i] = e[2]
        #     end
        #     @date Main.restartmachines()
        #     @date Main.restartcuda()
        # else
        #     error("Unknown parser "*args["parser"])
        # end
#macro meminfo() :(gc(); run(`nvidia-smi`); run(`ps auxww`|>`grep julia`); run(`free`)) end
#macro meminfo() :(run(`ps auxww`|>`grep julia`); run(`free`)) end
#macro meminfo() :(nothing) end

# function mycat(a::Array)
#     # concat elements of a along their last dimension
#     mx = size(a[1])[1:end-1]
#     nx = 0
#     for x in a
#         @assert size(x)[1:end-1] == mx
#         nx += size(x)[end]
#     end
#     b = similar(a[1], tuple(mx..., nx))
#     nb = 0
#     for x in a
#         copy!(b, nb+1, x, 1, length(x))
#         nb += length(x)
#     end
#     return b
# end

        # # gc() does not work, so we need these:
        # @date sleep(5)
        # @date Main.restartmachines()
        # @date sleep(5)
        # @date Main.restartcuda()
        # @date (@everywhere require("KUparser"))
        # @date (@everywhere gc())


text:
using ADMPS
using ADMPS: num_grad,Zofβ,logoverlap,Z,obs_env,magofβ,eneofβ,overlap,onestep,isingβc,init_mps
using CUDA
using KrylovKit
using LinearAlgebra: svd, norm
using LineSearches, Optim
using OMEinsum
using Plots
using Random
using Test
using Zygote

@testset "gradient with $atype{$dtype}" for atype in [Array], dtype in [ComplexF64]
    Random.seed!(100)
    D,χ = 2,5
    β = 0.4
    model = Ising(β)
    M = atype{dtype}(model_tensor(model))
    mps = init_mps(D = D, χ = χ)
    Au = mps
    Ad = mps
    ff(Ad) = logoverlap(Au, Ad, M)
    # @show logoverlap(Au, mps, M),Zygote.gradient(ff,Ad)
    gradzygote = first(Zygote.gradient(mps) do x
        logoverlap(Au, x, M)
    end)
    gradnum = num_grad(mps, δ=1e-4) do x
        logoverlap(Au, x, M)
    end
    @test gradzygote ≈ gradnum atol=1e-5
end

@testset "onestep optimize Ad with $atype{$dtype}" for atype in [Array], dtype in [ComplexF64]
    seed_number = 100
    β = 0.8
    χ = 10
    maxiter = 5
    infolder, outfolder = "./data/", "./data/"
    model = Ising(β)
    M = atype{dtype}(model_tensor(model))
    Random.seed!(seed_number)

    Ad = onestep(M; infolder = infolder*"$(model)/", outfolder = outfolder*"$(model)/", χ = χ, verbose= true, savefile = false)
    @test Ad !== nothing
end

@testset "oneside optimize mps" for atype in [Array], dtype in [ComplexF64]
    seed_number = 100
    β = 0.8
    D,χ = 2,20
    mapsteps = 20
    infolder, outfolder = "./data/", "./data/"

    model = Ising(β)
    M = atype{dtype}(model_tensor(model))
    Random.seed!(seed_number)

    Au, Ad = optimisemps(M; infolder = infolder*"$(model)/", outfolder = outfolder*"$(model)/", χ = χ, mapsteps = mapsteps, verbose= true, updown = false)

    env = obs_env(M,Au,Ad)
    @test magnetisation(env,model) ≈ magofβ(model) atol=1e-6
    @test energy(env,model) ≈ eneofβ(model) atol=1e-6
    # @show energy(env,model)+1.414213779415974 # β = isingβc
end

@testset "twoside optimize mps" for atype in [Array], dtype in [ComplexF64]
    seed_number = 100
    β = 0.8
    D,χ = 2,20
    mapsteps = 20
    infolder, outfolder = "./data/", "./data/"

    model = Ising(β)
    M = atype{dtype}(model_tensor(model))
    Random.seed!(seed_number)

    Au, Ad = optimisemps(M; infolder = infolder*"$(model)/", outfolder = outfolder*"$(model)/", χ = χ, mapsteps = mapsteps, verbose= true, downfromup = true)

    env = obs_env(M,Au,Ad)
    @test magnetisation(env,model) ≈ magofβ(model) atol=1e-6
    @test energy(env,model) ≈ eneofβ(model) atol=1e-6
end


text:
# Testing:
#
#  - computation of sufficient statistics
#  - distribution fitting (i.e. estimation)
#

using Distributions
using Test, Random, LinearAlgebra

n0 = 100
N = 10^5

w = rand(n0)

# DiscreteUniform

x = rand(DiscreteUniform(10, 15), n0)
d = fit(DiscreteUniform, x)
@test isa(d, DiscreteUniform)
@test minimum(d) == minimum(x)
@test maximum(d) == maximum(x)

d = fit(DiscreteUniform, rand(DiscreteUniform(10, 15), N))
@test minimum(d) == 10
@test maximum(d) == 15

# Bernoulli

x = rand(Bernoulli(0.7), n0)

ss = suffstats(Bernoulli, x)
@test isa(ss, Distributions.BernoulliStats)
@test ss.cnt0 == n0 - count(t->t != 0, x)
@test ss.cnt1 == count(t->t != 0, x)

ss = suffstats(Bernoulli, x, w)
@test isa(ss, Distributions.BernoulliStats)
@test ss.cnt0 ≈ sum(w[x .== 0])
@test ss.cnt1 ≈ sum(w[x .== 1])

d = fit(Bernoulli, x)
p = count(t->t != 0, x) / n0
@test isa(d, Bernoulli)
@test mean(d) ≈ p

d = fit(Bernoulli, x, w)
p = sum(w[x .== 1]) / sum(w)
@test isa(d, Bernoulli)
@test mean(d) ≈ p

d = fit(Bernoulli, rand(Bernoulli(0.7), N))
@test isa(d, Bernoulli)
@test isapprox(mean(d), 0.7, atol=0.01)

# Beta

d = fit(Beta, rand(Beta(1.3, 3.7), N))
@test isa(d, Beta)
@test isapprox(d.α, 1.3, atol=0.1)
@test isapprox(d.β, 3.7, atol=0.1)

# Binomial

x = rand(Binomial(100, 0.3), n0)

ss = suffstats(Binomial, (100, x))
@test isa(ss, Distributions.BinomialStats)
@test ss.ns ≈ sum(x)
@test ss.ne == n0
@test ss.n == 100

ss = suffstats(Binomial, (100, x), w)
@test isa(ss, Distributions.BinomialStats)
@test ss.ns ≈ dot(Float64[xx for xx in x], w)
@test ss.ne ≈ sum(w)
@test ss.n == 100

d = fit(Binomial, (100, x))
@test isa(d, Binomial)
@test ntrials(d) == 100
@test succprob(d) ≈ sum(x) / (n0 * 100)

d = fit(Binomial, (100, x), w)
@test isa(d, Binomial)
@test ntrials(d) == 100
@test succprob(d) ≈ dot(x, w) / (sum(w) * 100)

d = fit(Binomial, 100, rand(Binomial(100, 0.3), N))
@test isa(d, Binomial)
@test ntrials(d) == 100
@test isapprox(succprob(d), 0.3, atol=0.01)

# Categorical

p = [0.2, 0.5, 0.3]
x = rand(Categorical(p), n0)

ss = suffstats(Categorical, (3, x))
h = Float64[count(v->v == i, x) for i = 1 : 3]
@test isa(ss, Distributions.CategoricalStats)
@test ss.h ≈ h

d = fit(Categorical, (3, x))
@test isa(d, Categorical)
@test d.K == 3
@test probs(d) ≈ h / sum(h)

d2 = fit(Categorical, x)
@test isa(d2, Categorical)
@test probs(d2) == probs(d)

ss = suffstats(Categorical, (3, x), w)
h = Float64[sum(w[x .== i]) for i = 1 : 3]
@test isa(ss, Distributions.CategoricalStats)
@test ss.h ≈ h

d = fit(Categorical, (3, x), w)
@test isa(d, Categorical)
@test probs(d) ≈ h / sum(h)

d = fit(Categorical, suffstats(Categorical, 3, x, w))
@test isa(d, Categorical)
@test probs(d) ≈ (h / sum(h))

d = fit(Categorical, rand(Categorical(p), N))
@test isa(d, Categorical)
@test isapprox(probs(d), p, atol=0.01)

# Cauchy

@test fit(Cauchy, collect(-4.0:4.0)) === Cauchy(0.0, 2.0)

# Exponential

x = rand(Exponential(0.5), n0)

ss = suffstats(Exponential, x)
@test isa(ss, Distributions.ExponentialStats)
@test ss.sx ≈ sum(x)
@test ss.sw == n0

ss = suffstats(Exponential, x, w)
@test isa(ss, Distributions.ExponentialStats)
@test ss.sx ≈ dot(x, w)
@test ss.sw == sum(w)

d = fit(Exponential, x)
@test isa(d, Exponential)
@test scale(d) ≈ mean(x)

d = fit(Exponential, x, w)
@test isa(d, Exponential)
@test scale(d) ≈ dot(x, w) / sum(w)

d = fit(Exponential, rand(Exponential(0.5), N))
@test isa(d, Exponential)
@test isapprox(scale(d), 0.5, atol=0.01)

# Normal

μ = 11.3
σ = 3.2

x = rand(Normal(μ, σ), n0)

ss = suffstats(Normal, x)
@test isa(ss, Distributions.NormalStats)
@test ss.s  ≈ sum(x)
@test ss.m  ≈ mean(x)
@test ss.s2 ≈ sum((x .- ss.m).^2)
@test ss.tw ≈ n0

ss = suffstats(Normal, x, w)
@test isa(ss, Distributions.NormalStats)
@test ss.s  ≈ dot(x, w)
@test ss.m  ≈ dot(x, w) / sum(w)
@test ss.s2 ≈ dot((x .- ss.m).^2, w)
@test ss.tw ≈ sum(w)

d = fit(Normal, x)
@test isa(d, Normal)
@test d.μ ≈ mean(x)
@test d.σ ≈ sqrt(mean((x .- d.μ).^2))

d = fit(Normal, x, w)
@test isa(d, Normal)
@test d.μ ≈ dot(x, w) / sum(w)
@test d.σ ≈ sqrt(dot((x .- d.μ).^2, w) / sum(w))

d = fit(Normal, rand(Normal(μ, σ), N))
@test isa(d, Normal)
@test isapprox(d.μ, μ, atol=0.1)
@test isapprox(d.σ, σ, atol=0.1)

import Distributions.NormalKnownMu, Distributions.NormalKnownSigma

ss = suffstats(NormalKnownMu(μ), x)
@test isa(ss, Distributions.NormalKnownMuStats)
@test ss.μ == μ
@test ss.s2 ≈ sum(abs2.(x .- μ))
@test ss.tw ≈ n0

ss = suffstats(NormalKnownMu(μ), x, w)
@test isa(ss, Distributions.NormalKnownMuStats)
@test ss.μ == μ
@test ss.s2 ≈ dot((x .- μ).^2, w)
@test ss.tw ≈ sum(w)

d = fit_mle(Normal, x; mu=μ)
@test isa(d, Normal)
@test d.μ == μ
@test d.σ ≈ sqrt(mean((x .- d.μ).^2))

d = fit_mle(Normal, x, w; mu=μ)
@test isa(d, Normal)
@test d.μ == μ
@test d.σ ≈ sqrt(dot((x .- d.μ).^2, w) / sum(w))

ss = suffstats(NormalKnownSigma(σ), x)
@test isa(ss, Distributions.NormalKnownSigmaStats)
@test ss.σ == σ
@test ss.sx ≈ sum(x)
@test ss.tw ≈ n0

ss = suffstats(NormalKnownSigma(σ), x, w)
@test isa(ss, Distributions.NormalKnownSigmaStats)
@test ss.σ == σ
@test ss.sx ≈ dot(x, w)
@test ss.tw ≈ sum(w)

d = fit_mle(Normal, x; sigma=σ)
@test isa(d, Normal)
@test d.σ == σ
@test d.μ ≈ mean(x)

d = fit_mle(Normal, x, w; sigma=σ)
@test isa(d, Normal)
@test d.σ == σ
@test d.μ ≈ dot(x, w) / sum(w)

# Uniform

x = rand(Uniform(1.2, 5.8), n0)
d = fit(Uniform, x)
@test isa(d, Uniform)
@test 1.2 <= minimum(d) <= maximum(d) <= 5.8
@test minimum(d) == minimum(x)
@test maximum(d) == maximum(x)

d = fit(Uniform, rand(Uniform(1.2, 5.8), N))
@test 1.2 <= minimum(d) <= maximum(d) <= 5.8
@test isapprox(minimum(d), 1.2, atol=0.02)
@test isapprox(maximum(d), 5.8, atol=0.02)

# Gamma

x = rand(Gamma(3.9, 2.1), n0)

ss = suffstats(Gamma, x)
@test isa(ss, Distributions.GammaStats)
@test ss.sx    ≈ sum(x)
@test ss.slogx ≈ sum(log.(x))
@test ss.tw    ≈ n0

ss = suffstats(Gamma, x, w)
@test isa(ss, Distributions.GammaStats)
@test ss.sx    ≈ dot(x, w)
@test ss.slogx ≈ dot(log.(x), w)
@test ss.tw    ≈ sum(w)

d = fit(Gamma, rand(Gamma(3.9, 2.1), N))
@test isa(d, Gamma)
@test isapprox(shape(d), 3.9, atol=0.1)
@test isapprox(scale(d), 2.1, atol=0.2)

# Geometric

x = rand(Geometric(0.3), n0)

ss = suffstats(Geometric, x)
@test isa(ss, Distributions.GeometricStats)
@test ss.sx ≈ sum(x)
@test ss.tw ≈ n0

ss = suffstats(Geometric, x, w)
@test isa(ss, Distributions.GeometricStats)
@test ss.sx ≈ dot(x, w)
@test ss.tw ≈ sum(w)

d = fit(Geometric, x)
@test isa(d, Geometric)
@test succprob(d) ≈ inv(1. + mean(x))

d = fit(Geometric, x, w)
@test isa(d, Geometric)
@test succprob(d) ≈ inv(1. + dot(x, w) / sum(w))

d = fit(Geometric, rand(Geometric(0.3), N))
@test isa(d, Geometric)
@test isapprox(succprob(d), 0.3, atol=0.01)

# Laplace

d = fit(Laplace, rand(Laplace(5.0, 3.0), N))
@test isa(d, Laplace)
@test isapprox(location(d), 5.0, atol=0.1)
@test isapprox(scale(d)   , 3.0, atol=0.2)

# Pareto

x = rand(Pareto(3., 7.), N)
d = fit(Pareto, x)

@test isa(d, Pareto)
@test isapprox(shape(d), 3., atol=0.1)
@test isapprox(scale(d), 7., atol=0.1)

# Poisson

x = rand(Poisson(8.2), n0)

ss = suffstats(Poisson, x)
@test isa(ss, Distributions.PoissonStats)
@test ss.sx ≈ sum(x)
@test ss.tw ≈ n0

ss = suffstats(Poisson, x, w)
@test isa(ss, Distributions.PoissonStats)
@test ss.sx ≈ dot(x, w)
@test ss.tw ≈ sum(w)

d = fit(Poisson, x)
@test isa(d, Poisson)
@test mean(d) ≈ mean(x)

d = fit(Poisson, x, w)
@test isa(d, Poisson)
@test mean(d) ≈ dot(Float64[xx for xx in x], w) / sum(w)

d = fit(Poisson, rand(Poisson(8.2), N))
@test isa(d, Poisson)
@test isapprox(mean(d), 8.2, atol=0.2)


text:
𝑃, 𝑂 = spacecraft_rendezvous(T=200.0, orderT=10, orderQ=2, abs_tol=1e-28, max_steps=5000);

𝑂jets = Options(:orderT=>orderT, :orderQ=>orderQ, :abs_tol=>abs_tol, :max_steps=>max_steps)

sol = solve(𝑃, 𝑂, TMJets(𝑂jets), LazyDiscretePost(:check_invariant_intersection=>true))

# first mode
#p = IVP(𝑃.s.modes[1], 𝑃.x0[1][2])
#solve(p, 𝑂, op=TMJets(𝑂jets))
# compute projection onto the plot variables (project_reachset doesn't currently
# work for a limitation in the hybrid solve, which does not preserve the type of opC)

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
@with_kw mutable struct Battery
    capacity_power::Float64 = 0.0
    capacity_energy::Float64 = 0.0
    soc_init::Float64 = 0.5
    # power::Float64 = 0.0, making the assumption that power is 1/4th of the energy
    efficiency_dch::Float64 = 0.9
    efficiency_ch::Float64 = 0.9
    CAPEXˢ::Float64 = 0.0 #per kWh energy capacity
end

@with_kw mutable struct PV
    period::Period = Period()
    load_factor_PV::Array{Float64,2} = zeros((period.number_of_days, period.number_of_t_steps))
    P_max::Float64 = 0.0
    CAPEXᴾⱽ::Float64 = 0.0 #per kW capacity
end


text:
# Note: the path is IRRELEVANT OF cur working dir
include("../utils/test_module/mod0.jl")

# a. access available by simple `include`
canBeImportedBySimpleInclude()
#: hola! I am available by simple `include(mylocation)`

# b. not accessible by simple `include`
# specificFuction()
#: ERROR: LoadError: UndefVarError: specificFuction not defined

# c. access the not accessible
SpecificModule.specificFuction()
#: Ima specific func!

# d. access not accessible var
println(SpecificModule.specific_var)
#: 007


text:
"""
Simple evaluators implement a single objective function, and use the genotype
of a species when performing the evaluation (as opposed to using some later
stage of development).
"""
type SimpleEvaluator <: Evaluator
  evaluator::Function
  stage::AbstractString
  threads::Int
end

"""
Provides a definition for a simple evaluator.
"""
type SimpleEvaluatorDefinition <: EvaluatorDefinition
  evaluator::Function
  stage::AbstractString
  threads::Int
end

"""
Composes a simple evaluator from its definition.
"""
function compose!(e::SimpleEvaluatorDefinition, p::Population)
  if e.stage == ""
    e.stage = p.demes[1].species.genotype.label
  end
  SimpleEvaluator(e.evaluator, e.stage, e.threads)
end

"""
Simple evaluators implement a single objective function, and use the genotype
of a species when performing the evaluation (as opposed to using some later
stage of development).

"""
simple(f::Function; stage = "", threads = 1) =
  SimpleEvaluatorDefinition(f, stage, max(1, threads))

"""
Evaluates a provided phenome according to a simple evaluation method.
"""
function evaluate{T}(ev::SimpleEvaluator, sch::FitnessScheme, phenome::T)
  println("Called simple evaluator")
  ev.evaluator(sch, phenome)
end


text:
using GeometricEquations
using Requires

import ..Diagnostics: compute_invariant_error

export lotka_volterra_2d_ode,  lotka_volterra_2d_dae,
       lotka_volterra_2d_pode, lotka_volterra_2d_pdae,
       lotka_volterra_2d_iode, lotka_volterra_2d_idae,
       lotka_volterra_2d_lode, lotka_volterra_2d_ldae,
       lotka_volterra_2d_hode, lotka_volterra_2d_hdae,
       lotka_volterra_2d_dg,   lotka_volterra_2d_slrk,
       lotka_volterra_2d_idae_spark

export lotka_volterra_2d_ode_poincare_invariant_1st,
       lotka_volterra_2d_iode_poincare_invariant_1st

export compute_energy_error

const Δt = 0.1
const nt = 1000

const parameters = (a₁=-1.0, a₂=-1.0, b₁=1.0, b₂=2.0)
const reference_solution = [0.34520671417508214, 1.1140647339845946]

const q₀ = [2.0, 1.0]
const v₀=[v₁(0, q₀, parameters), v₂(0, q₀, parameters)]

function f_loop(s)
   rx = 0.2
   ry = 0.3
   x0 = 1.0
   y0 = 1.0

   xs = x0 + rx*cos(2π*s)
   ys = y0 + ry*sin(2π*s)

   qs = [xs, ys]

   return qs
end

function f_loop(i, n)
   f_loop(i/n)
end

function initial_conditions_loop(n)
   q₀ = zeros(2, n)

   for i in axes(q₀,2)
       q₀[:,i] .= f_loop(i, n)
   end

   return q₀
end

compute_energy_error(t,q,params) = compute_invariant_error(t,q, (t,q) -> hamiltonian(t,q,params))

"Creates an ODE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_ode(q₀=q₀; params=parameters)
    ODE(lotka_volterra_2d_v, q₀; parameters=params, invariants=(h=hamiltonian,))
end

"Creates a Hamiltonian ODE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_hode(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
    HODE(lotka_volterra_2d_v, lotka_volterra_2d_f, hamiltonian_pode, q₀, p₀;
         parameters=params)
end

"Creates an implicit ODE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_iode(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
    IODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
         lotka_volterra_2d_g, q₀, p₀;
         parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_2d_v)
end

"Creates a partitioned ODE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_pode(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
    PODE(lotka_volterra_2d_v, lotka_volterra_2d_f, q₀, p₀;
         parameters=params, invariants=(h=hamiltonian_pode,))
end

"Creates a variational ODE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_lode(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
    LODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
         lotka_volterra_2d_g, lagrangian, lotka_volterra_2d_ω, q₀, p₀;
         parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_2d_v)
end

"Creates a DAE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_dae(q₀=vcat(q₀,v₀), λ₀=zero(q₀); params=parameters)
    DAE(lotka_volterra_2d_v_dae, lotka_volterra_2d_u_dae, lotka_volterra_2d_ϕ_dae, q₀, λ₀;
        parameters=params, invariants=(h=hamiltonian,))
end

"Creates a Hamiltonian DAE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_hdae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀); params=parameters)
    HDAE(lotka_volterra_2d_v, lotka_volterra_2d_f,
         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
         lotka_volterra_2d_u̅, lotka_volterra_2d_g̅, lotka_volterra_2d_ψ,
         hamiltonian_pode, q₀, p₀, λ₀; parameters=params)
end

"Creates an implicit DAE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_idae(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
    IDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
         q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
         v̄=lotka_volterra_2d_v)
end

"Creates an implicit DAE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_idae_spark(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
    IDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f_ham,
         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
         q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
         v̄=lotka_volterra_2d_v, f̄=lotka_volterra_2d_f)
end

"Creates a partitioned DAE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_pdae(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
    PDAE(lotka_volterra_2d_v_ham, lotka_volterra_2d_f_ham,
         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
         q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_pode,),
         v̄=lotka_volterra_2d_v, f̄=lotka_volterra_2d_f)
end

"Creates a variational DAE object for the Lotka-Volterra 2D model."
function lotka_volterra_2d_ldae(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
    LDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f_ham,
         lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
         lotka_volterra_2d_u̅, lotka_volterra_2d_g̅, lotka_volterra_2d_ψ_lode,
         lagrangian, lotka_volterra_2d_ω,
         q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
         v̄=lotka_volterra_2d_v, f̄=lotka_volterra_2d_f)
end

"Creates a variational DAE object for the Lotka-Volterra 2D model for use with SLRK integrators."
function lotka_volterra_2d_slrk(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀), λ₀=zero(q₀); params=parameters)
    LDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
            lotka_volterra_2d_u, lotka_volterra_2d_g, lotka_volterra_2d_ϕ,
            lotka_volterra_2d_u̅, lotka_volterra_2d_g̅, lotka_volterra_2d_ψ,
            lagrangian, lotka_volterra_2d_ω,
            q₀, p₀, λ₀; parameters=params, invariants=(h=hamiltonian_iode,),
            v̄=lotka_volterra_2d_v, f̄=lotka_volterra_2d_f)
end

"Creates an implicit ODE object for the Lotka-Volterra 2D model for use with DG integrators."
function lotka_volterra_2d_dg(q₀=q₀, p₀=lotka_volterra_2d_pᵢ(q₀); params=parameters)
    IODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f, lotka_volterra_2d_g,
         q₀, p₀; parameters=params, invariants=(h=hamiltonian_iode,), v̄=lotka_volterra_2d_v)
end

function lotka_volterra_2d_loop_ode(n)
   lotka_volterra_2d_ode(initial_conditions_loop(n))
end

function lotka_volterra_2d_loop_iode(n)
   lotka_volterra_2d_iode(initial_conditions_loop(n))
end

function __init__()
    @require PoincareInvariants = "26663084-47d3-540f-bd97-40ca743aafa4" begin

        function lotka_volterra_2d_ode_poincare_invariant_1st(Δt, nloop, ntime, nsave, DT=Float64)
            PoincareInvariant1st(lotka_volterra_2d_ode, f_loop, ϑ, Δt, 2, nloop, ntime, nsave, DT)
        end

        function lotka_volterra_2d_iode_poincare_invariant_1st(Δt, nloop, ntime, nsave, DT=Float64)
            PoincareInvariant1st(lotka_volterra_2d_iode, f_loop, ϑ, Δt, 2, nloop, ntime, nsave, DT)
        end

    end
end


text:
# First generate some random data and test functions in Base on it
const NVALS_complex = 1000

const input_complex = Dict(
    t=>[ (randindomain(t, NVALS_complex, domain),) for (_, _, domain) in base_unary_complex ]
        for t in (ComplexF32, ComplexF64)
)

const fns_complex = [x[1:2] for x in base_unary_complex]

# output = Dict(
#     t=>[ fns[i](input[t][i]...) for i = 1:length(fns) ]
#         for t in (ComplexF32, ComplexF64)
# )

@testset "Definitions and Comparison with Base for Complex" begin
  for t in (ComplexF32, ComplexF64), i = 1:length(fns_complex)
    inp = input_complex[t][i]
    mod, fn = fns_complex[i]
    base_fn = getproperty(mod, fn)
    vml_fn = getproperty(IntelVectorMath, fn)
    vml_fn! = getproperty(IntelVectorMath, Symbol(fn, "!"))

    Test.@test parentmodule(vml_fn) == IntelVectorMath

    # Test.test_approx_eq(output[t][i], fn(input[t][i]...), "Base $t $fn", "IntelVectorMath $t $fn")
    baseres = base_fn.(inp...)
    Test.@test vml_fn(inp...) ≈ base_fn.(inp...)

    if inp == 1
      if fn != :abs && fn != :angle
        vml_fn!(inp[1])
        Test.@test inp[1] ≈ baseres
      end
    elseif length(inp) == 2
      out = similar(inp[1])
      vml_fn!(out, inp...)
      Test.@test out ≈ baseres
    end

  end

end