text:
const LOG2π = log(2.0π)
const INIT_V = 1.0

struct RELBrownianTree <: AbstractTreeDataBlock
	ind::Int
	n_chars::Int
	llh::Vector{Float64}
end

struct RELBrownianNode <: AbstractNodeDataBlock
	treeblock::RELBrownianTree
	xprune::Vector{Float64}
	llh::Vector{Float64}
end

mutable struct RELBrownianBranch <: AbstractBranchDataBlock
	treeblock::RELBrownianTree
	v::Float64
	vprune::Float64
	δv::Float64
end

#=

	Model initialisators

=#

function init_model!(
	m::RELBrownianTree,
	br::AbstractBranch,
	copylength::Bool=true;
	i::Int=0
	)
	init_v = copylength ? br.length : INIT_V
	mm = RELBrownianBranch(m, init_v, init_v, 0.0)

	if i == 0
		push!(br.datablocks, mm)
	else
		br.datablocks[i] = mm
	end

	return nothing
end

function init_model!(
	m::RELBrownianTree,
	p::AbstractNode,
	x::SpeciesDataMatrix{Float64};
	i::Int=0
	)
	if p.species > 0 && p.species ∈ x.dir
		xprune = x[p.species, :]
	else
		xprune = fill(NaN, m.n_chars)
	end

	mm = RELBrownianNode(m, xprune, fill(NaN, m.n_chars))

	if i == 0
		push!(p.datablocks, mm)
	else
		p.datablocks[i] = mm
	end

	return nothing
end

function _init_model!(
	m::RELBrownianTree,
	p::AbstractNode,
	q::AbstractNode,
	x::SpeciesDataMatrix{Float64},
	copylength::Bool=true;
	i::Int=0
	)

	for link in q.links
		link.to == p && continue
		_init_model!(m, q, link.to, x, copylength, i=i)
		init_model!(m, link.branch, copylength, i=i)
	end

	init_model!(m, q, x)

	return nothing
end

function init_model!(
	tree::AbstractTree,
	x::SpeciesDataMatrix{Float64},
	copylength::Bool=true;
	i::Int=0
	)

	ii = i == 0 ? lastindex(tree.datablocks) + 1 : i
	k = size(x, 2)
	m = RELBrownianTree(ii, k, fill(NaN, k))

	for link in tree.origin.links
		_init_model!(m, tree.origin, link.to, x, copylength, i=i)
		init_model!(m, link.branch, copylength, i=i)
	end
	init_model!(m, tree.origin, x)

	if i == 0
		push!(tree.datablocks, m)
	else
		@inbounds tree.datablocks[i] = m
	end

	return ii
end

#=

	Compute likelihood

=#

"""
Compute the log-likelihood of a node based on its two children
"""
@inline function _calc_llh_2c(
	p::AbstractNode,
	q::AbstractNode,
	i::Int
	)
	childlinks = filter(x -> x.to ≠ p, q.links)

	x₁ = childlinks[1].to.datablocks[i].xprune
	v₁ = childlinks[1].branch.datablocks[i].vprune
	x₂ = childlinks[2].to.datablocks[i].xprune
	v₂ = childlinks[2].branch.datablocks[i].vprune

	# Following Felsenstein (1981), Eqn. 9:
	return @. -0.5 * ((log(v₁ + v₂) + LOG2π) + ((x₁ - x₂)^2) / (v₁ + v₂))
end

"""
Compute the log-likelihood of a node based on its three children
"""
@inline function _calc_llh_3c(p::AbstractNode, i::Int)
	x₁ = p.links[1].to.datablocks[i].xprune
	v₁ = p.links[1].branch.datablocks[i].vprune
	x₂ = p.links[2].to.datablocks[i].xprune
	v₂ = p.links[2].branch.datablocks[i].vprune
	x₃ = p.links[3].to.datablocks[i].xprune
	v₃ = p.links[3].branch.datablocks[i].vprune
	Σv₁₂ = v₁ + v₂

	# Following Felsenstein (1981), Eqn. A1-1:
	llh = @. log(Σv₁₂) + ((x₁ - x₂)^2) / Σv₁₂
	@. llh += log(v₃ + v₁ * v₂ / Σv₁₂)
	@. llh += ((x₃ - (v₂ * x₁ + v₁ * x₂) / Σv₁₂)^2) / (v₃ + v₁ * v₂ / Σv₁₂)
	@. llh *= -0.5
	@. llh -= LOG2π

	return llh
end

function _calc_llh!(
	p::AbstractNode,
	q::AbstractNode,
	i::Int,
	accllh::Vector{Float64};
	prune::Bool=true
	)
	for link in q.links
		link.to == p && continue
		istip(link.to) && continue
		accllh .= _calc_llh!(q, link.to, i, accllh)
	end
	prune && nodeprune!(p, q, i)
	@inbounds q.datablocks[i].llh .= _calc_llh_2c(p, q, i)
	@inbounds accllh .+= q.datablocks[i].llh

	return accllh
end

function calc_llh!(tree::AbstractTree, i::Int; prune::Bool=true)
	k = tree.datablocks[i].n_chars
	accllh = zeros(k)
	for link in tree.origin.links
		istip(link.to) && continue
		accllh .= _calc_llh!(tree.origin, link.to, i, accllh, prune=prune)
	end

	(2 ≤ n_neighbour(tree.origin) ≤ 3) ||
		throw(Phylodendron2.InvalidTopology("the origin of the tree must be of degree 2 or 3"))

	if n_neighbour(tree.origin) > 2
		@inbounds tree.origin.datablocks[i].llh .=
			_calc_llh_3c(tree.origin, i)
	else
		@inbounds tree.origin.datablocks[i].llh .=
			_calc_llh_2c(tree.origin, tree.origin, i)
	end
	@inbounds tree.datablocks[i].llh .= accllh .+ tree.origin.datablocks[i].llh

	return sum(tree.datablocks[i].llh)
end

"""
    phylip_llh(m::RELBrownianTree)

Return the log-likelihood of the model `m` as computed by Phylip.

Phylip computes log-likelihoods of trees differently because it ommits a constant term from the likelihoods of every node. The term is (n - 1)/2 * k * ln(2π), where k is the number of characters and n the number of tips.
"""
function phylip_llh(tree::AbstractTree, i::Int)::Float64
	m = tree.datablocks[i]::RELBrownianTree

	return sum(m.llh) + (n_tip(tree) - 1) / 2 * m.n_chars * LOG2π
end

#=

	Prune

=#

function nodeprune!(
	p::AbstractNode,
	q::AbstractNode,
	i::Int,
	calc_llh::Bool=false
	)
	if istip(q)
		br = getbranch(p, q)
		br.datablocks[i].vprune = br.datablocks[i].v

		return nothing
	end
	child₁, child₂ = filter(x -> x ≠ p, neighbours(q))
	br₀ = getbranch(p, q)
	br₁ = getbranch(q, child₁)
	br₂ = getbranch(q, child₂)

	x₁ = child₁.datablocks[i].xprune
	x₂ = child₂.datablocks[i].xprune
	v₁ = br₁.datablocks[i].vprune
	v₂ = br₂.datablocks[i].vprune
	Σv = v₁ + v₂

	# Following Felsenstein (1981), Eqn. 12:
	@inbounds @. q.datablocks[i].xprune = (x₁ * v₂ + x₂ * v₁) / Σv

	# Following Felsenstein (1981), Eqn. 13:
    @inbounds br₀.datablocks[i].δv = v₁ * v₂ / Σv
	@inbounds br₀.datablocks[i].vprune =
		br₀.datablocks[i].v + br₀.datablocks[i].δv

	return nothing
end

function subtreeprune!(
	p::AbstractNode,
	q::AbstractNode,
	i::Int;
	calc_llh::Bool=false
	)
	for link in q.links
		link.to == p && continue
		subtreeprune!(q, link.to, i)
	end

	nodeprune!(p, q, i)

	return nothing
end

function treeprune!(tree::AbstractTree, i::Int)
	for link in tree.origin.links
		subtreeprune!(tree.origin, link.to, i)
	end

	return nothing
end

#=

	Optimise branch lengths

=#

function optimv_3c!(p::AbstractNode, i::Int)
	child₁, child₂, child₃ = neighbours(p)

	br₁ = getbranch(p, child₁)
	br₂ = getbranch(p, child₂)
	br₃ = getbranch(p, child₃)

	x₁ = child₁.datablocks[i].xprune
	x₂ = child₂.datablocks[i].xprune
	x₃ = child₃.datablocks[i].xprune

	k = length(x₁)

	kv̂₁ = sum((x₁ .- x₂) .* (x₁ .- x₃))
    kv̂₂ = sum((x₂ .- x₁) .* (x₂ .- x₃))
    kv̂₃ = sum((x₃ .- x₁) .* (x₃ .- x₂))

	if kv̂₁ < 0.0
        kv̂₁ = 0.0
        kv̂₂ = sum((x₁ .- x₂).^2)
        kv̂₃ = sum((x₁ .- x₃).^2)
    elseif kv̂₂ < 0.0
        kv̂₁ = sum((x₂ .- x₁).^2)
        kv̂₂ = 0.0
        kv̂₃ = sum((x₂ .- x₃).^2)
    elseif kv̂₃ < 0.0
        kv̂₁ = sum((x₃ .- x₁).^2)
        kv̂₂ = sum((x₃ .- x₂).^2)
        kv̂₃ = 0.0
    end

    @inbounds br₁.datablocks[i].vprune = kv̂₁ / k
    @inbounds br₂.datablocks[i].vprune = kv̂₂ / k
    @inbounds br₃.datablocks[i].vprune = kv̂₃ / k

	@inbounds br₁.datablocks[i].vprune =
		max(br₁.datablocks[i].vprune, br₁.datablocks[i].δv)
	@inbounds br₂.datablocks[i].vprune =
        max(br₂.datablocks[i].vprune, br₂.datablocks[i].δv)
	@inbounds br₃.datablocks[i].vprune =
        max(br₃.datablocks[i].vprune, br₃.datablocks[i].δv)

	@inbounds br₁.datablocks[i].v =
		br₁.datablocks[i].vprune - br₁.datablocks[i].δv
	@inbounds br₂.datablocks[i].v =
		br₂.datablocks[i].vprune - br₂.datablocks[i].δv
	@inbounds br₃.datablocks[i].v =
		br₃.datablocks[i].vprune - br₃.datablocks[i].δv

	return nothing
end

function optimv!(tree::AbstractTree, i::Int; niter = 5)
	for _ in 1:niter
		@assert n_neighbour(tree.origin) == 3 "the origin node must have three neighbours."
		treeprune!(tree, i)
		optimv_3c!(tree.origin, i)
		p = tree.origin
		for q in preorder_vector(tree.origin)[2:end]
			istip(q) && continue
			#=
            The following loop works as a shortcut. Instead of prunning the entire tree anew, it only prunes the neighbours of the nodes along the path between the current node (`q`) and the node from the previous iteration (`p`) (Felsenstein 1981, p. 1238).

            Finding the path between the two nodes is itself relatively costly, so this is likely only an improvement in big trees.
			=#
			for r in nodepath(p, q)[2:end]
				for link in r.links
					nodeprune!(r, link.to, i)
				end
			end
			optimv_3c!(q, i)
			p = q
		end
	end

	return nothing
end

"""
	v2brlength!(tree::AbstractTree, m::RELBrownianTree)

Set the branch lengths of a tree from the branch length parameters (v) of a Brownian model.
"""
function v2brlength!(tree::AbstractTree, i::Int)
	@assert tree.datablocks[i] isa RELBrownianTree "datablock " * string(i) * " does not contain a Brownian model"
	for (p, q) in preorder(tree)[2:end]
		v = getbranch(p, q).datablocks[i].v
		brlength!(p, q, v)
	end

	return nothing
end


text:
"""
    AbstractPrecondition

Abstract type for preconditioners of interval linear systems.
"""
abstract type AbstractPrecondition end

"""
    NoPrecondition <: AbstractPrecondition

Type of the trivial preconditioner which does nothing.

### Notes

- An object of type `NoPrecondition` is a function with method

        (np::NoPrecondition)(A::AbstractMatrix{T},
                             b::AbstractVector{T}) where {T<:Interval}

### Example

```jldoctest
julia> A = [2..4 -2..1; -1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> b = [-2..2, -2..2]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-2, 2]

julia> np = NoPrecondition()
NoPrecondition()

julia> np(A, b)
(Interval{Float64}[[2, 4] [-2, 1]; [-1, 2] [2, 4]], Interval{Float64}[[-2, 2], [-2, 2]])
```
"""
struct NoPrecondition <: AbstractPrecondition end

(np::NoPrecondition)(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T<:Interval} = A, b

"""
    InverseMidpoint <: AbstractPrecondition

Preconditioner that preconditions the linear system ``Ax=b`` with ``A_c^{-1}``,
where ``A_c`` is the midpoint matrix of ``A``.

### Notes

- An object of type `InverseMidpoint` is a function with method

        (imp::InverseMidpoint)(A::AbstractMatrix{T},
                               b::AbstractVector{T}) where {T<:Interval}

### Examples

```jldoctest
julia> A = [2..4 -2..1; -1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> b = [-2..2, -2..2]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-2, 2]

julia> imp = InverseMidpoint()
InverseMidpoint()

julia> imp(A, b)
(Interval{Float64}[[0.594594, 1.40541] [-0.540541, 0.540541]; [-0.540541, 0.540541] [0.594594, 1.40541]], Interval{Float64}[[-0.756757, 0.756757], [-0.756757, 0.756757]])
```
"""
struct InverseMidpoint <: AbstractPrecondition end

function (imp::InverseMidpoint)(A::AbstractMatrix{T},
                                b::AbstractVector{T}) where {T<:Interval}
    R = inv(mid.(A))
    return R*A, R*b
end

"""
    InverseDiagonalMidpoint <: AbstractPrecondition

Preconditioner that preconditions the linear system ``Ax=b`` with the diagonal matrix of
``A_c^{-1}``, where ``A_c`` is the midpoint matrix of ``A``.

### Notes

- An object of type `InverseDiagonalMidpoint` is a function with method

        (idmp::InverseDiagonalMidpoint)(A::AbstractMatrix{T},
                                        b::AbstractVector{T}) where {T<:Interval}

### Example

```jldoctest
julia> A = [2..4 -2..1; -1..2 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-2, 1]
 [-1, 2]   [2, 4]

julia> b = [-2..2, -2..2]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-2, 2]

julia> idmp = InverseDiagonalMidpoint()
InverseDiagonalMidpoint()

julia> idmp(A, b)
(Interval{Float64}[[0.666666, 1.33334] [-0.666667, 0.333334]; [-0.333334, 0.666667] [0.666666, 1.33334]], Interval{Float64}[[-0.666667, 0.666667], [-0.666667, 0.666667]])
```
"""
struct InverseDiagonalMidpoint <: AbstractPrecondition end

function (idmp::InverseDiagonalMidpoint)(A::AbstractMatrix{T},
                                         b::AbstractVector{T}) where {T<:Interval}
    R = inv(Diagonal(mid.(A)))
    return R*A, R*b
end

ext: jl


text:
export L2Loss, Regularization, LogLikeLoss, prior_loss

struct Regularization{L,P} <: DiffEqBase.DECostFunction
  λ::L
  penalty::P
end
Regularization(λ) = Regularization{typeof(λ),typeof(L2Penalty())}(λ,L2Penalty())

function (f::Regularization)(p)
  f.λ*value(f.penalty, p)
end

function prior_loss(prior,p)
  ll = 0.0
  if eltype(prior) <: UnivariateDistribution
    for i in 1:length(prior)
      ll -= logpdf(prior[i],p[i])
    end
  else
    ll -= logpdf(prior,p)
  end
  ll
end

struct L2Loss{T,D,U,W} <: DiffEqBase.DECostFunction
  t::T
  data::D
  differ_weight::U
  data_weight::W
end

function (f::L2Loss)(sol::DiffEqBase.DESolution)
  data = f.data
  weight = f.data_weight
  diff_weight = f.differ_weight
  fill_length = length(f.t)-length(sol)
  for i in 1:fill_length
    if eltype(sol.u) <: Number
      push!(sol.u,Inf)
    else
      push!(sol.u,fill(Inf,size(sol[1])))
    end
  end
  sumsq = 0.0

  if weight == nothing
    @inbounds for i in 2:length(sol)
      for j in 1:length(sol[i])
        sumsq +=(data[j,i] - sol[j,i])^2
      end
      if diff_weight != nothing
          for j in 1:length(sol[i])
            if typeof(diff_weight) <: Real
              sumsq += diff_weight*((data[j,i] - data[j,i-1] - sol[j,i] + sol[j,i-1])^2)
            else
             sumsq += diff_weight[j,i]*((data[j,i] - data[j,i-1] - sol[j,i] + sol[j,i-1])^2)
            end
          end
      end
    end
  else
    @inbounds for i in 2:length(sol)
      if typeof(weight) <: Real
        for j in 1:length(sol[i])
          sumsq = sumsq + ((data[j,i] - sol[j,i])^2)*weight
        end
      else
        for j in 1:length(sol[i])
          sumsq = sumsq + ((data[j,i] - sol[j,i])^2)*weight[j,i]
        end
      end
      if diff_weight != nothing
        for j in 1:length(sol[i])
          if typeof(diff_weight) <: Real
            sumsq += diff_weight*((data[j,i] - data[j,i-1] - sol[j,i] + sol[j,i-1])^2)
          else
            sumsq += diff_weight[j,i]*((data[j,i] - data[j,i-1] - sol[j,i] + sol[j,i-1])^2)
          end
        end
      end
    end
  end
  sumsq
end
L2Loss(t,data;differ_weight=nothing,data_weight=nothing) = L2Loss(t,data,differ_weight,data_weight)

function (f::L2Loss)(sol::DiffEqBase.AbstractMonteCarloSolution)
  mean(f.(sol.u))
end

struct LogLikeLoss{T,D} <: DiffEqBase.DECostFunction
  t::T
  data_distributions::D
  diff_distributions::L where L<:Union{Nothing,D}
  weight
end
LogLikeLoss(t,data_distributions) = LogLikeLoss(t,data_distributions,nothing,nothing)
LogLikeLoss(t,data_distributions,diff_distributions) = LogLikeLoss(t,data_distributions,diff_distributions,1)

function (f::LogLikeLoss)(sol::DESolution)
  distributions = f.data_distributions
  fill_length = length(f.t)-length(sol)
  for i in 1:fill_length
    push!(sol.u,fill(Inf,size(sol[1])))
  end
  ll = 0.0

  if eltype(distributions) <: UnivariateDistribution
    for j in 1:length(f.t), i in 1:length(sol[1][1])
      # i is the number of time points
      # j is the size of the system
      # corresponds to distributions[i,j]
      ll -= logpdf(distributions[i,j],sol[i,j])
    end
  else # MultivariateDistribution
    for j in 1:length(f.t), i in 1:length(sol[1][1])
      # i is the number of time points
      # j is the size of the system
      # corresponds to distributions[i,j]
      ll -= logpdf(distributions[i],sol[i])
    end
  end

  if f.diff_distributions != nothing
    distributions = f.diff_distributions
    diff_data = sol.u[2:end] - sol.u[1:end-1]
    fill_length = length(f.t)-length(diff_data)
    for i in 1:fill_length
      push!(diff_data,fill(Inf,size(sol[1])))
    end
    fdll = 0
    if eltype(distributions) <: UnivariateDistribution
      for j in 1:length(f.t)-1, i in 1:length(diff_data[1])
        fdll -= logpdf(distributions[j,i],diff_data[j][i])
      end
    else
      for j in 1:length(f.t)-1
        fdll -= logpdf(distributions[j],diff_data[j])
      end
    end
    ll += f.weight*fdll
  end

  ll
end

function (f::LogLikeLoss)(sol::DiffEqBase.AbstractMonteCarloSolution)
  distributions = f.data_distributions
  for s in sol
    fill_length = length(f.t)-length(s)
    for i in 1:fill_length
      push!(s.u,fill(Inf,size(s[1])))
    end
  end
  ll = 0.0
  if eltype(distributions) <: UnivariateDistribution
    for j in 1:length(f.t), i in 1:length(sol[1][1])
      # i is the number of time points
      # j is the size of the system
      # corresponds to distributions[i,j]
      vals = [s[i,j] for s in sol]
      ll -= loglikelihood(distributions[i,j],vals)
    end
  else
    for j in 1:length(f.t)
      # i is the number of time points
      # j is the size of the system
      # corresponds to distributions[i,j]
      vals = [s[i,j] for i in 1:length(sol[1][1]), s in sol]
      ll -= loglikelihood(distributions[j],vals)
    end
  end

  if f.diff_distributions != nothing
    distributions = f.diff_distributions
    fdll = 0
    if eltype(distributions) <: UnivariateDistribution
      for j in 2:length(f.t), i in 1:length(sol[1][1])
        vals = [s[i,j] - s[i,j-1] for s in sol]
        fdll -= logpdf(distributions[j-1,i],vals)[1]
      end
    else
      for j in 2:length(f.t)
        vals = [s[i,j] - s[i,j-1] for i in 1:length(sol[1][1]), s in sol]
        fdll -= logpdf(distributions[j-1],vals)[1]
      end
    end
    ll += f.weight*fdll
  end

  ll
end


text:
# TODO: Adapt python model below to new Julia format
# ########### Silicon Sellmeier model for waveguide phase matching calculations
#
# def n_silicon_sym():
#     """This function creates a symbolic representation (using SymPy) of the
#     Sellmeier Equation model for the index of refraction.
#     Equation form is based on Osgood, Panoiu, Dadap, et al.
#     "Engineering nonlinearities in nanoscale optical systems: physics and
# 	applications in dispersion-engineered silicon nanophotonic wires"
# 	Advances in Optics and Photonics Vol. 1, Issue 1, pp. 162-235 (2009)
#     https://doi.org/10.1364/AOP.1.000162
#     valid from 1.2-?? um
#     Thermo-optic coefficients from
#     This model is then exported to other functions that use it and its
#     derivatives to return index, group index and GVD values as a function
#     of temperature and wavelength.
#     Variable units are lm in [um] and T in [deg C]
#     """
#     eps = 11.6858
#     A = 0.939816 # um^2
#     B = 8.10461e-3
#     lm1 = 1.1071 # um^2
#     dn_dT = 0.0e-5 # 1/degK
#     lm, T = sp.symbols('lm T')
#     T0 = 24.5 # reference temperature in [Deg C]
#     n_sym = sp.sqrt( eps  +  A / lm**2 + ( B * lm**2 ) / ( lm**2 - lm1 ) ) # + dn_dT * ( T - T0 )
#     # n_sym = sp.sqrt(   A0  + ( B1 * lm**2 ) / ( lm**2 - C1 ) + ( B2 * lm**2 ) / ( lm**2 - C2 ) ) + dn_dT * ( T - T0 )
#     return lm, T, n_sym
#
# def n_silicon(lm_in,T_in):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the index of refraction of silicon. Equation form is based on
#     "Broadband mid-infrared frequency comb generation in a silicon microresonator"
#     by Luke et al., Optics Letters Vol. 40, Issue 21, pp. 4823-4826 (2015)
#     Variable [units] passed to symbolic equation are lm in [um] and T in [deg C]
#     """
#     lm_um = np.array([Q_(lm_in).to(u.um).magnitude]).flatten()
#     T_C = np.array([Q_(T_in).to(u.degC).magnitude]).flatten()
#     lm, T, n_sym = n_silicon_sym()
#     n = sp.lambdify([lm,T],n_sym,'numpy')
#     output = np.zeros((T_C.size,lm_um.size))
#     for T_idx, TT in enumerate(T_C):
#         output[T_idx,:] = n(lm_um, T_C[T_idx])
#     return output
#
# def n_g_silicon(lm_in,T_in):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the group index of refraction of silicon. Equation form is based on
#     "Broadband mid-infrared frequency comb generation in a silicon microresonator"
#     by Luke et al., Optics Letters Vol. 40, Issue 21, pp. 4823-4826 (2015)
#     Variable [units] passed to symbolic equation are lm in [um] and T in [deg C]
#     """
#     lm_um = Q_(lm_in).to(u.um).magnitude
#     T_C = Q_(T_in).to(u.degC).magnitude
#     lm, T, n_sym = n_silicon_sym()
#     n_sym_prime = sp.diff(n_sym,lm)
#     n_g_sym = n_sym - lm*n_sym_prime
#     n_g = sp.lambdify([lm,T],n_g_sym,'numpy')
#     return n_g(lm_um, T_C)
#
#
# def gvd_silicon(lm_in,T_in):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the group velocity dispersion of silicon. Equation form is based on
#     "Broadband mid-infrared frequency comb generation in a silicon microresonator"
#     by Luke et al., Optics Letters Vol. 40, Issue 21, pp. 4823-4826 (2015)
#     Variable [units] passed to symbolic equation are lm in [um] and T in [deg C]
#     """
#     lm_um = Q_(lm_in).to(u.um).magnitude
#     T_C = Q_(T_in).to(u.degC).magnitude
#     lm, T, n_sym = n_silicon_sym()
#     n_sym_double_prime = sp.diff(n_sym,lm,lm)
#     c = Q_(3e8,'m/s') # unitful definition of speed of light
#     gvd_sym_no_prefactors = (lm**3)*n_sym_double_prime # symbolic gvd without unitful prefactors, to be made unitful below
#     gvd_no_prefactors = sp.lambdify([lm,T],gvd_sym_no_prefactors,'numpy') # numerical gvd without unitful prefactors
#     gvd = (1 / (2 * np.pi * (c**2))) * Q_(gvd_no_prefactors(lm_um,T_C),'um')
#     return gvd.to('fs**2 / mm')


text:
# Build script for Armadillo: C++ library for linear algebra.  It contains a
# library libarmadillo.so that is a wrapper for OpenBLAS, and header files.

using BinaryBuilder

name = "armadillo"
version = v"9.850.1"
sources = [
    ArchiveSource("http://sourceforge.net/projects/arma/files/armadillo-9.850.1.tar.xz",
                  "d4c389b9597a5731500ad7a2656c11a6031757aaaadbcafdea5cc8ac0fd2c01f")
]

script = raw"""
cd ${WORKSPACE}/srcdir/armadillo-*/
mkdir build && cd build

FLAGS=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}
       -DCMAKE_INSTALL_PREFIX=${prefix}
       -DCMAKE_BUILD_TYPE=Release
       -DBUILD_SHARED_LIBS=ON)

if [[ "${nbits}" == 64 ]] && [[ "${target}" != aarch64* ]]; then
    FLAGS+=(-Dopenblas_LIBRARY="${libdir}/libopenblas64_.${dlext}")
    # Force Armadillo's CMake configuration to accept OpenBLAS as a LAPACK
    # replacement.
    FLAGS+=(-DLAPACK_LIBRARY="${libdir}/libopenblas64_.${dlext}")

    SYMB_DEFS=()
    for sym in sasum dasum snrm2 dnrm2 sdot ddot sgemv dgemv cgemv zgemv sgemm dgemm cgemm zgemm ssyrk dsyrk cherk zherk; do
        SYMB_DEFS+=("-D${sym}=${sym}_64")
    done
    if [[ "${target}" == *-apple-* ]] || [[ "${target}" == *-mingw* ]]; then
        FLAGS+=(-DALLOW_OPENBLAS_MACOS=ON)
    fi

    for sym in cgbcon cgbsv cgbsvx cgbtrf cgbtrs cgecon cgees cgeev cgeevx cgehrd cgels cgelsd cgemm cgemv cgeqrf cgesdd cgesv cgesvd cgesvx cgetrf cgetri cgetrs cgges cggev cgtsv cgtsvx cheev cheevd cherk clangb clange clanhe clansy cpbtrf cpocon cposv cposvx cpotrf cpotri cpotrs ctrcon ctrsyl ctrtri ctrtrs cungqr dasum ddot dgbcon dgbsv dgbsvx dgbtrf dgbtrs dgecon dgees dgeev dgeevx dgehrd dgels dgelsd dgemm dgemv dgeqrf dgesdd dgesv dgesvd dgesvx dgetrf dgetri dgetrs dgges dggev dgtsv dgtsvx dlahqr dlangb dlange dlansy dlarnv dnrm2 dorgqr dpbtrf dpocon dposv dposvx dpotrf dpotri dpotrs dstedc dsyev dsyevd dsyrk dtrcon dtrevc dtrsyl dtrtri dtrtrs ilaenv sasum sdot sgbcon sgbsv sgbsvx sgbtrf sgbtrs sgecon sgees sgeev sgeevx sgehrd sgels sgelsd sgemm sgemv sgeqrf sgesdd sgesv sgesvd sgesvx sgetrf sgetri sgetrs sgges sggev sgtsv sgtsvx slahqr slangb slange slansy slarnv snrm2 sorgqr spbtrf spocon sposv sposvx spotrf spotri spotrs sstedc ssyev ssyevd ssyrk strcon strevc strsyl strtri strtrs zgbcon zgbsv zgbsvx zgbtrf zgbtrs zgecon zgees zgeev zgeevx zgehrd zgels zgelsd zgemm zgemv zgeqrf zgesdd zgesv zgesvd zgesvx zgetrf zgetri zgetrs zgges zggev zgtsv zgtsvx zheev zheevd zherk zlangb zlange zlanhe zlansy zpbtrf zpocon zposv zposvx zpotrf zpotri zpotrs ztrcon ztrsyl ztrtri ztrtrs zungqr; do
        SYMB_DEFS+=("-D${sym}=${sym}_64")
    done

    export CXXFLAGS="${SYMB_DEFS[@]}"

    # Force the configuration parameter ARMA_BLAS_LONG to be true, as in our
    # setting 64-bit systems are going to need a 64-bit integer to be used for
    # Armadillo's `blas_int` type.
    sed -i 's|// #define ARMA_BLAS_LONG$|#define ARMA_BLAS_LONG|' ../include/armadillo_bits/config.hpp.cmake
    cat ../include/armadillo_bits/config.hpp.cmake
else
    # Force Armadillo's CMake configuration to accept OpenBLAS as a LAPACK
    # replacement.
    FLAGS+=(-DLAPACK_LIBRARY="${libdir}/libopenblas.${dlext}")
fi

cmake .. "${FLAGS[@]}"
make -j${nproc}
make install

# Armadillo links against a _very_ specific version of OpenBLAS on macOS by
# default:
if [[ ${target} == *apple* ]]; then
    # Figure out what version it probably latched on to:
    OPENBLAS_LINK=$(otool -L ${libdir}/libarmadillo.dylib | grep libopenblas64_ | awk '{ print $1 }')
    install_name_tool -change ${OPENBLAS_LINK} @rpath/libopenblas64_.dylib ${libdir}/libarmadillo.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
platforms = supported_platforms()

# The products that we will ensure are always built.
products = [
    LibraryProduct("libarmadillo", :libarmadillo)
]

# Dependencies that must be installed before this package can be built.
dependencies = [
    Dependency("OpenBLAS_jll")
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function derivative!(df::AbstractArray{<:Number}, f, x::Union{Number,AbstractArray{<:Number}}, fx::AbstractArray{<:Number}, integrator, grad_config)
    alg = unwrap_alg(integrator, true)
    tmp = length(x) # We calculate derivtive for all elements in gradient
    if alg_autodiff(alg)
        T = if standardtag(alg)
          typeof(ForwardDiff.Tag(OrdinaryDiffEqTag(),eltype(df)))
        else
          typeof(ForwardDiff.Tag(f,eltype(df)))
        end

        xdual = Dual{T,eltype(df),1}(convert(eltype(df),x),ForwardDiff.Partials((one(eltype(df)),)))
        f(grad_config,xdual)
        df .= first.(ForwardDiff.partials.(grad_config))
        integrator.destats.nf += 1
    else
        FiniteDiff.finite_difference_gradient!(df, f, x, grad_config, dir = diffdir(integrator))
        fdtype = alg_difftype(alg)
        if fdtype == Val{:forward} || fdtype == Val{:central}
            tmp *= 2
            if eltype(df)<:Complex
              tmp *= 2
            end
        end
        integrator.destats.nf += tmp
    end
    nothing
end

function derivative(f, x::Union{Number,AbstractArray{<:Number}},
                    integrator)
    local d
    tmp = length(x) # We calculate derivtive for all elements in gradient
    alg = unwrap_alg(integrator, true)
    if alg_autodiff(alg)
      integrator.destats.nf += 1
      d = ForwardDiff.derivative(f, x)
    else
      d = FiniteDiff.finite_difference_derivative(f, x, alg_difftype(alg), dir = diffdir(integrator))
      if alg_difftype(alg) === Val{:central} || alg_difftype(alg) === Val{:forward}
          tmp *= 2
      end
      integrator.destats.nf += tmp
      d
    end
end

jacobian_autodiff(f, x, odefun, alg) = (ForwardDiff.derivative(f,x),1, alg)
function jacobian_autodiff(f, x::AbstractArray, odefun, alg)
  jac_prototype = odefun.jac_prototype
  sparsity,colorvec = sparsity_colorvec(odefun,x)
  maxcolor = maximum(colorvec)
  chunk_size = get_chunksize(alg)===Val(0) ? nothing : get_chunksize(alg)
  num_of_chunks = chunk_size===nothing ? Int(ceil(maxcolor / getsize(ForwardDiff.pickchunksize(maxcolor)))) :
                                        Int(ceil(maxcolor / _unwrap_val(chunk_size)))
  (forwarddiff_color_jacobian(f,x,colorvec = colorvec, sparsity = sparsity,
                              jac_prototype = jac_prototype, chunksize=chunk_size),
   num_of_chunks)
end

function _nfcount(N,::Type{diff_type}) where diff_type
  if diff_type===Val{:complex}
    tmp = N
  elseif diff_type===Val{:forward}
    tmp = N + 1
  else
    tmp = 2N
  end
  tmp
end

jacobian_finitediff(f, x, ::Type{diff_type}, dir, colorvec, sparsity, jac_prototype) where {diff_type} =
    (FiniteDiff.finite_difference_derivative(f, x, diff_type, eltype(x), dir = dir),2)
function jacobian_finitediff(f, x::AbstractArray, ::Type{diff_type}, dir, colorvec, sparsity, jac_prototype) where {diff_type}
  f_in = diff_type === Val{:forward} ? f(x) : similar(x)
  ret_eltype = eltype(f_in)
  J = FiniteDiff.finite_difference_jacobian(f, x, diff_type, ret_eltype, f_in,
                                            dir = dir, colorvec = colorvec, sparsity = sparsity, jac_prototype = jac_prototype)
  return J, _nfcount(maximum(colorvec),diff_type)
end
function jacobian(f, x, integrator)
    alg = unwrap_alg(integrator, true)
    local tmp
    if alg_autodiff(alg)
      J, tmp = jacobian_autodiff(f, x, integrator.f, alg)
    else
      jac_prototype = integrator.f.jac_prototype
      sparsity,colorvec = sparsity_colorvec(integrator.f,x)
      dir = diffdir(integrator)
      J, tmp = jacobian_finitediff(f, x, alg_difftype(alg), dir, colorvec, sparsity, jac_prototype)
    end
    integrator.destats.nf += tmp
    J
end

jacobian_finitediff_forward!(J,f,x,jac_config,forwardcache,integrator)=
  (FiniteDiff.finite_difference_jacobian!(J,f,x,jac_config,forwardcache,
    dir=diffdir(integrator));maximum(jac_config.colorvec))
jacobian_finitediff!(J,f,x,jac_config,integrator)=
  (FiniteDiff.finite_difference_jacobian!(J,f,x,jac_config,
    dir=diffdir(integrator));2*maximum(jac_config.colorvec))

function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, fx::AbstractArray{<:Number}, integrator::DiffEqBase.DEIntegrator, jac_config)
    alg = unwrap_alg(integrator, true)
    if alg_autodiff(alg)
      forwarddiff_color_jacobian!(J,f,x,jac_config)
      integrator.destats.nf += 1
    else
      isforward = alg_difftype(alg) === Val{:forward}
      if isforward
        forwardcache = get_tmp_cache(integrator, alg, unwrap_cache(integrator, true))[2]
        f(forwardcache, x)
        integrator.destats.nf += 1
        tmp=jacobian_finitediff_forward!(J, f, x, jac_config, forwardcache, integrator)
      else # not forward difference
        tmp=jacobian_finitediff!(J, f, x, jac_config, integrator)
      end
      integrator.destats.nf += tmp
    end
    nothing
end

function DiffEqBase.build_jac_config(alg,f::F1,uf::F2,du1,uprev,u,tmp,du2,::Val{transform}=Val(true)) where {transform,F1,F2}
  haslinsolve = hasfield(typeof(alg),:linsolve)

  if !DiffEqBase.has_jac(f) && # No Jacobian if has analytical solution
    (transform || !DiffEqBase.has_Wfact(f)) && # No Jacobian if has_Wfact and Wfact is the one that's used
    (!transform || !DiffEqBase.has_Wfact_t(f)) && # No Jacobian has_Wfact and Wfact_t is the one that's used
    ((concrete_jac(alg) === nothing && (!haslinsolve || (haslinsolve && # No Jacobian if linsolve doesn't want it
    (alg.linsolve === nothing || LinearSolve.needs_concrete_A(alg.linsolve))))) ||
      (concrete_jac(alg) !== nothing && concrete_jac(alg))) # Jacobian if explicitly asked for
    jac_prototype = f.jac_prototype
    sparsity,colorvec = sparsity_colorvec(f,u)
    if alg_autodiff(alg)
      _chunksize = get_chunksize(alg)===Val(0) ? nothing : get_chunksize(alg) # SparseDiffEq uses different convection...

      T = if standardtag(alg)
        typeof(ForwardDiff.Tag(OrdinaryDiffEqTag(),eltype(u)))
      else
        typeof(ForwardDiff.Tag(uf,eltype(u)))
      end
      jac_config = ForwardColorJacCache(uf,uprev,_chunksize;colorvec=colorvec,sparsity=sparsity,tag=T)
    else
      if alg_difftype(alg) !== Val{:complex}
        jac_config = FiniteDiff.JacobianCache(tmp,du1,du2,alg_difftype(alg),colorvec=colorvec,sparsity=sparsity)
      else
        jac_config = FiniteDiff.JacobianCache(Complex{eltype(tmp)}.(tmp),Complex{eltype(du1)}.(du1),nothing,alg_difftype(alg),eltype(u),colorvec=colorvec,sparsity=sparsity)
      end
    end
  else
    jac_config = nothing
  end
  jac_config
end

get_chunksize(jac_config::ForwardDiff.JacobianConfig{T,V,N,D}) where {T,V,N,D} = Val(N) # don't degrade compile time information to runtime information

function DiffEqBase.resize_jac_config!(jac_config::SparseDiffTools.ForwardColorJacCache, i)
  resize!(jac_config.fx, i)
  resize!(jac_config.dx, i)
  resize!(jac_config.t, i)
  ps = SparseDiffTools.adapt.(DiffEqBase.parameterless_type(jac_config.dx),
                 SparseDiffTools.generate_chunked_partials(jac_config.dx,
                 1:length(jac_config.dx),Val(ForwardDiff.npartials(jac_config.t[1]))))
  resize!(jac_config.p, length(ps))
  jac_config.p .= ps
end

function DiffEqBase.resize_jac_config!(jac_config::FiniteDiff.JacobianCache, i)
  resize!(jac_config, i)
  jac_config
end

function resize_grad_config!(grad_config::AbstractArray, i)
  resize!(grad_config, i)
  grad_config
end

function resize_grad_config!(grad_config::ForwardDiff.DerivativeConfig, i)
  resize!(grad_config.duals, i)
  grad_config
end

function resize_grad_config!(grad_config::FiniteDiff.GradientCache, i)
  @unpack fx, c1, c2 = grad_config
  fx !== nothing && resize!(fx, i)
  c1 !== nothing && resize!(c1, i)
  c2 !== nothing && resize!(c2, i)
  grad_config
end

function build_grad_config(alg,f::F1,tf::F2,du1,t) where {F1,F2}
  if !DiffEqBase.has_tgrad(f)
    if alg_autodiff(alg)

      T = if standardtag(alg)
        typeof(ForwardDiff.Tag(OrdinaryDiffEqTag(),eltype(du1)))
      else
        typeof(ForwardDiff.Tag(f,eltype(du1)))
      end

      if du1 isa Array
        dualt = Dual{T,eltype(du1),1}(first(du1)*t, ForwardDiff.Partials((one(eltype(du1)),)))
        grad_config = similar(du1,typeof(dualt))
        fill!(grad_config,false)
      else
        grad_config = ArrayInterface.restructure(du1,Dual{T,eltype(du1),1}.(du1, (ForwardDiff.Partials((one(eltype(du1)),)),)) .* false)
      end
    else
      grad_config = FiniteDiff.GradientCache(du1,t,alg_difftype(alg))
    end
  else
    grad_config = nothing
  end
  grad_config
end

function sparsity_colorvec(f,x)
  sparsity = f.sparsity
  colorvec = DiffEqBase.has_colorvec(f) ? f.colorvec :
              (isnothing(sparsity) ? (1:length(x)) : matrix_colors(sparsity))
  sparsity,colorvec
end


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

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

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

add_assert_scope(:NfOrdQuoRing)

set_assert_level(:NfOrdQuoRing, 0)

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

elem_type(::Type{NfOrdQuoRing}) = NfOrdQuoRingElem

elem_type(::NfOrdQuoRing) = NfOrdQuoRingElem

base_ring(Q::NfOrdQuoRing) = Q.base_ring

ideal(Q::NfOrdQuoRing) = Q.ideal

basis_mat(Q::NfOrdQuoRing) = Q.basis_mat

parent(x::NfOrdQuoRingElem) = x.parent

parent_type(::Type{NfOrdQuoRingElem}) = NfOrdQuoRing

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

hash(x::NfOrdQuoRingElem, h::UInt) = hash(x.elem, h)

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

needs_parentheses(::NfOrdQuoRingElem) = true

isnegative(::NfOrdQuoRingElem) = false

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

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

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

Base.deepcopy_internal(x::NfOrdQuoRingElem, dict::ObjectIdDict) =
        NfOrdQuoRingElem(parent(x), Base.deepcopy_internal(x.elem, dict))

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

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

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

function show(io::IO, x::NfOrdQuoRingElem)
  print(io, "$(x.elem)")
end

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

function (O::NfOrdQuoRing)(x::NfOrdElem)
  parent(x) != base_ring(O) && error("Cannot coerce element into the quotient ring")
  return NfOrdQuoRingElem(O, x)
end

function (Q::NfOrdQuoRing)(x::Integer)
  return Q(base_ring(Q)(x))
end

function (Q::NfOrdQuoRing)(x::fmpz)
  return Q(base_ring(Q)(x))
end

################################################################################
#
#  Quotient function
#
# (and standart helpers)
#
################################################################################
doc"""
    quo(O::NfOrd, I::NfOrdIdl) -> NfOrdQuoRing, Map
> The quotient ring $O/I$ as a ring together with the section $M: O/I \to O$.
> The pointwise inverse of $M$ is the canonical projection $O\to O/I$.
"""
function quo(O::NfOrd, I::NfOrdIdl)
  # We should check that I is not zero
  Q = NfOrdQuoRing(O, I)
  f = NfOrdQuoMap(O, Q)
  return Q, f
end

doc"""
    ResidueRing(O::NfOrd, I::NfOrdIdl) -> NfOrdQuoRing
> The quotient ring $O$ modulo $I$ as a new ring.
"""
Nemo.ResidueRing(O::NfOrd, I::NfOrdIdl) = NfOrdQuoRing(O, I)

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

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

function +(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  parent(x) != parent(y) && error("Elements must have same parents")
  return parent(x)(x.elem + y.elem)
end

function -(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  parent(x) != parent(y) && error("Elements must have same parents")
  return parent(x)(x.elem - y.elem)
end

function -(x::NfOrdQuoRingElem)
  return parent(x)(-x.elem)
end

function *(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  parent(x) != parent(y) && error("Elements must have same parents")
  return parent(x)(x.elem * y.elem)
end

function mul!(z::NfOrdQuoRingElem, x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  mul!(z.elem, x.elem, y.elem)
  mod!(z.elem, parent(z))
  return z
end

function *(x::Integer, y::NfOrdQuoRingElem)
  return parent(y)(x * y.elem)
end

*(x::NfOrdQuoRingElem, y::Integer) = y*x

function *(x::fmpz, y::NfOrdQuoRingElem)
  return parent(x)(x * y.elem)
end

*(x::NfOrdQuoRingElem, y::fmpz) = y*x

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

#^(a::NfOrdQuoRingElem, f::Integer) = a^fmpz(f)

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

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

iszero(x::NfOrdQuoRingElem) = iszero(x.elem)

isone(x::NfOrdQuoRingElem) = isone(x.elem)

function one(Q::NfOrdQuoRing)
  return Q(one(Q.base_ring))
end

function zero(Q::NfOrdQuoRing)
  return Q(zero(Q.base_ring))
end

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

==(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem) = x.elem == y.elem

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

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

function isdivisible(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  parent(x) != parent(y) && error("Elements must have same parents")

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

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

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

  V = R.tmp_div
  A = representation_matrix(y.elem)
  B = parent(x).basis_mat

  V[1, 1] = 1

  a = elem_in_basis(x.elem, Val{false})

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

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

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

  hnf_modular_eldiv!(V, minimum(parent(x).ideal))

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

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

  ccall((:fmpz_mat_zero, :libflint), Void, (Ptr{fmpz_mat}, ), &V)

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

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

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

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

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

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

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

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

  return q
end

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

function inv(x::NfOrdQuoRingElem)
  return divexact(one(parent(x)), x)
end

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

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

  U = parent(x).tmp_euc

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

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

  hnf_modular_eldiv!(U, parent(x).ideal.minimum)

  z = fmpz(1)

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

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

  return z
end

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

function divrem(x::NfOrdQuoRingElem, y::NfOrdQuoRingElem)
  #q = rand(parent(x))
  #r = x - q*y
  #e = euclid(y)

  # This should be only one case and don't do the try/catch crap
  # Write a proper _is_divisible function

  #if e == 1
  #  q = x*inv(y)
  #  r = x - q*y
  #  @hassert :NfOrdQuoRing 1 iszero(x - q*y)
  #  @hassert :NfOrdQuoRing 1 euclid(r) < e
  #  return q, r
  #end

  #try q = divexact(x, y)
  #  r = x - q*y
  #  @hassert :NfOrdQuoRing 1 iszero(x - q*y)
  #  @hassert :NfOrdQuoRing 1 euclid(r) < e
  #  return q, r
  #catch
  #end

  b, q = isdivisible(x, y)
  if b
    return q, zero(parent(x))
  end

  e = euclid(y)

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

  @hassert :NfOrdQuoRing 1 euclid(r) < e
end

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

function rand(Q::NfOrdQuoRing)
  A = basis_mat(Q)
  B = basis(base_ring(Q))
  z = rand(fmpz(1):A[1,1]) * B[1]

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

  return Q(z)
end

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

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

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

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

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

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

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

  f = NfOrdQuoMap(O, parent(x))

  return f(I)
end

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

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

  d = degree(O)

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

  g = gcd(x, y)

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

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

  M_I = Q.basis_mat

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

  a = elem_in_basis(Q(O(1)).elem, Val{false})

  V = parent(x).tmp_xxgcd

  V[1, 1] = 1

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

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

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

  #U = V

  hnf_modular_eldiv!(V, minimum(Q.ideal))::fmpz_mat

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

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

  ccall((:fmpz_mat_zero, :libflint), Void, (Ptr{fmpz_mat}, ), &V)

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

function (M::Generic.MatSpace{NfOrdQuoRingElem})(x::Generic.Mat{NfOrdElem})
  z = map(base_ring(M), x.entries)::Array{NfOrdQuoRingElem, 2}
  return M(z)::Generic.Mat{NfOrdQuoRingElem}
end
################################################################################
#
#  Hensel Lifting
#
################################################################################
function lift(R::NmodPolyRing, a::fq_nmod)
  f = R()
  for i=0:degree(parent(a))-1
    setcoeff!(f, i, _get_coeff_raw(a, i))
  end
  return f
end

function (Zx::FmpzPolyRing)(a::nf_elem)
  b = Zx()
  @assert denominator(a) == 1
  for i=0:a.elem_length
    setcoeff!(b, i, numerator(coeff(a, i)))
  end
  return b
end

#missing:
#  assertion on the lifting
#  filtering of the roots (maybe we don't want all?)
#  find all/ one/ few roots
#  break down in various modules:
#  - find powers of ideal
#  - find LLL/HNF basis
#  - lifting?
#  use ResRing(Poly) instead of doing % pgg?
#  an open variant where k is increased until we have a root?

function _hensel(f::Generic.Poly{nf_elem}, p::Int, k::Int; max_roots::Int = degree(f))
  k = max(k, 2)
  #assumes f squarefree
  #assumes constant_coefficient(f) != 0
  ZX, X = FlintZZ["X"]
  Rp = ResidueRing(FlintZZ, p, cached=false)
  Rpt, t = PolynomialRing(Rp, "t", cached=false)
  K = base_ring(f)

  gp = Rpt(K.pol)
  #to avoid embarrasment...
  @assert !iszero(discriminant(gp))

  #find the prime ideal - as I don't want to use orders, this is
  #fun (computing a max order just for this is wasteful)
  #fun fact: if g = prod g_i mod p^k, then P_i^k = <p^k, g_i>
  #so instead of powering, and simplify and such, lets write it down
  lp = factor(gp).fac
  if length(lp) > 1
    g1 = lift(ZX, first(keys(lp)))
    gg = hensel_lift(ZX(K.pol), g1, fmpz(p), k)
  else
    gg = ZX(K.pol)
  end
  # now for all i<= k, <p^i, K(gg)+p^i> is a normal presentation for
  #                                     the prime ideal power.
  #later we'll get the HNF matrix for selected powers as well

  #set up the mod p data:
  #need FiniteField as I need to factor (roots)
  S = FqNmodFiniteField(first(keys(lp)), :z, false)
  ST, T = PolynomialRing(S,"T", cached=false)
  fp = ST([S(Rpt(coeff(f, i))) for i=0:degree(f)])
  rt = roots(fp)
  #we're going to lift the roots - and for efficiency 1/f'(rt) as well:
  fps = derivative(fp)
  irt = [inv(fps(x)) for x= rt]

  # this is in the finite field, but now I need it in the res. ring.
  # in ZX for ease of transport.
  RT = [lift(ZX, lift(Rpt, x)) for x = rt]
  IRT = [lift(ZX, lift(Rpt, x)) for x = irt]

  #the den: ala Kronnecker:
  den = ZX(derivative(K.pol)(gen(K)))
  iden = inv(derivative(K.pol)(gen(K)))

  #for the result, to check for stabilising

  res = nf_elem[K(0) for x = RT]
  rs = nf_elem[]

  #we'll be working in (Z/p^k)[t]/gg

  #an "optimal" lifting chain:
  pr = [k]
  while k>1
    k = div(k+1, 2)
    push!(pr, k)
  end
  pr = reverse(pr)

  #lets start:

  for i=2:length(pr)
    pp = fmpz(p)^pr[i]
    Q = ResidueRing(FlintZZ, pp, cached=false)
    Qt, t = PolynomialRing(Q,"t", cached=false)

    #possibly this should be done with max precision and then adjusted down
    #the poly mod P^??
    fp = [Qt(ZX(coeff(f, k))) for k=0:degree(f)]

    #we need to evaluate fp and fp' at the roots (later)
    #given that we evaluate "by hand" we don't need polynomials

    pgg = Qt(gg) #we'll do the reductions by hand - possibly not optimal

    #the lattice for reco:
    n = degree(K)
    M = zero_matrix(FlintZZ, n, n)
    for j=1:degree(pgg)
      M[j,j] = pp
    end
    pt = t^(degree(pgg)-1)
    for j=degree(pgg)+1:n
      pt = (pt*t) % pgg
      M[j,j] = 1
      for k=0:degree(pt)
        M[j, k+1] = -lift(coeff(pt, k))
      end
    end
    #this is (or should be) the HNF basis for P^??
    M = lll(M)
    Mi, d = pseudo_inv(M)

    for j=1:length(RT)
      if RT[j] == 0
        continue
      end
      #to eval fp and the derivative, we pre-compute the powers of the
      #evaluation point - to save on large products...

      pow = [Qt(1), Qt(RT[j])]
      while length(pow) <= degree(f)+1
        push!(pow, pow[2]*pow[end] % pgg)
      end

      eval_f = sum(pow[k] * fp[k] for k=1:length(fp)) % pgg
      eval_fs = sum(pow[k-1] *(k-1)*fp[k] for k=2:length(fp)) % pgg

      #double lift:
      #IRT = invmod(fp'(rt), p^k)
      # using x -> x(2-xy) to compute the inverse of y
      IRT[j] = lift(ZX, Qt(IRT[j])*(Qt(2-IRT[j]*eval_fs) % pgg) %pgg)
      #RT = rt mod p^k normal Newton
      # using x -> x-fp(x)//fp'(x) = x-fp(x) * IRT
      RT[j] = lift(ZX, Qt(pow[2] - eval_f*Qt(IRT[j])) % pgg)

      #before the reconstruction, we need to scale by den
      cf = lift(ZX, Qt(RT[j]*den) % pgg)
      ve = matrix(FlintZZ, 1, n, [coeff(cf, k) for k=0:n-1])
      _ve = ve*Mi
      mu = matrix(FlintZZ, 1, n,  [ round(_ve[1, k]//d) for k=1:n])
      ve = ve - mu*M
      z = ZX()
      for k=1:n
        setcoeff!(z, k-1, ve[1, k])
      end
      zz = K(z)*iden
      if res[j] == zz || i == length(pr)
        if f(zz) == 0
          push!(rs, zz)
          if length(rs) >= max_roots
            return rs
          end
          RT[j] = ZX(0)
        end
      else
        res[j] = zz
      end
    end
  end
  return rs
end

function ispower(a::Nemo.fq_nmod, m::Int)
  s = size(parent(a))
  if gcd(s-1, m) == 1
    return true, a^invmod(m, s-1)
  end
  St, t = PolynomialRing(parent(a), "t", cached=false)
  f = t^m-a
  rt = roots(f)
  if length(rt) > 0
    return true, rt[1]
  else
    return false, a
  end
end

function has_primitive_root_1(K::Nemo.FqNmodFiniteField, m::Int)
  @assert m>0
  if m == 1
    return K(1)
  end
  if size(K) % m != 1
    return false, K(1)
  end
  if m == 2
    return K(-1)
  end
  fm = factor(m)
  while true
    g = rand(K)
    if iszero(g)
      continue
    end
    if any(x -> isone(g^div(m, x)), keys(fm))
      continue
    end
    return true, g^div(size(K)-1, m)
  end
end

#as above, but for f = x^m-a, so m-th root of a

function _hensel(a::nf_elem, m::Int, p::Int, k::Int; max_roots::Int = m)
  k = max(k, 2)
  #@assert denominator(a) == 1
  #well, actually, denominator(a, maximal_order)==1 would be sufficient, but
  #hard to test...
  ZX, X = FlintZZ["X"]
  Rp = ResidueRing(FlintZZ, p, cached = false)
  Rpt, t = PolynomialRing(Rp, "t", cached = false)
  K = parent(a)

  gp = Rpt(K.pol)
  #to avoid embarrasment...
  @assert !iszero(discriminant(gp))

  #find the prime ideal - as I don't want to use orders, this is
  #fun (computing a max order just for this is wasteful)
  #fun fact: if g = prod g_i mod p^k, then P_i^k = <p^k, g_i>
  #so instead of powering, and simplify and such, lets write it down
  lp = factor(gp).fac
  if length(lp) > 1
    g1 = lift(ZX, first(keys(lp)))
    gg = hensel_lift(ZX(K.pol), g1, fmpz(p), k)
  else
    gg = ZX(K.pol)
  end

  # now for all i<= k, <p^i, K(gg)+p^i> is a normal presentation for
  #                                     the prime ideal power.
  #later we'll get the HNF matrix for selected powers as well

  #set up the mod p data:
  #need FiniteField as I need to factor (roots)
  S = FqNmodFiniteField(first(keys(lp)), :z, false)
  ST, T = PolynomialRing(S, "T", cached=false)
  fp = T^m - S(Rpt(a))
  rt = roots(fp)
  #OK, we're doing s.th. different
#  b = a^(m-1) #this dominates everything!!!, maybe it should be done in the
#              #residue ring as well.
  rt = [inv(x^(m-1)) for x = rt]

  # this is in the finite field, but now I need it in the res. ring.
  # in ZX for ease of transport.
  RT = [lift(ZX, lift(Rpt, x)) for x = rt]

  #the den: ala Kronnecker:
  den = ZX(derivative(K.pol)(gen(K)))
  iden = inv(derivative(K.pol)(gen(K)))

  #for the result, to check for stabilising

  res = nf_elem[K(0) for x = RT]
  rs = nf_elem[]

  #we'll be working in (Z/p^k)[t]/gg

  #an "optimal" lifting chain:
  pr = [k]
  while k>1
    k = div(k+1, 2)
    push!(pr, k)
  end
  pr = reverse(pr)

  #lets start:

  for i=2:length(pr)
    pp = fmpz(p)^pr[i]
    Q = ResidueRing(FlintZZ, pp, cached=false)
    Qt, t = PolynomialRing(Q, "t", cached = false)

    #possibly this should be done with max precision and then adjusted down
    #the poly mod P^??

    pgg = Qt(gg) #we'll do the reductions by hand - possibly not optimal

    #the lattice for reco:
    n = degree(K)
    M = zero_matrix(FlintZZ, n, n)
    for j=1:degree(pgg)
      M[j,j] = pp
    end
    pt = t^(degree(pgg)-1)
    for j=degree(pgg)+1:n
      pt = (pt*t) % pgg
      M[j,j] = 1
      for k=0:degree(pt)
        M[j, k+1] = -lift(coeff(pt, k))
      end
    end
    #this is (or should be) the HNF basis for P^??
    M = lll(M)
    Mi, d = pseudo_inv(M)

    ap = Qt(ZX(a))
    bp = powmod(ap, m-1, pgg)
    minv = invmod(fmpz(m), pp)

    for j=1:length(RT)
      if RT[j] == 0
        continue
      end

      RTjp = Qt(RT[j])
      RT[j] = lift(ZX, (RTjp*(1+minv) - bp*minv* powmod(RTjp, m+1, pgg)) % pgg)

      #before the reconstruction, we need to scale by den and by a
      cf = lift(ZX, (Qt(RT[j]*den) % pgg)*ap % pgg)
      ve = matrix(FlintZZ, 1, n, [coeff(cf, k) for k=0:n-1])
      _ve = ve*Mi
      mu = matrix(FlintZZ, 1, n,  [ round(_ve[1, k]//d) for k=1:n])
      ve = ve - mu*M
      z = ZX()
      for k=1:n
        setcoeff!(z, k-1, ve[1, k])
      end
      zz = K(z)*iden
      if res[j] == zz || i == length(pr)
        if zz^m == a
          push!(rs, zz)
          if length(rs) >= max_roots
            return rs
          end
          RT[j] = ZX(0)
        end
      else
        res[j] = zz
      end
    end
  end
  return rs
end

## Hensel lifting of roots
##todo: redo for equation order using the kronnecker basis (co-different)
##      to avoid denominators
##    : use double iteration to avoid division
#     : use exponent chain to not overshoot (too much) in lifting
#     : common case is f == K.pol. In this case we known a sharp T2-bound
#     : for torsion polynomials, embed torsion units faithully and
#     : lift only one root of maximal order
#mostly done for nf_elem, see next function
#missing: different (better) handling for x^n-a
# possibly based on https://scicomp.stackexchange.com/questions/18846/newton-iteration-for-cube-root-without-division
#
# g = x^-n -a^(n-1), and r = root(g), then r = a^(-(n-1)/n), hence
#                                         ar = a^(-(n-1)/n + 1) = a^(1/n)
# Here, the Newton-iteration is easy: (b = a^(n-1))
# g = x^-n - b, => g' = -nx^(-n-1)
#Newton: r -> r-g(r)/g'(r) = r - (r^(-n) -b)/(-n r^(-n-1))
#                          = r - r/(-n) - b/n r^(n+1)
#                          = r(1+1/n) - b/n r^(n+1)
# no division! and no double iteration
# however, for reconstruction: we don't want r (which is non-integral)
# but ar which is...

function _roots_hensel(f::Generic.Poly{NfOrdElem}, max_roots::Int = degree(f))
  # f must be squarefree
  # I should check that
  O = base_ring(f)
  n = degree(O)
  deg = degree(f)

  # First we find a prime ideal such that f is squarefree modulo P
  # (The discriminant of f has only finitely many divisors).

  p = degree(f)+1

  fder = derivative(f)

  found_prime = false

  # Dummy variable
  Q = NfOrdIdl(O)
  pi_F = NfOrdToFqNmodMor()
  rt = Array{fq_nmod, 1}()
  di = discriminant(O)

  while !found_prime
    p = next_prime(p)

    if di % p == 0
      continue
    end

    lP = prime_decomposition(O, p)

    for P in lP
      F, pi_F = ResidueField(O, P[1])

      fmodP = pi_F(f)

      if !issquarefree(fmodP)
        continue
      end

      rt = roots(fmodP)
      Q = P[1]
      found_prime = true
      break
    end
  end

  #println("prime $p $(Q.gen_two)")

  fmodQ = pi_F(f)

  # compute the lifting exponent a la Friedrich-Fieker

  R = ArbField(64, false)
  (r1, r2) = signature(O)

  bound_root = [ R(0) for i in 1:(r1 + r2) ]

  CC = AcbField(64, false)
  CCt, t = PolynomialRing(CC, "t", cached=false)
  conjugates_of_coeffs = [ conjugates_arb(coeff(f, i), 32) for i in 0:degree(f) ]

  for i in 1:r1
    g = CCt([ conjugates_of_coeffs[j + 1][i] for j in 0:degree(f) ])
    bound_root[i] = roots_upper_bound(g)
  end

  for i in 1:r2
    g = CCt([ conjugates_of_coeffs[j + 1][r1 + i] for j in 0:degree(f) ])
    bound_root[r1 + i] = roots_upper_bound(g)
  end

  ss = _lifting_expo(p, Q.splitting_type[2], order(Q), bound_root)
  s = Int(ss)

  #println(length(lin_factor))
  #println("lifting bound: ", s)

  rts = NfOrdElem[]
  #println("aiming for $max_roots roots")
  for j in 1:length(rt)
  #  println("have $(length(roots)) found, want $max_roots")
    if length(rts) >= max_roots
      break
    end

    zero_mod_Q = rt[j]

    #println(zero_mod_Q)

    #println(0, " ", pi_F\(zero_mod_Q))

    # The following should be a uniformizing element
    Q_pi = Q.gen_two

    @hassert :NfOrdQuoRing 1 fmodQ(zero_mod_Q) == 0

    # This is the first step

    I = Q^2
    R, pi_R = quo(O, I)

    t1 = divexact(pi_R(subst(f, pi_F\zero_mod_Q)), pi_R(Q_pi))
    t2 = -inv(pi_R(subst(fder,pi_F\zero_mod_Q)))
    new_a = pi_R\(pi_R(pi_F\zero_mod_Q) + t1*t2*pi_R(Q_pi))
    #return pi_R(f)

    #println(1, " ", new_a)

    old_a = new_a

    RR, pi_RR = R, pi_R

    # reconstruction
    B = basis_mat(I)
    L = lll(B)

    rhs = MatrixSpace(FlintZZ, degree(O), 1, false)(elem_in_basis(new_a)'')
    lhs = transpose(L)

    X, d = solve_rational(lhs, rhs)

    zz = [ fmpq(BigInt(X[l, 1])//BigInt(d) - round(BigInt(X[l, 1])//BigInt(d))) for l in 1:degree(O)]

    cden = denominator(zz[1])

    for l in 2:degree(O)
      cden = lcm(cden, denominator(zz[l]))
    end

    zz_num = [ numerator(cden*zz[l]) for l in 1:degree(O) ]

    v = matrix(FlintZZ, 1, degree(O), zz_num)

    w = v*L

    # There is no slower function

    reconstructed_new = O(fmpz[ divexact(w[1, l], cden) for l in 1:degree(O) ])
    # end reconstruction

    reconstructed_old = reconstructed_new

    stabilized = 0

    i = 1
    while 2^i < s
      if reconstructed_new == reconstructed_old
        stabilized = stabilized + 1
      else
        stabilized = 0
      end

      if stabilized >= 2
        if subst(f, reconstructed_new) == 0
          #push!(rts, reconstructed_new)
          break
        else
          stabilized = 0
        end
      end

      if i == s
        if iszero(subst(f, reconstructed_new))
          push!(rts, reconstructed_new)
        end
        break
      end

      reconstructed_old = reconstructed_new
      old_a = new_a
      R, pi_R = RR, pi_RR
      I = simplify(I^2)

      # From Q^i -> Q^(i+1)
      # From Q^(2^i) -> Q^(2^i+1)

      RR, pi_RR = quo(O, I)
      #t1 = divexact(pi_RR(subst(f, old_a)), pi_RR(Q_pi)^(2^i))
      #t2 = pi_RR(pi_F\(-inv(pi_F(subst(fder, old_a)))))
      #new_a = pi_RR\(pi_RR(old_a) + t1*t2*pi_RR((Q_pi))^(2^i))
      new_a = old_a - subst(f, old_a) * (pi_R\(divexact(one(R), pi_R(subst(fder, old_a)))))
      new_a = pi_RR\(pi_RR(new_a))
      #println(i + 1, " ", new_a)

      # Try to reconstruct:

      B = basis_mat(I)
      L = lll(B)

      rhs = MatrixSpace(FlintZZ, degree(O), 1, false)(elem_in_basis(new_a)'')
      lhs = transpose(L)

      X, d = solve_rational(lhs, rhs)

      zz = [ fmpq(BigInt(X[l, 1])//BigInt(d) - round(BigInt(X[l, 1])//BigInt(d))) for l in 1:degree(O)]

      cden = denominator(zz[1])

      for l in 2:degree(O)
        cden = lcm(cden, denominator(zz[l]))
      end

      zz_num = [ numerator(cden*zz[l]) for l in 1:degree(O) ]

      v = matrix(FlintZZ, 1, degree(O), zz_num)

      w = v*L

      # There is no slower function

      reconstructed_new = O(fmpz[ divexact(w[1, l], cden) for l in 1:degree(O) ])

      #println(i + 1, " ", reconstructed_new)

      i = i + 1
    end
    if f(reconstructed_new) == 0
      push!(rts, reconstructed_new)
    end
  end
  return rts
end

function _lifting_expo(p::Int, deg_p::Int, O::NfOrd, bnd::Array{arb, 1})
  # compute the lifting exponent a la Friedrich-Fieker
  #bnd has upper bounds on |x^{(i)}| 1<= i <= r1+r2 as arbs
  #we're using a prime ideal above p of intertia degree deg_p

  (c1, c2) = norm_change_const(O)
  r1, r2 = signature(O)
  R = parent(bnd[1])
  bd = R(0)
  n = degree(O)
  #so   |x|_mink  <= c_1 |x|_coeff
  #and  |x|_coeff <= c_2 |x|_mink

  for i in 1:r1
    bd += bnd[i]^2
  end

  for i=1:r2
    bd += 2*bnd[i+r1]^2
  end

  boundt2 = max(bd, R(1))

  #CF: there is a prob, in the paper wrt LLL bounds on |x| or |x|^2....
  boundk = R(n)*log(R(c1)*R(c2)*boundt2*exp((R(n*(n-1))//2 + 2)*log(R(2)))//n)//(2*deg_p*log(R(p)))

  ss = abs_upper_bound(boundk, fmpz)
  return ss
end

function _roots_hensel(f::Generic.Poly{nf_elem}, max_roots::Int = degree(f))
  # f must be squarefree
  # I should check that
  K = base_ring(f)
  if iszero(constant_coefficient(f))
    rs = [K(0)]
    f = div(f, gen(parent(f)))
  else
    rs = nf_elem[]
  end
  if length(rs) >= max_roots
    return rs
  end
  n = degree(K)
  deg = degree(f)

  # First we find a prime ideal such that f is squarefree modulo P
  # (The discriminant of f has only finitely many divisors).

  p = degree(f)+1
  deg_p = 1

  while true
    p = next_prime(p)

    Rp = ResidueRing(FlintZZ, p, cached=false)
    Rpt, t = PolynomialRing(Rp, "t", cached=false)
    gp = Rpt(K.pol)
    if iszero(discriminant(gp))
      continue
    end

    lp = factor(gp).fac
    #set up the mod p data:
    #need FiniteField as I need to factor (roots)
    deg_p = degree(first(keys(lp)))
    S = FqNmodFiniteField(first(keys(lp)), :z, false)
    ST, T = PolynomialRing(S, "T", cached=false)
    fp = ST([S(Rpt(coeff(f, i))) for i=0:degree(f)])
    if !issquarefree(fp)
      continue
    end
    rt = roots(fp)
    if length(rt) == 0
      return nf_elem[]
    end
    break
  end

  #println("prime $p ")

  # compute the lifting exponent a la Friedrich-Fieker

  R = ArbField(64, false)

  (r1, r2) = signature(K)

  # for Kronnecker:
  gsa = derivative(K.pol)(gen(K))
  gsa_con = conjugates_arb(gsa, 32)

  bound_root = [ R(0) for i in 1:(r1 + r2) ]

  CC = AcbField(64, false)
  CCt, t = PolynomialRing(CC, "t", cached=false)
  conjugates_of_coeffs = [ conjugates_arb(coeff(f, i), 32) for i in 0:degree(f) ]

  for i in 1:r1
    g = CCt([ conjugates_of_coeffs[j + 1][i] for j in 0:degree(f) ])
    bound_root[i] = roots_upper_bound(g) * abs(gsa_con[i])
  end

  for i in 1:r2
    g = CCt([ conjugates_of_coeffs[j + 1][r1 + i] for j in 0:degree(f) ])
    bound_root[r1 + i] = roots_upper_bound(g) * abs(gsa_con[r1+i])
  end

  ss = _lifting_expo(p, deg_p, EquationOrder(K), bound_root)

  rts = _hensel(f, p, Int(ss), max_roots = max_roots - length(rs))
  return vcat(rs, rts)
end

function _roots_hensel(a::nf_elem, m::Int, max_roots::Int = m)
  # f must be squarefree
  # I should check that
  K = parent(a)
  n = degree(K)

  # First we find a prime ideal such that f is squarefree modulo P
  # (The discriminant of f has only finitely many divisors).

  p = m+1
  deg_p = 1

  while true
    p = next_prime(p)

    Rp = ResidueRing(FlintZZ, p, cached=false)
    Rpt, t = PolynomialRing(Rp, "t", cached = false)
    gp = Rpt(K.pol)
    if iszero(discriminant(gp))
      continue
    end

    lp = factor(gp).fac
    #set up the mod p data:
    #need FiniteField as I need to factor (roots)
    deg_p = degree(first(keys(lp)))
    S = FqNmodFiniteField(first(keys(lp)), :z, false)
    ST, T = PolynomialRing(S, "T", cached=false)
    fp = T^m - S(Rpt(a))
    if !issquarefree(fp)
      continue
    end
    rt = roots(fp)
    if length(rt) == 0
      return nf_elem[]
    end
    break
  end

  #println("prime $p ")

  # compute the lifting exponent a la Friedrich-Fieker

  R = ArbField(64, false)
  (r1, r2) = signature(K)

  # for Kronnecker:
  gsa = derivative(K.pol)(gen(K))
  gsa_con = conjugates_arb(gsa, 32)
  a_con = conjugates_arb(a, 32)

  bound_root = [ R(0) for i in 1:(r1 + r2) ]

  for i in 1:r1+r2
    bound_root[i] = root(abs(a_con[i]), m) * abs(gsa_con[i])
  end

  ss = _lifting_expo(p, deg_p, EquationOrder(K), bound_root)

  rts = _hensel(a, m, p, Int(ss), max_roots = max_roots)
  return rts
end
#identical to hasroot - which one do we want?
function ispower(a::NfOrdElem, n::Int)
  Ox, x = PolynomialRing(parent(a), "x", cached=false)
  f = x^n - a
  r = _roots_hensel(f, 1)

  if length(r) == 0
    return false, parent(a)(0)
  else
    return true, r[1]
  end
end

function probablity(O::NfOrdQuoRing)
  p = 1.0
  I = O.ideal
  f = factor(norm(I))
  for (P, e) in f
    if valuation(I, P) > 0
      p = p * (1 - 1/Float64(norm(P)))
    end
  end
  return p
end

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

doc"""
***
    group_structure(Q::NfOrdQuoRing) -> GrpAbFinGenSnf

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


text:
@inline zeroforcing(args...) = 0

"""
    regularize_forcing(forcing, field, field_name, model_field_names)

"Regularizes" or "adds information" to user-defined forcing objects that are passed to
model constructors. `regularize_forcing` is called inside `model_forcing`.

We need `regularize_forcing` because it is only until `model_forcing` is called that
the fields (and field locations) of various forcing functions are available. The `field`
can be used to infer the location at which the forcing is applied, or to add a field
dependency to a special forcing object, as for `Relxation`.
"""
regularize_forcing(forcing, field, field_name, model_field_names) = forcing # fallback

"""
    regularize_forcing(forcing::Function, field, field_name, model_field_names)

Wrap `forcing` in a `ContinuousForcing` at the location of `field`.
"""
function regularize_forcing(forcing::Function, field, field_name, model_field_names)
    LX, LY, LZ = location(field)
    return ContinuousForcing{LX, LY, LZ}(forcing)
end

regularize_forcing(::Nothing, field::AbstractField, field_name, model_field_names) = zeroforcing

"""
    model_forcing(model_fields; forcings...)

Return a `NamedTuple` of forcing functions for each field in `model_fields`, wrapping
forcing functions in `ContinuousForcing`s and ensuring that `ContinuousForcing`s are
located correctly for each field.
"""
function model_forcing(model_fields; forcings...)

    model_field_names = keys(model_fields)

    regularized_forcings = Tuple(
        field_name in keys(forcings) ?
            regularize_forcing(forcings[field_name], field, field_name, model_field_names) :
            regularize_forcing(nothing, field, field_name, model_field_names)
        for (field_name, field) in pairs(model_fields)
    )

    specified_forcings = NamedTuple{model_field_names}(regularized_forcings)

    return specified_forcings
end


text:
# this script basially only handles `ARGS`
let distributed = Base.PkgId(Base.UUID("8ba89e20-285c-5b6f-9357-94700520ee1b"), "Distributed")
    if haskey(Base.loaded_modules, distributed) && (Distributed = Base.loaded_modules[distributed]).nprocs() > 1
        Distributed.remotecall_eval(Main, 1:Distributed.nprocs(), :(pushfirst!(LOAD_PATH, joinpath($(@__DIR__), "..", "packages"))))
        using VSCodeServer
        Distributed.remotecall_eval(Main, 1:Distributed.nprocs(), :(popfirst!(LOAD_PATH)))
    else
        pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
        using VSCodeServer
        popfirst!(LOAD_PATH)
    end
end

let
    args = [popfirst!(Base.ARGS) for _ in 1:6]
    # load Revise ?
    if "USE_REVISE=true" in args
        try
            @static if VERSION ≥ v"1.5"
                using Revise
            else
                @eval using Revise
                Revise.async_steal_repl_backend()
            end
        catch err
        end
    end

    atreplinit() do repl
        VSCodeServer.toggle_plot_pane(nothing, (;enable="USE_PLOTPANE=true" in args))
        VSCodeServer.toggle_progress(nothing, (;enable="USE_PROGRESS=true" in args))
    end

    conn_pipeline, telemetry_pipeline = args[1:2]
    VSCodeServer.serve(conn_pipeline; is_dev="DEBUG_MODE=true" in args, crashreporting_pipename=telemetry_pipeline)
end

max_stars_repo_name: xgdgsc/julia-vscode

max_issues_repo_name: xgdgsc/julia-vscode

max_forks_repo_name: xgdgsc/julia-vscode


text:
@testset "Dumb aggregate" begin

    sph1 = GRAYTR.Sphere()
    sph2 = GRAYTR.Sphere(translation(0, -2, 0))

    mat = GRAYTR.Lambert(GRAYTR.SingleLine(532.0, 1.0))

    P1 = GRAYTR.GeometricPrimitive(sph1, mat, nothing, 0)
    P2 = GRAYTR.GeometricPrimitive(sph2, mat, nothing, 0)

    Agg = GRAYTR.DumbAggregate([P1, P2])

    @test typeof(Agg) <: GRAYTR.DumbAggregate

    @test GRAYTR.world_bounds(Agg) == GRAYTR.BoundingBox(Point3(-1,-3,-1), Point3(1,1,1))

    R = Ray(Point3(5, 0, 0), Vector3(-1, 0, 0))
    @test GRAYTR.intersectP(R, Agg)
    isect = GRAYTR.intersect(R, Agg)
    @test isect.tmin ≈ 4.0
    @test isect.geometry.p ≈ Point3(1, 0, 0)
    @test isect.geometry.n ≈ Normal3(1, 0, 0)

    R = Ray(Point3(5, 2, 0), Vector3(-1, 0, 0))
    @test !GRAYTR.intersectP(R, Agg)
    isect = GRAYTR.intersect(R, Agg)
    @test isect == nothing

    R = Ray(Point3(5, 0.9*sqrt(2), 0.9*sqrt(2)), Vector3(-1, 0, 0))
    @test !GRAYTR.intersectP(R, Agg)
    isect = GRAYTR.intersect(R, Agg)
    @test isect == nothing

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# Auxiliary functions
function gcd(a::Integer, b::Integer)
    if a > b
        return gcd(b, a)
    else
        c = b % a
        if c == 0
            return a
            else return gcd(a, c)
        end
    end
end

function gcd(a::Vector{Integer})
    m = a[1]
    for i in 1:length(a)
        m = gcd(m, a[i])
    end
    m
end

lcm(a::Integer, b::Integer) = round(Int, a * b / gcd(a, b))

function lcm(a::AbstractVector{<:Integer})
    m = a[1]
    for i in 1:length(a)
        m = lcm(m, a[i])
    end
m
end

"""
    sieve(n)
Sieve of Eratosthenes, finding all prime numbers up to `n`, see https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
"""
function sieve(n)
    is_prime=trues(n)
        for i in 2:ceil(Int,√n)
            if is_prime[i]
                for j in i^2:i:n
                    if j %i==0
                        is_prime[j]=false
                    end
                end
            end
        end
    return (x for x in 2:n if is_prime[x])
end

function num_of_factors(n)
    n> 1 || return 1
    key=0
    val=0
    s=1
    k=2
    while n!=k
        if n %k==0
        key=k
        val+=1
            n=round(Int,n/k)
        else
        s*=val+1
        val=0
            k+= k==2 ? 1 : 2
        end
    end
    n==key ? val+=1 : val=1
    s*=val+1
end

function is_abundant(n)
    n>=12 || return false
    s=1
    i=2
    while i*i <=n
        a= n%i
        b=div(n,i)
        if a==0
            if i != b
                s+=i+b
            else
                s+=i
            end
        end
        i+=1
        if s > n
            return true
        end
    end
    return false
end

"""
    cyclelen(x::Integer)
length of Reciprocal cycles
"""
function cyclelen(x::Integer)
    rems=zeros(x)
    loc=1
    @assert x>=1
    r=1
    d,r=divrem(r,x)
    while r!=0
        rems[loc]=r
        r*=10
        d,r=divrem(r,x)
        t=findfirst(x->x==r,rems[1:loc]) #第一个相同位置
        loc+=1 # 当前位置
        if !isnothing(t)
            return loc-t
        end
    end
    return 0
end

"""
    primecounter(f::Function)
return the amount of  prime numbers  generated by `f`.
"""
function primecounter(f::Function)
    x=0
    out=f(x)
    cnt=0
    while isprime(out)
        cnt+=1
        x+=1
        out=f(x)
    end
    return cnt
end

function isprime(n::Integer)
    n<= 1 && return false
    if iseven(n) && n!=2
        return false
    end
    if n==3
        return true
    end
    i=3
    while i*i <=n
        if n%i ==0
            return false
        end
        i+=2
    end
    return true
end

"""
    cycles(n)
return circular numbers generated by `n`, eg. 123 gives [123,231,312].
"""
function cycles(n)
    @assert n>=0
    if n <10
        return [n]
    end
    digs=digits(n)
    d=length(digs)
    cyc=zeros(Int,d)
    x=[10^i for i in 0:d-1]
    for i in 1:d
        cyc[i]=sum(x .* digs[vcat(i:d,1:i-1)])
    end
    return cyc
end

function istruncatableprime(n)
    @assert n >10
    if !isprime(n)
        return false
    else
        str=string(n)
        flag=true
        for i in 1:length(str)-1
            a=parse(Int,str[1:i])
            b=parse(Int,str[i+1:end])
            if !isprime(a) || !isprime(b)
                flag=false
                break
            end
        end
        return flag
    end
end


text:
# Concrete algorithm implementations.
include("support/helper.jl")
include("support/resample.jl")
# @suppress_err begin
  include("support/distributions.jl")
# end
include("support/hmc_core.jl")
include("support/adapt.jl")
include("support/init.jl")
include("support/stan-interface.jl")
include("hmcda.jl")
include("nuts.jl")
include("sghmc.jl")
include("sgld.jl")
include("hmc.jl")
include("mh.jl")
include("is.jl")
include("smc.jl")
include("pgibbs.jl")
include("pmmh.jl")
include("ipmcmc.jl")
include("gibbs.jl")

## Fallback functions

# utility funcs for querying sampler information
require_gradient(spl::Sampler) = false
require_particles(spl::Sampler) = false

assume(spl::Sampler, dist::Distribution) =
  error("[assume]: unmanaged inference algorithm: $(typeof(spl))")

observe(spl::Sampler, weight::Float64) =
  error("[observe]: unmanaged inference algorithm: $(typeof(spl))")

## Default definitions for assume, observe, when sampler = nothing.
assume(spl::Nothing, dist::Distribution, vn::VarName, vi::VarInfo) = begin
  if haskey(vi, vn)
    r = vi[vn]
  else
    r = init(dist)
    push!(vi, vn, r, dist, 0)
  end
  # NOTE: The importance weight is not correctly computed here because
  #       r is genereated from some uniform distribution which is different from the prior
  # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))

  r, logpdf_with_trans(dist, r, istrans(vi, vn))
end

assume(spl::Nothing, dists::Vector{T}, vn::VarName, var::Any, vi::VarInfo) where T<:Distribution = begin
  @assert length(dists) == 1 "[assume] Turing only support vectorizing i.i.d distribution"
  dist = dists[1]
  n = size(var)[end]

  vns = map(i -> copybyindex(vn, "[$i]"), 1:n)

  if haskey(vi, vns[1])
    rs = vi[vns]
  else
    rs = init(dist, n)

    if isa(dist, UnivariateDistribution) || isa(dist, MatrixDistribution)
      for i = 1:n
        push!(vi, vns[i], rs[i], dist, 0)
      end
      @assert size(var) == size(rs) "[assume] variable and random number dimension unmatched"
      var = rs
    elseif isa(dist, MultivariateDistribution)
      for i = 1:n
        push!(vi, vns[i], rs[:,i], dist, 0)
      end
      if isa(var, Vector)
        @assert length(var) == size(rs)[2] "[assume] variable and random number dimension unmatched"
        for i = 1:n
          var[i] = rs[:,i]
        end
      elseif isa(var, Matrix)
        @assert size(var) == size(rs) "[assume] variable and random number dimension unmatched"
        var = rs
      else
        error("[Turing] unsupported variable container")
      end
    end
  end

  # acclogp!(vi, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1]))))

  var, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1])))
end

observe(spl::Nothing, dist::Distribution, value::Any, vi::VarInfo) = begin
  vi.num_produce += 1
  @debug "dist = $dist"
  @debug "value = $value"

  # acclogp!(vi, logpdf(dist, value))
  logpdf(dist, value)
end

observe(spl::Nothing, dists::Vector{T}, value::Any, vi::VarInfo) where T<:Distribution = begin
  @assert length(dists) == 1 "[observe] Turing only support vectorizing i.i.d distribution"
  dist = dists[1]
  @assert isa(dist, UnivariateDistribution) || isa(dist, MultivariateDistribution) "[observe] vectorizing matrix distribution is not supported"
  if isa(dist, UnivariateDistribution)  # only univariate distributions support broadcast operation (logpdf.) by Distributions.jl
    # acclogp!(vi, sum(logpdf.(Ref(dist), value)))
    sum(logpdf.(Ref(dist), value))
  else
    # acclogp!(vi, sum(logpdf(dist, value)))
    sum(logpdf(dist, value))
  end
end


text:
"""
    struct PreprocessorOptions <: Any
Hold a pointer to a `clang::PreprocessorOptions` object.
"""
struct PreprocessorOptions
    ptr::CXPreprocessorOptions
end

Base.unsafe_convert(::Type{CXPreprocessorOptions}, x::PreprocessorOptions) = x.ptr
Base.cconvert(::Type{CXPreprocessorOptions}, x::PreprocessorOptions) = x


text:
#=
@TODO use our own UI components for this:
atom-ide-ui is already deprecated, ugly, not fully functional, and and...
=#

handle("datatip") do data
  @destruct [
    word,
    mod || "Main",
    path || "",
    column || 1,
    row || 1,
    startRow || 0,
    context || ""
  ] = data
  datatip(word, mod, path, column, row, startRow, context)
end

function datatip(word, mod, path, column = 1, row = 1, startrow = 0, context = "")
  if isdebugging() && (ddt = JunoDebugger.datatip(word, path, row, column)) !== nothing
    return Dict(:error => false, :strings => ddt)
  end

  ldt = localdatatip(word, column, row, startrow, context)
  isempty(ldt) || return datatip(ldt)

  tdt = globaldatatip(mod, word)
  tdt !== nothing && return Dict(:error => false, :strings => tdt)

  return Dict(:error => true) # nothing hits
end

datatip(dt::Vector{Dict{Symbol, Any}}) = Dict(:error => false, :strings => dt)
datatip(dt::Int) = Dict(:error => false, :line => dt)
datatip(dt::Vector{Int}) = datatip(dt[1])

function localdatatip(word, column, row, startrow, context)
  word = first(split(word, '.')) # always ignore dot accessors
  position = row - startrow
  ls = locals(context, position, column)
  filter!(l -> l.name == word && l.line < position, ls)
  # there should be zero or one element in `ls`
  return localdatatip.(ls, word, startrow)
end

function localdatatip(l, word, startrow)
  verbatim = l.verbatim
  return if verbatim == word # when `word` is an argument or such
    startrow + l.line - 1
  else
    Dict(:type => :snippet, :value => verbatim)
  end
end

function globaldatatip(mod, word)
  docs = @errs getdocs(mod, word)
  docs isa EvalError && return nothing

  # don't show verbose stuff
  docstr = replace(string(docs), nodoc_regex => "")
  occursin(nobinding_regex, docstr) && return nothing

  datatip = []

  val = getfield′(getmodule(mod), word)
  processval!(val, docstr, datatip)

  processdoc!(docs, docstr, datatip)

  ml = methods(val)
  processmethods!(ml, datatip)

  return datatip
end

# adapted from `REPL.summarize(binding::Binding, sig)`
const nodoc_regex = r"^No documentation found."
const nobinding_regex = r"Binding `.*` does not exist.\n"

function processdoc!(docs, docstr, datatip)
  # Separate code blocks from the other markdown texts in order to render them
  # as code snippets in atom-ide-ui's datatip service.
  # Setting up functions to deconstruct each `Markdown.MD.content` into
  # an appropriate markdown string would be more robust.
  texts = split(docstr, codeblock_regex)
  codes = searchcodeblocks(docs)

  pushmarkdown!(datatip, texts[1])
  for (code, text) in zip(codes, texts[2:end])
    pushsnippet!(datatip, code)
    pushmarkdown!(datatip, text)
  end
end

# will match code blocks within a markdown text
const codeblock_regex = r"```((?!```).)*?```"s

# extract only code blocks from Markdown.MD
function searchcodeblocks(docs)
  codeblocks = []
  searchcodeblocks(docs, codeblocks)
  codeblocks
end
function searchcodeblocks(docs, codeblocks)
  for content in docs.content
    if content isa Markdown.Code
      push!(codeblocks, content.code)
    elseif content isa Markdown.MD
      searchcodeblocks(content, codeblocks)
    end
  end
end

processval!(@nospecialize(val), docstr, datatip) = begin
  valstr = @> repr(MIME("text/plain"), val, context = :limit => true) strlimit(1000, " ...")
  occursin(valstr, docstr) || pushsnippet!(datatip, valstr)
end
processval!(val::Function, docstr, datatip) = begin
  valstr = string(val) # this would get rid of the unnecessary module prefix
  occursin(valstr, docstr) || pushsnippet!(datatip, valstr)
end
processval!(::Undefined, docstr, datatip) = nothing

function processmethods!(ml, datatip)
  ms = collect(ml)
  isempty(ms) && return

  substr = s"\g<sig> in \g<mod> at \g<loc>"
  msstr = map(ms) do m
    s = replace(string(m), methodloc_regex => substr)
    "
$s
"
  end |> join

  name = ms[1].name
  pushmarkdown!(datatip, "$name has **$(length(ms))** methods:
$(msstr)
")
end

function pushmarkdown!(datatip, markdown)
  (markdown == "" || markdown == "\n") && return
  push!(datatip, Dict(:type => :markdown, :value => markdown))
end

function pushsnippet!(datatip, snippet)
  push!(datatip, Dict(:type => :snippet, :value => snippet))
end


text:
using AnyAscii
using Test

@testset "AnyAscii.jl" begin
	check(s, expected) = @test anyascii(s) == expected

	check("", "")
	check("\x00\x01\t\n\x1f ~\x7f", "\x00\x01\t\n\x1f ~\x7f")
	check("sample", "sample")

	check(Char(0x0080), "")
	check(Char(0x00ff), "y")
	check(Char(0xe000), "")
	check(Char(0xfdff), "")
	check(Char(0x000e0020), " ")
	check(Char(0x000e007e), "~")
	check(Char(0x000f0000), "")
	check(Char(0x000f0001), "")
	check(Char(0x0010ffff), "")

	check("René François Lacôte", "Rene Francois Lacote")
	check("Blöße", "Blosse")
	check("Trần Hưng Đạo", "Tran Hung Dao")
	check("Nærøy", "Naeroy")
	check("Φειδιππίδης", "Feidippidis")
	check("Δημήτρης Φωτόπουλος", "Dimitris Fotopoylos")
	check("Борис Николаевич Ельцин", "Boris Nikolaevich El'tsin")
	check("Володимир Горбулін", "Volodimir Gorbulin")
	check("Търговище", "T'rgovishche")
	check("深圳", "ShenZhen")
	check("深水埗", "ShenShuiBu")
	check("화성시", "HwaSeongSi")
	check("華城市", "HuaChengShi")
	check("さいたま", "saitama")
	check("埼玉県", "QiYuXian")
	check("ደብረ ዘይት", "debre zeyt")
	check("ደቀምሓረ", "dek'emhare")
	check("دمنهور", "dmnhwr")
	check("Աբովյան", "Abovyan")
	check("სამტრედია", "samt'redia")
	check("אברהם הלוי פרנקל", "'vrhm hlvy frnkl")
	check("⠠⠎⠁⠽⠀⠭⠀⠁⠛", "+say x ag")
	check("ময়মনসিংহ", "mymnsimh")
	check("ထန်တလန်", "thntln")
	check("પોરબંદર", "porbmdr")
	check("महासमुंद", "mhasmumd")
	check("ಬೆಂಗಳೂರು", "bemgluru")
	check("សៀមរាប", "siemrab")
	check("ສະຫວັນນະເຂດ", "sahvannaekhd")
	check("കളമശ്ശേരി", "klmsseri")
	check("ଗଜପତି", "gjpti")
	check("ਜਲੰਧਰ", "jlmdhr")
	check("රත්නපුර", "rtnpur")
	check("கன்னியாகுமரி", "knniyakumri")
	check("శ్రీకాకుళం", "srikakulm")
	check("สงขลา", "sngkhla")
	check("😎 👑 🍎", ":sunglasses: :crown: :apple:")
	check("☆ ♯ ♰ ⚄ ⛌", "* # + 5 X")
	check("№ ℳ ⅋ ⅍", "No M & A/S")

	check("トヨタ", "toyota")
	check("ߞߐߣߊߞߙߌ߫", "konakri")
	check("𐬰𐬀𐬭𐬀𐬚𐬎𐬱𐬙𐬭𐬀", "zarathushtra")
	check("ⵜⵉⴼⵉⵏⴰⵖ", "tifinagh")
	check("𐍅𐌿𐌻𐍆𐌹𐌻𐌰", "wulfila")
	check("ދިވެހި", "dhivehi")
end