text:
using FeynmanKacParticleFilters, ExactWrightFisher

function multinomial_logpotential(obs_vector::AbstractArray{T, 1}) where T <: Real
    function pot(state::AbstractArray{T, 1})  where T <: Real
        return logpdf(Multinomial(sum(obs_vector), state), obs_vector)
    end
    return pot
end

function multinomial_logpotential(obs_vector::AbstractArray{T, 2}) where T <: Real
    if (size(obs_vector, 1) == 1) || size(obs_vector, 2) == 1
        return multinomial_logpotential(vec(obs_vector))
    else
        error("must write the multinomial potential function for multiple observations")
    end
end

function create_transition_kernels_WF(data, α_vec::AbstractArray{U, 1}) where {T <: Real, U <: Real}
    sα = sum(α_vec)
    K = length(α_vec)
    function create_Mt(Δt)
        return state -> Wright_Fisher_K_dim_transition_with_t005_approx(state, Δt, α_vec, sα)
    end

    prior = Dirichlet(repeat([1.], K))
    return FeynmanKacParticleFilters.create_transition_kernels(data, create_Mt, prior)
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
"""
    NumaAllocators

Extends `ArrayAllocators` to allocate memory on specific NUMA nodes.

# Examples

```julia
using NumaAllocators

Array{UInt8}(numa(0), 100)
Array{UInt8}(NumaAllocator(1), 100)
```
"""
module NumaAllocators

@static if VERSION >= v"1.3"
    using NUMA_jll
end

using ArrayAllocators: AbstractArrayAllocator

export NumaAllocator, numa
export current_numa_node, highest_numa_node

abstract type AbstractNumaAllocator{B} <: AbstractArrayAllocator{B} end

include("Windows.jl")
include("LibNUMA.jl")

@static if Sys.iswindows()
    import .Windows: WinNumaAllocator

    const NumaAllocator = WinNumaAllocator
elseif ( VERSION >= v"1.6" && NUMA_jll.is_available() ) || isdefined(NUMA_jll, :libnuma)
    import .LibNUMA: LibNumaAllocator

    const NumaAllocator = LibNumaAllocator
end

if @isdefined(NumaAllocator)
    numa(node) = NumaAllocator(node)
end

"""
    NumaAllocator(node)

Cross-platform NUMA allocator

# Example

```jldoctest
julia> using NumaAllocators

julia> Array{UInt8}(NumaAllocator(0), 32, 32);
```
"""
NumaAllocator

"""
    numa(node)

Create a `NumaAllocator` on NUMA node `node`. Short hand for [`NumaAllocator`](@ref) constructor.

# Example

```jldoctest
julia> using NumaAllocators

julia> Array{UInt8}(numa(0), 32, 32);
```
"""
numa

"""
    current_numa_node()::Int

Returns the current NUMA node as an Int
"""
function current_numa_node()::Int
    @static if Sys.iswindows()
        return Int(Windows.GetNumaProcessorNode())
    elseif Sys.islinux()
        return Int(LibNUMA.numa_node_of_cpu())
    else
        error("current_numa_node is not implemented for this platform")
    end
end

"""
    highest_numa_node()::Int

Returns the highest NUMA node as an Int
"""
function highest_numa_node()::Int
    @static if Sys.iswindows()
        return Int(Windows.GetNumaHighestNodeNumber())
    elseif Sys.islinux()
        return Int(LibNUMA.numa_max_node())
    else
        error("highest_numa_node is not implemented for this platform")
    end
end

end

ext: jl


text: @show binomial(5, 3)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Missings
using FWF
using DataStreams
using DataFrames

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end
if VERSION < v"0.7.0-DEV.2575"
    using Base.Dates
else
    using Dates
end

dir = joinpath(dirname(@__FILE__),"testfiles/")

include("FWF.jl")
include("Source.jl")
include("parsefields.jl")
include("io.jl")
include("scan.jl")


text:
module TestObjectiveFunctionize

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_solve_singlevariable_obj()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Objective.Functionize{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [1.0], MOI.FEASIBLE_POINT),
    )
    MOI.Test.test_objective_ObjectiveFunction_VariableIndex(
        bridged_mock,
        MOI.Test.Config(),
    )
    @test MOI.get(mock, MOI.ObjectiveFunctionType()) ==
          MOI.ScalarAffineFunction{Float64}
    @test MOI.get(bridged_mock, MOI.ObjectiveFunctionType()) ==
          MOI.VariableIndex
    @test MOI.get(mock, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    @test MOI.get(bridged_mock, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    vis = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
    func = vis[1]
    @test MOI.get(
        mock,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    ) ≈ convert(MOI.ScalarAffineFunction{Float64}, func)
    @test MOI.get(bridged_mock, MOI.ObjectiveFunction{MOI.VariableIndex}()) ==
          func
    MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(mock, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    @test MOI.get(bridged_mock, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    _test_delete_objective(bridged_mock, 1, tuple())
    return
end

function test_solve_result_index()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Objective.Functionize{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0]),
            MOI.FEASIBLE_POINT,
            (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [1.0],
        ),
    )
    return MOI.Test.test_solve_result_index(bridged_mock, MOI.Test.Config())
end

end  # module

TestObjectiveFunctionize.runtests()


text:
# Source term

function init_source(param::Param)
   (;GridSize, FullSize, nVar) = param

   source_GV = Array{Float64,4}(undef, GridSize..., nVar)
   div_G = Array{Float64,3}(undef, GridSize...)

   U = Array{Float64,4}(undef, FullSize..., 3)

   return source_GV, U, div_G
end

function calc_source!(param::Param, state_GV, source_GV, U, div_G)
   (;nI, nJ, nK, nG, iMin, iMax, jMin, jMax, kMin, kMax) = param

   # Compute ∇⋅B using central difference
   @views divergence!(param, state_GV[:,:,:,B_], div_G)

   source_GV[:,:,:,Rho_] .= 0.0

   @inbounds for k = 1:nK, j = 1:nJ, i = 1:nI
      source_GV[i,j,k,Mx_] = -state_GV[i+nG,j+nG,k+nG,Bx_]*div_G[i,j,k]
      source_GV[i,j,k,My_] = -state_GV[i+nG,j+nG,k+nG,By_]*div_G[i,j,k]
      source_GV[i,j,k,Mz_] = -state_GV[i+nG,j+nG,k+nG,Bz_]*div_G[i,j,k]
      source_GV[i,j,k,Bx_] =  state_GV[i+nG,j+nG,k+nG,Ux_]*div_G[i,j,k]
      source_GV[i,j,k,By_] =  state_GV[i+nG,j+nG,k+nG,Uy_]*div_G[i,j,k]
      source_GV[i,j,k,Bz_] =  state_GV[i+nG,j+nG,k+nG,Uz_]*div_G[i,j,k]
   end

   # Calculate divergence of U
   for k = 1:size(U,3), j = 1:size(U,2), i = 1:size(U,1)
      U[i,j,k,1] = state_GV[i,j,k,Mx_] / state_GV[i,j,k,Rho_]
      U[i,j,k,2] = state_GV[i,j,k,My_] / state_GV[i,j,k,Rho_]
      U[i,j,k,3] = state_GV[i,j,k,Mz_] / state_GV[i,j,k,Rho_]
   end
   divergence!(param, U, div_G)

   if !param.UseConservative
      @inbounds for k = 1:nK, j = 1:nJ, i = 1:nI
         source_GV[i,j,k,P_] = -(γ-1.0)*state_GV[i+nG,j+nG,k+nG,P_]*div_G[i,j,k]
      end
   else
      @inbounds for k = kMin:kMax, j = jMin:jMax, i = iMin:iMax
         ub = state_GV[i,j,k,Ux_]*state_GV[i,j,k,Bx_] +
            state_GV[i,j,k,Uy_]*state_GV[i,j,k,By_] +
            state_GV[i,j,k,Uz_]*state_GV[i,j,k,Bz_]
         source_GV[i-nG,j-nG,k-nG,E_] = -ub*div_G[i-nG,j-nG,k-nG]
      end
   end

   return source_GV
end


text:
#
# Basic two-body orbit simularion written in julia
#
using PyPlot

# Universal gravitation constant
const G = 6.67e-11

#
# Read speed, position and mass of first body
#
linia = readline()
pos = split(strip(linia), " ")
linia = readline()
spd = split(strip(linia), " ")
mass = strip(readline())

# Initialize and load position and velocity vectors
r1 = zeros(3,1)
v1 = zeros(3,1)
for i in 1:3
  r1[i] = float64(pos[i])
  v1[i] = float64(spd[i])
end
# Initialize mass
m1 = float64(mass)

#
# Read speed, position and mass of second body
#
linia = readline()
pos = split(strip(linia), " ")
linia = readline()
spd = split(strip(linia), " ")
mass = strip(readline())

# Initialize and load position and velocity vectors
r2 = zeros(3,1)
v2 = zeros(3,1)
for i in 1:3
  r2[i] = float64(pos[i])
  v2[i] = float64(spd[i])
end
# Initialize mass
m2 = float64(mass)

# Time increment (in seconds) and number of steps
# will be a day.
linia = readline()
times = split(strip(linia), " ")
nintv = int32(times[2])
dt = float64(times[1])

# Newton's formula for gravitational force
function grav(r1,r2,m1,m2)
  r = norm(r1-r2)
  f1 = G * m1 * m2 * (r2 - r1) / r^3
  f2 = G * m1 * m2 * (r1 - r2) / r^3
  return f1,f2
end

# Newton's formula for acceleration
function accel(f,m)
  return f/m
end

# Apply a change (accel or velocity) to a vector for an interval (dt)
function applydt(vector, change, dt)
  return vector + change * dt
end

# Matrices to contain the (x,y) coordinates computed for each iteration
# for bodies 1 and 2
C1 = zeros(nintv ,2)
C2 = zeros(nintv ,2)

# Compute the orbit position for each step
for i in 1:nintv
  C1[i,1] = r1[1]
  C1[i,2] = r1[2]
  C2[i,1] = r2[1]
  C2[i,2] = r2[2]

  f1,f2 = grav(r1,r2,m1,m2)
  a1 = accel(f1,m1)
  a2 = accel(f2,m2)
  v1 = applydt(v1,a1,dt)
  v2 = applydt(v2,a2,dt)
  r1 = applydt(r1,v1,dt)
  r2 = applydt(r2,v2,dt)
  println(i," r1=",r1,",r2=",r2)
end

# Plot both orbits
scatter(C1[1:nintv,1],C1[1:nintv,2],marker=".")
scatter(C2[1:nintv,1],C2[1:nintv,2],marker=",")

max_stars_repo_name: jguillaumes/orbit-examples

max_issues_repo_name: jguillaumes/orbit-examples

max_forks_repo_name: jguillaumes/orbit-examples


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct RouteTablePropertiesFormat <: SwaggerModel
    routes::Any # spec type: Union{ Nothing, Vector{Route} } # spec name: routes
    subnets::Any # spec type: Union{ Nothing, Vector{Subnet} } # spec name: subnets
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function RouteTablePropertiesFormat(;routes=nothing, subnets=nothing, provisioningState=nothing)
        o = new()
        validate_property(RouteTablePropertiesFormat, Symbol("routes"), routes)
        setfield!(o, Symbol("routes"), routes)
        validate_property(RouteTablePropertiesFormat, Symbol("subnets"), subnets)
        setfield!(o, Symbol("subnets"), subnets)
        validate_property(RouteTablePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type RouteTablePropertiesFormat

const _property_map_RouteTablePropertiesFormat = Dict{Symbol,Symbol}(Symbol("routes")=>Symbol("routes"), Symbol("subnets")=>Symbol("subnets"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_RouteTablePropertiesFormat = Dict{Symbol,String}(Symbol("routes")=>"Vector{Route}", Symbol("subnets")=>"Vector{Subnet}", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ RouteTablePropertiesFormat }) = collect(keys(_property_map_RouteTablePropertiesFormat))
Swagger.property_type(::Type{ RouteTablePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_RouteTablePropertiesFormat[name]))}
Swagger.field_name(::Type{ RouteTablePropertiesFormat }, property_name::Symbol) =  _property_map_RouteTablePropertiesFormat[property_name]

function check_required(o::RouteTablePropertiesFormat)
    true
end

function validate_property(::Type{ RouteTablePropertiesFormat }, name::Symbol, val)
end


text:
using SwitchTimeOpt
using Base.Test
using Ipopt

include("linsys.jl")
linsystest()

include("nonlinsys.jl")
nonlinsystest()


text:
export coordinates, absolute_coordinates, absolute_norm, absolute_tr,
       absolute_minpoly

################################################################################
#
#
#
################################################################################

isunit(a::NumFieldElem) = !iszero(a)

canonical_unit(a::NumFieldElem) = a

################################################################################
#
#  Base case for dot products
#
################################################################################

dot(x::fmpz, y::NumFieldElem) = x * y

dot(x::Integer, y::NumFieldElem) = x * y

dot(x::NumFieldElem, y::Integer) = x * y

function dot(a::Vector{<: NumFieldElem}, b::Vector{fmpz})
  d = zero(parent(a[1]))
  t = zero(d)
  for i=1:length(a)
    mul!(t, a[i], b[i])
    add!(d, d, t)
  end
  return d
end

function dot(a::Vector{NfAbsNSElem}, b::Vector{fmpz})
  Qxy = parent(a[1].data)
  d = zero(Qxy)
  t = zero(Qxy)
  for i = 1:length(a)
    if !iszero(b[i])
      mul!(t, a[i].data, b[i])
      add!(d, d, t)
    end
  end
  return parent(a[1])(d)
end

################################################################################
#
#  Parent
#
################################################################################

# Covered by the general ring interface

################################################################################
#
#  Integrality
#
################################################################################

isintegral(a::fmpq) = isone(denominator(a))

@doc doc"""
    isintegral(a::NumFieldElem) -> Bool

Returns whether $a$ is integral, that is, whether the minimal polynomial of $a$
has integral coefficients.
"""
function isintegral(a::NumFieldElem)
  K = parent(a)
  if ismaximal_order_known(K)
    OK = maximal_order(K)
    return a in OK
  end
  f = minpoly(a)
  for i in 0:(degree(f) - 1)
    if !isintegral(coeff(f, i))
      return false
    end
  end
  return true
end

################################################################################
#
#  Random elements from arrays of number field elements
#
################################################################################

## rand(::Vector{NumFieldElem}, ::UnitRange)

@doc doc"""
    rand([rng::AbstractRNG], b::Vector{NumFieldElem}, r::UnitRange) -> NumFieldElem
    rand([rng::AbstractRNG], make(F::NumField, b::Vector{NumFieldElem}, r::UnitRange)) -> NumFieldElem

A random linear combination of elements in `b`, with parent `F` and coefficients in `r`.
"""
rand(b::Vector{<: NumFieldElem}, r::UnitRange) = rand(Random.GLOBAL_RNG, b, r)

function rand(rng::AbstractRNG, b::Vector{<: NumFieldElem}, r::UnitRange)
  length(b) == 0 && error("Array must not be empty")
  return rand(rng, make(parent(b[1]), b, r))
end

function rand(rng::AbstractRNG,
              sp::SamplerTrivial{<:Make3{<:NumFieldElem{T},<:NumField{T},
                                         <:Vector{<:NumFieldElem{T}},<:UnitRange}}) where {T}
  return rand!(rng, zero(sp[][1]), sp)
end

## rand(::Vector{<: NumFieldElem}, ::UnitRange, ::Int)

@doc doc"""
    rand([rng::AbstractRNG], b::Vector{NumFieldElem}, r::UnitRange, terms::Int) -> NumFieldElem
    rand([rng::AbstractRNG],
         make(F::NumField, b::Vector{NumFieldElem}, r::UnitRange, terms::Int)) -> NumFieldElem

A random linear combination (with repetitions) of `terms` elements of `b`
with parent `F` and coefficients in `r`.
"""
function rand(b::Vector{<: NumFieldElem}, r::UnitRange, terms::Int)
  return rand(Random.GLOBAL_RNG, b, r, terms)
end

function rand(rng::AbstractRNG, b::Vector{<: NumFieldElem}, r::UnitRange, terms::Int)
  length(b) == 0 && error("Array must not be empty")
  return rand(rng, make(parent(b[1]), b, r, terms))
end

function rand(rng::AbstractRNG,
              sp::SamplerTrivial{<:Make4{<:NumFieldElem{T},<:NumField{T},
                                         <:Vector{<:NumFieldElem{T}},<:UnitRange,Int}}) where {T}
  return rand!(rng, zero(sp[][1]), sp)
end

## rand!(::NumFieldElem, ::Vector{NumFieldElem}, ::UnitRange, terms::Int)

@doc doc"""
    rand!([rng::AbstractRNG], c::NumFieldElem, b::Vector{NumFieldElem},
          r::UnitRange, terms::Int) -> NumFieldElem
    rand!([rng::AbstractRNG], c::NumFieldElem,
          make(F::NumField, b::Vector{NumFieldElem}, r::UnitRange, terms::Int)) -> NumFieldElem

Sets `c` to a random linear combination (with repetitions) of \code{terms}
elements of `b` with coefficients in `r`. The element `c` is returned.
In the second form, `F` must be the parent of `c`.
"""
function rand!(c::T, b::Vector{T}, r::UnitRange, terms::Int) where {T <: NumFieldElem}
  rand!(Random.GLOBAL_RNG, c, b, r, terms)
end

function rand!(rng::AbstractRNG, c::T, b::Vector{T}, r::UnitRange, terms::Int) where {T <: NumFieldElem}
  return rand!(rng, c, make(parent(c), b, r, terms))
end

function rand!(rng::AbstractRNG, c::NumFieldElem{T},
               sp::SamplerTrivial{<:Make4{<:NumFieldElem{T},<:NumField{T},
                                          <:Vector{<:NumFieldElem{T}},<:UnitRange,Int}}) where {T}
  F, b, r, terms = sp[][1:end]

  length(b) == 0 && error("Array must not be empty")
  0 < terms <= length(b) || error("Number of terms should be at least 1 and cannot exceed the length of the array")
  F == parent(c) || error("incompatible parent and element")

  t = zero(F)

  mul!(c, rand(rng, b), rand(rng, r))
  for i = 2:terms
    mul!(t, rand(rng, b), rand(rng, r))
    add!(c, t, c)
  end

  return c
end

## rand!(::NumFieldElem, ::Vector{NumFieldElem}, ::UnitRange)

@doc doc"""
    rand!(c::NumFieldElem, b::Vector{NumFieldElem}, r::UnitRange) -> NumFieldElem
    rand!(c::NumFieldElem, make(F::NumField, b::Vector{NumFieldElem}, r::UnitRange)) -> NumFieldElem

Sets `c` to a random linear combination of elements in `b` with coefficients
in `r`. The element `c` is returned.
In the second form, `F` must be the parent of `c`.
"""
function rand!(c::T, b::Vector{T}, r::UnitRange) where {T <: NumFieldElem}
  return rand!(GLOBAL_RNG, c, b, r)
end

function rand!(rng::AbstractRNG, c::T, b::Vector{T}, r::UnitRange) where {T <: NumFieldElem}
  return rand!(rng, c, make(parent(c), b, r))
end

RandomExtensions.maketype(F::NumField{T}, ::Vector{<:NumFieldElem{T}}, ::UnitRange, ::Int...) where {T} = elem_type(F)

function rand!(rng::AbstractRNG, c::NumFieldElem{T},
               sp::SamplerTrivial{<:Make3{<:NumFieldElem{T},<:NumField{T},
                                          <:Vector{<:NumFieldElem{T}},<:UnitRange}}) where {T}
  F, b, r = sp[][1:end]

  length(b) == 0 && error("Array must not be empty")
  F == parent(c) || error("incompatible parent and element")

  mul!(c, b[1], rand(rng, r))
  t = zero(F)

  for i = 2:length(b)
    t = mul!(t, b[i], rand(rng, r))
    c = add!(c, t, c)
  end

  return c
end

################################################################################
#
#  Basis matrix
#
################################################################################

@doc doc"""
    basis_matrix(v::Vector{NumFieldElem}) -> Mat

Given a vector $v$ of $n$ elements of a number field $K$ of degree $d$, this
function returns an $n x d$ matrix with entries in the base field of $K$, where
row $i$ contains the coefficients of $v[i]$ with respect of the canonical
basis of $K$.
"""
basis_matrix(v::Vector{<: NumFieldElem})

################################################################################
#
#  Characteristic polynomial
#
################################################################################

@doc doc"""
    charpoly(a::NumFieldElem) -> PolyElem

Given a number field element $a$ of a number field $K$, this function returns
the characteristic polynomial of $a$ over the base field of $K$.
"""
charpoly(::NumFieldElem)

@doc doc"""
    absolute_charpoly(a::NumFieldElem) -> PolyElem

Given a number field element $a$ of a number field $K$, this function returns
the characteristic polynomial of $a$ over the rationals $\mathbf{Q}$.
"""
absolute_charpoly(::NumFieldElem)

@doc doc"""
    minpoly(a::NumFieldElem) -> PolyElem

Given a number field element $a$ of a number field $K$, this function returns
the minimal polynomial of $a$ over the base field of $K$.
"""
minpoly(::NumFieldElem)

@doc doc"""
    absolute_minpoly(a::NumFieldElem) -> PolyElem

Given a number field element $a$ of a number field $K$, this function returns
the minimal polynomial of $a$ over the rationals $\mathbf{Q}$.
"""
absolute_minpoly(::NumFieldElem)

################################################################################
#
#  Powering with fmpz
#
################################################################################

function ^(x::NumFieldElem, y::fmpz)
  if fits(Int, y)
    return x^Int(y)
  end

  return _power(x, y)
end

# We test once if it fits, otherwise we would have to check for every ^-call
function _power(x::NumFieldElem, y::fmpz)
  res = parent(x)()
  if y < 0
    res = _power(inv(x), -y)
  elseif y == 0
    res = parent(x)(1)
  elseif y == 1
    res = deepcopy(x)
  elseif mod(y, 2) == 0
    z = _power(x, div(y, 2))
    res = z*z
  else
    res = _power(x, y - 1) * x
  end
  return res
end

################################################################################
#
#  SimpleNumFieldElem
#
################################################################################

@doc doc"""
    coeff(a::SimpleNumFieldElem, i::Int) -> FieldElem

Given a number field element `a` of a simple number field extension `L/K`, this
function returns the `i`-th coefficient of `a`, when expanded in the canonical
power basis of `L`. The result is an element of `K`.
"""
coeff(::SimpleNumFieldElem, ::Int)

# copy does not do anything (so far), this is only for compatibility with coefficients(::AbsAlgAssElem)

@doc doc"""
    coefficients(a::SimpleNumFieldElem, i::Int) -> Vector{FieldElem}

Given a number field element `a` of a simple number field extension `L/K`, this
function returns the coefficients of `a`, when expanded in the canonical
power basis of `L`.
"""
function coefficients(a::SimpleNumFieldElem; copy::Bool = false)
  return [ coeff(a, i) for i = 0:degree(parent(a)) - 1 ]
end

@doc doc"""
    (L::SimpleNumField)(c::Vector) -> SimpleNumFieldElem

Given a simple number field extension `L/K` and a vector of elements that are
coercible to `K`, this functions returns the element of `L` with coefficients
`c`.
"""
(K::SimpleNumField)(c::Vector)

function (K::AnticNumberField)(c::Vector{fmpq})
  @assert length(c) == degree(K)
  return K(parent(K.pol)(c))
end

function (L::NfRel{T})(a::Vector{T}) where T
  z = NfRelElem{T}(parent(L.pol)(a))
  z.parent = L
  return z
end

for F in [AnticNumberField, NfRel]
  @eval begin
    function (L::$F)(a::Vector)
      length(a) != degree(L) && throw(error("Vector must have length ($(length(a))) equal to the degree ($(degree(L)))"))
      K = base_field(L)
      aprom = Vector{elem_type(K)}(undef, degree(L))
      for i in 1:degree(L)
        aprom[i] = K(a[i])
      end
      return L(aprom)
    end
  end
end

################################################################################
#
#  Norm and trace
#
################################################################################

@doc doc"""
    tr(a::NumFieldElem) -> NumFieldElem

Returns the trace of an element $a$ of a number field extension $L/K$. This
will be an element of $K$.
"""
tr(::NumFieldElem)

@doc doc"""
    norm(a::NumFieldElem) -> NumFieldElem

Returns the norm of an element $a$ of a number field extension $L/K$. This
will be an element of $K$.
"""
norm(::NumFieldElem)

function _elem_tr_to(a, k::T) where {T}
  if parent(a) isa T
    @assert parent(a) == k
    return a
  else
    return _elem_tr_to(tr(a), k)
  end
end

function _elem_norm_to(a, k::T) where {T}
  if parent(a) isa T
    @assert parent(a) == k
    return a
  else
    return _elem_norm_to(norm(a), k)
  end
end

@doc doc"""
    tr(a::NumFieldElem, k::NumField) -> NumFieldElem

Returns the trace of an element $a$ of a number field $L$ with respect to
a subfield $k$ of $L$. This will be an element of $k$.
"""
function tr(a::NumFieldElem, k::NumField)
  return _elem_tr_to(a, k)
end

tr(a::NumFieldElem, ::FlintRationalField) = _elem_tr_to(a, FlintQQ)

@doc doc"""
    norm(a::NumFieldElem, k::NumField) -> NumFieldElem

Returns the norm of an element $a$ of a number field $L$ with respect to
a subfield $k$ of $L$. This will be an element of $k$.
"""
function norm(a::NumFieldElem, k::NumField)
  _elem_norm_to(a, k)
end

norm(a::NumFieldElem, ::FlintRationalField) = _elem_norm_to(a, FlintQQ)

@doc doc"""
    absolute_tr(a::NumFieldElem) -> fmpq

Given a number field element $a$, returns the absolute trace of $a$.
"""
function absolute_tr(a::T) where T <: NumFieldElem
  return absolute_tr(tr(a))
end

absolute_tr(a::nf_elem) = tr(a)

absolute_tr(a::NfAbsNSElem) = tr(a)

absolute_tr(x::fmpq) = x

@doc doc"""
    absolute_norm(a::NumFieldElem) -> fmpq

Given a number field element $a$, returns the absolute norm of $a$.
"""
function absolute_norm(a::T) where T <: NumFieldElem
  return absolute_norm(norm(a))
end

absolute_norm(a::T) where T <: Union{nf_elem, NfAbsNSElem} = norm(a)

absolute_norm(a::fmpq) = a

################################################################################
#
#  Polynomials
#
################################################################################

@doc doc"""
    norm(f::PolyElem{<:NumFieldElem}) -> PolyElem

Returns the norm of $f$, that is, the product of all conjugates of $f$ taken
coefficientwise.
"""
function norm(f::PolyElem{<: NumFieldElem})
  K = base_ring(f)
  P = polynomial_to_power_sums(f, degree(f)*degree(K))
  PQ = elem_type(base_field(K))[tr(x) for x in P]
  return power_sums_to_polynomial(PQ)
end

function norm(f::PolyElem{<:NumFieldElem}, k::NumField)
  K = base_ring(f)
  P = polynomial_to_power_sums(f, degree(f)*degree(K))
  PQ = elem_type(base_field(K))[tr(x, k) for x in P]
  return power_sums_to_polynomial(PQ)
end

norm(a::fmpq_poly) = a

function absolute_norm(f::PolyElem{nf_elem})
  return norm(f)
end

function absolute_norm(f::PolyElem{<: NumFieldElem})
  return absolute_norm(norm(f))
end

function isirreducible(f::PolyElem{<: NumFieldElem})
  # TODO (easy): We can do better then this. First do a squarefree factorization
  lf = factor(f)
  return sum(values(lf.fac)) == 1
end

function AbstractAlgebra.factor(f::PolyElem{<: NumFieldElem})
  K = base_ring(f)
  Ka, mKa = absolute_simple_field(K)

  fKa = map_coefficients(inv(mKa), f)
  lf = factor(fKa)
  res = Fac(map_coefficients(mKa, lf.unit, parent = parent(f)), Dict(map_coefficients(mKa, k, parent = parent(f)) => v for (k,v) = lf.fac))

  return res
end

function roots(f::PolyElem{<: NumFieldElem})
  lf = factor(f)
  @assert degree(unit(lf)) == 0
  scale = inv(coeff(unit(lf), 0))
  return elem_type(base_ring(f))[-constant_coefficient(x)*scale for x = keys(lf.fac) if degree(x) == 1]
end

################################################################################
#
#
#
################################################################################

@doc doc"""
    representation_matrix(a::NumFieldElem) -> MatElem

Returns the representation matrix of $a$, that is, the matrix representing
multiplication with $a$ with respect to the canonical basis of the parent of $a$.
"""
representation_matrix(a::NumFieldElem)

################################################################################
#
#
#
################################################################################

@doc doc"""
    gen(L::SimpleNumField) -> NumFieldElem

Given a simple number field $L = K[x]/(x)$ over $K$, this functions returns the
class of $x$, which is the canonical primitive element of $L$ over $K$.
"""
gen(::SimpleNumField)

@doc doc"""
    gens(L::NonSimpleNumField) -> Vector{NumFieldElem}

Given a non-simple number field $L = K[x_1,\dotsc,x_n]/(f_1,\dotsc,f_n)$ over
$K$, this functions returns the list $\bar x_1,\dotsc,\bar x_n$.
"""
gens(::NonSimpleNumField)

################################################################################
#
#  Coordinates
#
################################################################################

@doc doc"""
    coordinates(x::NumFieldElem{T}) -> Vector{T}

Given an element $x$ in a number field $K$, this function returns the coordinates of $x$
with respect to the basis of $K$ (the output of the 'basis' function).
"""
coordinates(::NumFieldElem)

function coordinates(a::nf_elem)
  K = parent(a)
  v = Vector{fmpq}(undef, degree(K))
  for i = 1:length(v)
    v[i] = coeff(a, i-1)
  end
  return v
end

function coordinates(a::NfAbsNSElem)
  adata = data(a)
  K = parent(a)
  v = Vector{fmpq}(undef, degree(K))
  for i = 1:length(v)
    v[i] = fmpq(0)
  end
  for j in 1:length(adata)
    exps = exponent_vector(adata, j)
    k = monomial_to_index(K, exps)
    v[k] = coeff(adata, j)
  end
  return v
end

function coordinates(a::NfRelElem{T}) where T
  K = parent(a)
  k = base_field(K)
  v = Vector{T}(undef, degree(K))
  for i = 1:length(v)
    v[i] = coeff(data(a), i-1)
  end
  return v
end

function coordinates(a::NfRelNSElem{T}) where T
  K = parent(a)
  k = base_field(K)
  v = Vector{T}(undef, degree(K))
  for j = 1:length(v)
    v[j] = zero(k)
  end
  for j = 1:length(a.data)
    v[monomial_to_index(j, a)] = a.data.coeffs[j]
  end
  return v
end

################################################################################
#
#  Absolute coordinates
#
################################################################################

@doc doc"""
    absolute_coordinates(x::NumFieldElem{T}) -> Vector{T}

Given an element $x$ in a number field $K$, this function returns the coordinates of $x$
with respect to the basis of $K$ over the rationals (the output of the 'absolute_basis' function).
"""
absolute_coordinates(::NumFieldElem)

function absolute_coordinates(a::NumFieldElem{fmpq})
  return coordinates(a)
end

function absolute_coordinates(a::T) where T <: Union{NfRelElem, NfRelNSElem}
  K = parent(a)
  v = Vector{fmpq}(undef, absolute_degree(K))
  va = coordinates(a)
  ind = 1
  for i = 1:length(va)
    vi = absolute_coordinates(va[i])
    for j = 1:length(vi)
      v[ind] = vi[j]
      ind += 1
    end
  end
  return v
end

################################################################################
#
#  Denominator
#
################################################################################

function denominator!(z::fmpz, a::nf_elem)
   ccall((:nf_elem_get_den, libantic), Nothing,
         (Ref{fmpz}, Ref{nf_elem}, Ref{AnticNumberField}),
         z, a, a.parent)
   return z
end

################################################################################
#
#  Valuation
#
################################################################################

@doc Markdown.doc"""
    valuation(a::NumFieldElem, p::NfOrdIdl) -> fmpz

Computes the $\mathfrak p$-adic valuation of $a$, that is, the largest $i$
such that $a$ is contained in $\mathfrak p^i$.
"""
function valuation(::NumFieldElem, p) end

################################################################################
#
#   Support
#
################################################################################

function support(a::NumFieldElem{fmpq}, R::NfAbsOrd = maximal_order(parent(a)))
  @assert nf(R) == parent(a)
  return support(a * R)
end

function support(a::NumFieldElem, R::NfRelOrd = maximal_order(parent(a)))
  @assert nf(R) == parent(a)
  return support(a * R)
end

################################################################################
#
#   Absolute minpoly
#
################################################################################

function minpoly(a::T, ::FlintRationalField) where T <: Union{NfRelNSElem, NfRelElem}
  f = minpoly(a)
  n = absolute_norm(f)
  g = gcd(n, derivative(n))
  if isone(g)
    return n
  end
  n = divexact(n, g)
  return n
end

absolute_minpoly(a::nf_elem) = minpoly(a)
absolute_minpoly(a::NfAbsNS) = minpoly(a)

absolute_minpoly(a::T) where T <: Union{NfRelNSElem, NfRelElem} = minpoly(a, FlintQQ)


text:
basis(::Basis) = error("Not implemented!")
nodes(::Basis) = error("Not implemented!")

nbasis(::Basis) = error("Not implemented!")
nnodes(::Basis) = error("Not implemented!")

order(::Basis) = error("Not implemented!")
degree(::Basis) = error("Not implemented!")


text:
function reduceband(A)
    l,u = bandwidths(A)
    H = _BandedMatrix(A.data, ∞, l+u-1, 1)
    Q1,L1 = ql(H)
    D = Q1[1:l+u+1,1:1]'A[1:l+u+1,1:u-1]
    D, Q1, L1
end

function householderiterate(Q::AbstractMatrix{T}, n) where T
    ret = Matrix{T}(I,n,n)
    for j = n-1:-1:1
        ret[j:j+1,:] = Q*ret[j:j+1,:]
    end
    ret
end

function householderiterate(Q1::AbstractMatrix{T}, Q2::AbstractMatrix{T}, n) where T
    ret = Matrix{T}(I,n,n)
    for j = n-1:-1:1
        ret[j:j+1,:] = (iseven(j) ? Q1 : Q2)*ret[j:j+1,:]
    end
    ret
end

function householderiterate(Q3::AbstractMatrix{T}, Q1::AbstractMatrix{T}, Q2::AbstractMatrix{T}, n) where T
    ret = Matrix{T}(I,n,n)
    for j = n-1:-1:1
        if mod(j,3) == 0
            ret[j:j+1,:] = Q1*ret[j:j+1,:]
        elseif mod(j,3) == 1
            ret[j:j+1,:] = Q2*ret[j:j+1,:]
        else
            ret[j:j+1,:] = Q3*ret[j:j+1,:]
        end
    end
    ret
end

@testset "Householder combine" begin
    @testset "sign = -1" begin
        σ = -0.9998783597231515 + 0.01559697910943707im
        v = -0.18889614060709453 + 7.427756729736341e-19im
        τ = 1.9310951421717215 - 7.593434620701001e-18im
        Z,A,B = 2,2.1+0.01im,0.5
        a = [Z,A,B]
        n = 1_000; T = Tridiagonal(Fill(ComplexF64(Z),∞), Fill(ComplexF64(A),∞), Fill(ComplexF64(B),∞)); Qn,Ln = ql(BandedMatrix(T)[1:n,1:n]);
        @test T isa TriToeplitz
        @test InfToeplitz(T) isa InfToeplitz

        de = d,e =tail_de([Z,A,B])
        X =  [Z A B; 0 d e]
        F = ql_X!(copy(X))
        H = I - τ *[v,1]*[v,1]'
        @test H ≈ QLPackedQ([NaN v; NaN NaN], [0,τ])
        Qt = [σ 0; 0 1] * H
        @test Qt ≈ F.Q'
        Qt2 = [1 0; 0 -1] * Qt
        @test (householderiterate(Qt2, 11)')[1:10,1:10] ≈ Qn[1:10,1:10]
        @test Qt2*X  ≈ [1 0; 0 -1] * ql(X).L

        t, ω = Qn.τ[2], Qn.factors[1,2]
        Q̃ = QLPackedQ([NaN ω; NaN NaN], [0, t])
        @test I - t*[ω,1]*[ω,1]' ≈ Q̃
        @test I - t'*[ω,1]*[ω,1]' ≈ Q̃'
        @test Q̃'[1,2] ≈ Qt2[1,2]
        @test Q̃'[2,1] ≈ Qt2[2,1]
        # (1-τ)/z == Qt[2,2]
        z = (1-τ)/(1-t')
        s = sqrt(z)
        @test abs(z) ≈ 1
        @test [-s 0; 0 1/s] * Qt2 * [s 0 ; 0 -1/s]  ≈ Q̃' ≈ I - t'*[ω,1]*[ω,1]'

        # [-σ*(1-τ*abs2(v)) -σ*τ*v; τ*conj(v) 1-τ]
        # [-z*σ*(1-τ*abs2(v)) -σ*τ*v; τ*conj(v)*(1-τ)/z]

        @test (1-τ) ≈ (1-t')*z
        @test -z*σ*(1-τ*abs2(v)) ≈ 1-t'*abs2(ω)
        @test -σ*τ*v ≈ -t'*ω
        @test τ*v' ≈ -t'*ω'

        @test (1-τ) ≈ (1-t')*z

        @test -z*σ*(1-τ*abs2(v)) ≈ 1-σ*τ*v*ω'
        @test -z*σ*(1-τ*abs2(v))*t' ≈ t'-σ*τ*v*ω'*t'
        @test -z*σ*(1-τ*abs2(v))*t' ≈ t'+τ*conj(v)*σ*τ*v
        @test (-(1-t')*z*σ*(1-τ*abs2(v))-(1-t'))*t' ≈ (1-t')*τ*conj(v)*σ*τ*v
        @test -t'*(1-τ)*σ*(1-τ*abs2(v))-(1-t')*t' ≈ (1-t')*abs2(v)*σ*τ^2
        @test -t'*((1-τ)*σ*(1-τ*abs2(v))+1-abs2(v)*σ*τ^2) + t'*t' ≈ abs2(v)*σ*τ^2

        β = -(σ*(1-τ*abs2(v))-σ*τ+1); γ = -abs2(v)*σ*τ^2;
        @test t' ≈ (-β + sqrt(β^2 - 4γ))/2
        @test ω ≈ σ*τ*v/t'
        @test abs2((1-τ)/(1-t')) ≈ 1

        # Check other sign choice fails
        β2 =σ-τ*σ-σ*τ*abs2(v)-1; γ2 = abs2(v)*σ*τ^2;
        t2t = (-β2 + sqrt(β2^2 - 4γ2))/2
        @test !(abs2((1-τ)/(1-t2t)) ≈ 1)

        @test all(householderparams(ql_X!(copy(X))) .≈ (σ,τ,v))
        @test all(combine_two_Q(σ,τ,v) .≈ (-1,t,ω))

        X = [transpose(a); 0 d e]
        ql_X!(X)
        s,t,ω = combine_two_Q(σ,τ,v) # combined two Qs into one, these are the parameteris
        Q∞11 = 1 - ω*t*conj(ω)  # Q[1,1] without the callowing correction
        Q̃n = QLPackedQ(Qn.factors, [0; Qn.τ[2:end]])
        @test Q∞11 ≈ Q̃n[1,1]
        τ1 = 1 - (a[end-1] -t*ω * X[2,end-1])/(Q∞11 * de[end]) # Choose τ[1] so that (Q*L)[1,1] = A
        @test τ1 ≈ Qn.τ[1]

        M = BandedMatrix(0 => Fill(A,∞), -1 => Fill(Z,∞), 1 => Fill(B,∞))
        Q∞, L∞ = ql(M)
        @test Q∞[1:10,1:10] ≈ Qn[1:10,1:10]
        @test L∞[1:10,1:10] ≈ Ln[1:10,1:10]
        @test Q∞[1:10,1:11]*L∞[1:11,1:10] ≈ M[1:10,1:10]
    end

    @testset "sign = +1" begin
        n = 1000
        T = BandedMatrix(-2 => Fill(1,∞), 0 => Fill(-0.5-0.1im, ∞), 1 => Fill(0.25, ∞));  Qn,Ln = ql(T[1:n,1:n]);
        a = reverse(T.data.args[1])
        de = tail_de(a)
        X =  [transpose(a); 0 transpose(de)]
        F = ql_X!(copy(X))
        @test F.factors[1,1:3] ≈ de
        @test F.factors[2,:] ≈ Ln[4,1:4]
        (σ,τ,v) = householderparams(F)
        H = I - τ *[v,1]*[v,1]'
        @test H ≈ QLPackedQ([NaN v; NaN NaN], [0,τ])
        Qt = [σ 0; 0 1] * H
        @test Qt ≈ F.Q'
        @test ((diagm(0 => [-1; Ones(10)]) * householderiterate(Qt, 11))')[1:10,1:10] ≈ Qn[1:10,1:10]

        t, ω = Qn.τ[2], Qn.factors[1,2]
        Q̃ = QLPackedQ([NaN ω; NaN NaN], [0, t])

        z = (1-τ)/(1-t')
        @test z ≈ (1-τ)/conj(Q̃[2,2])
        s = sqrt(z)
        @test [s 0; 0 1/s] * Qt * [s 0 ; 0 1/s]  ≈ Q̃' ≈ I - t'*[ω,1]*[ω,1]'
        @test z*σ*(1-τ*abs2(v)) ≈ 1-t'*abs2(ω)
        @test τ*v' ≈ t'*ω'

        # this derivation reduces to quadratic equation for t' in terms of knowns v, σ, and τ
        @test z*σ*(1-τ*abs2(v)) ≈ 1-σ*τ*v*ω'               ≈  Q̃[1,1]'
        @test z*σ*(1-τ*abs2(v))*t' ≈ t'-σ*τ*v*ω'*t'         ≈  Q̃[1,1]' * t'
        @test z*σ*(1-τ*abs2(v))*t' ≈ t'-σ*τ*v*τ*v'          ≈ Q̃[1,1]' * t'
        @test (1-τ)*σ*(1-τ*abs2(v))*t' ≈ (1-t')*(t'-σ*τ*v*τ*v')          ≈ (1-t')*Q̃[1,1]' * t'
        @test (t')^2 + ((1-τ)*σ*(1-τ*abs2(v))-1-σ*τ*v*τ*v')*t' + σ*τ^2*abs2(v) ≈  0  atol=1E-14  # simplify to stabndard form

        β =σ-τ*σ-σ*τ*abs2(v)-1; γ = abs2(v)*σ*τ^2;
        @test t' ≈ (-β + sqrt(β^2 - 4γ))/2
        @test ω ≈ σ*τ*v/t'
        @test abs2((1-τ)/(1-t')) ≈ 1

        β1 = -(σ*(1-τ*abs2(v))-σ*τ+1); γ1 = -abs2(v)*σ*τ^2;
        t1t = (-β1 + sqrt(β1^2 - 4γ1))/2
        @test !(abs2((1-τ)/(1-t1t)) ≈ 1)

        @test all(householderparams(ql_X!(copy(X))) .≈ (σ,τ,v))
        @test all(combine_two_Q(σ,τ,v) .≈ (1,t,ω))
        Q∞11 = 1 - ω*t*conj(ω)  # Q[1,1] without the callowing correction
        Q̃n = QLPackedQ(Qn.factors, [0; Qn.τ[2:end]])
        @test Q∞11 ≈ Q̃n[1,1]
        @test  -t*ω ≈ Q̃n[1,2]
        @test de[end] ≈ -Ln[1,1]
        @test (Q̃n*diagm(0 => [1-Qn.τ[1]; Ones(n-1)]))[1:10,1:10] ≈ Qn[1:10,1:10]
        @test Qn[1,1] * Ln[1,1] + Qn[1,2] * Ln[2,1] ≈ T[1,1]
        @test Q∞11 * (Qn.τ[1]-1) * de[end] - t*ω * F.factors[2,end-1] ≈ T[1,1]

        τ1 = (a[end-1] +t*ω * F.factors[2,end-1])/(Q∞11 * de[end])+1 # Choose τ[1] so that (Q*L)[1,1] = A
        @test τ1 ≈ Qn.τ[1]

        Q, L = ql(T)
        @test L[1,1] ≈ Ln[1,1] ≈ -de[end]
        @test L[1:10,1:10] ≈  Ln[1:10,1:10]
        @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]
    end
    @testset "periodic 1" begin
        A = BandedMatrix(-2 => Fill(1/4,∞), 1 => Fill(1,∞))-im*I
        n = 1000; Qn,Ln = ql(A[1:n,1:n])
        a = reverse(A.data.args[1])
        de = tail_de(a)
        X =  [transpose(a); 0 transpose(de)]
        F = ql_X!(copy(X))
        @test F.factors[1,1:3] ≈ de
        @test F.factors[2,:] ≈ -Ln[4,1:4]
        @test F.factors[2,:] ≈ Ln[5,2:5]

        # @test F.factors[2,:] ≈ Ln[4,1:4]
        # @test F.factors[2,:] ≈ -Ln[5,2:5]

        (σ,τ,v) = householderparams(F)
        H = I - τ *[v,1]*[v,1]'
        @test H ≈ QLPackedQ([NaN v; NaN NaN], [0,τ])
        Qt = [σ 0; 0 1] * H
        @test Qt ≈ F.Q'
        @test  (householderiterate(Qt, 11)' * diagm(0 => [-1; (-1).^(1:10)]))[1:10,1:10] ≈ Qn[1:10,1:10]
        # @test  (householderiterate(Qt, 11)' * diagm(0 => [1; 1; (-1).^(1:9)]))[1:10,1:10] ≈ Qn[1:10,1:10]

        t1, ω1 = Qn.τ[2], Qn.factors[1,2]
        Q1 = QLPackedQ([NaN ω1; NaN NaN], [0, t1])
        t2, ω2 = Qn.τ[3], Qn.factors[2,3]
        Q2 = QLPackedQ([NaN ω2; NaN NaN], [0, t2])

        @test abs.(householderiterate(Q2', Q1', 11)')[1:10,1:10] ≈ abs.(Qn[1:10,1:10])

        Q̃n = QLPackedQ(Qn.factors, [0; Qn.τ[2:end]])
        @test Q̃n[1:10,1:10] ≈ (householderiterate(Q2', Q1', 11)')[1:10,1:10]

        S1 = [-1 0; 0 1] * Qt
        @test S1 ≈ [-σ 0; 0 1] * H
        S2 = [1 0; 0 -1] * Qt * [1 0; 0 -1]
        @test S2 ≈ [σ 0; 0 1] * (I - τ *[-v,1]*[-v,1]')
        @test (diagm(0 => [-1; Ones(10)]) * householderiterate(S1, S2, 11))'[1:10,1:10] ≈ Qn[1:10,1:10]
        # @test (diagm(0 => [1; -Ones(10)]) * householderiterate(S1, S2, 11))'[1:10,1:10] ≈ Qn[1:10,1:10]

        s1 = sign(S1[1,1]/Q2[1,1]')
        s2 = sign(Q2[2,2]')
        @test [s1 0; 0 1] * Q2' * [1 0; 0 -1/s2] ≈ S1
        @test [-s2 0; 0 1] * Q1' * [1 0 ; 0 1/s1] ≈ S2

        @test (1-τ)*(-s2) ≈ 1-t2'  # [2,2] entry
        @test (1-τ)*s1 ≈ 1-t1'
        @test 1/s1 * (-σ) * (1-τ*abs2(v)) ≈ 1-t2' * abs2(ω2) # [1,1] entry
        @test -1/s2 * (σ) * (1-τ*abs2(v)) ≈ 1-t1' * abs2(ω1) # [1,1] entry
        @test 1/s1 * (-s2) * (σ*τ*v) ≈ -t2' * ω2
        @test τ*v ≈ t2'*ω2'
        @test -1/s2 * (s1) * (σ*τ*v) ≈ -t1' * ω1
        @test τ*v ≈ -t1'*ω1'

        @test (-σ) * (1-τ*abs2(v)) ≈ s1 - s2 * (σ*τ*v) * ω2' # remove abs2
        @test -σ * (1-τ*abs2(v)) ≈ s2- (s1) * (σ*τ*v) * ω1'
        @test -σ * (1-τ*abs2(v)) * t2' ≈ s1 * t2' - s2 * (σ*τ*v) * τ*v  # remove abs2
        @test -σ * (1-τ*abs2(v)) * t1' ≈ s2 * t1' + (s1) * (σ*τ*v) * τ*v

        @test -σ * (1-τ*abs2(v)) * (1-(1-τ)*(-s2)) ≈ s1 * (1-(1-τ)*(-s2)) - s2 * (σ*τ*v) * τ*v  # remove abs2
        @test -σ * (1-τ*abs2(v)) * (1-(1-τ)*s1) ≈ s2 * (1-(1-τ)*s1 ) + (s1) * (σ*τ*v) * τ*v

        α = -σ * (1-τ*abs2(v))
        β =  (α*(1-τ)+σ*(τ*v)^2)
        @test α + β*s2 ≈ s1 + (1-τ)*s1*s2   # remove abs2
        @test α - β*s1 ≈ s2  - (1-τ)*s1*s2

        @test α + β*s2 - s1 ≈ s2 - α + β*s1
        @test 2α + (β-1)*s2  ≈  (β+1)*s1

        @test (β+1)*α - 2α  + (β^2+1)*s2 ≈   (1-τ)*2α*s2 + (1-τ)*(β-1)*s2^2   # remove abs2
        a,b,c = (1-τ)*(β-1), (1-τ)*2α-(β^2+1), 2α- (β+1)*α
        @test a*s2^2 + b*s2 ≈ -c

        @test s2 ≈ (-b - sqrt(b^2-4a*c))/(2a)
        @test s1 ≈ (2α + (β-1)*s2)/(β+1)
        @test t1' ≈ 1-(1-τ)*s1
        @test ω1' ≈ -τ*v/t1'
        @test t2' ≈ 1-(1-τ)*(-s2)
        @test ω2' ≈ τ*v / t2'

        X = F.factors
        m = size(X,2)
        data = Hcat([-X[1,end-1]; (-1).^(1:m-1) .* X[2,end-1:-1:1]],
                    ApplyArray(*, ((-1).^(1:m) .* X[2,end:-1:1]), ((-1).^(0:∞))'))
        L∞ = _BandedMatrix(data,∞,m-1,0)
        @test L∞[1:10,1:10] ≈ Ln[1:10,1:10]

        Q∞11 = 1 - ω1*t1*conj(ω1)  # Q[1,1] without the callowing correction
        Q̃n = QLPackedQ(Qn.factors, [0; Qn.τ[2:end]])
        @test Q∞11 ≈ Q̃n[1,1]

        data2 = Vcat(Hcat(0.0+0im, mortar(Fill([ω1 ω2],1,∞))), data)
        factors = _BandedMatrix(data2,∞,m-1,1)
        Q̃∞ = QLPackedQ(factors, Vcat(0.0+0im,mortar(Fill([t1,t2],∞))))
        @test Q̃∞[1:10,1:10] ≈ Q̃n[1:10,1:10]

        @test Q̃∞[1:10,1:11]*L∞[1:11,2:10] ≈ A[1:10,2:10]

        @test (1-Qn.τ[1])*Q̃∞[1,1] * L∞[1,1] + Q̃∞[1,2] * L∞[2,1] ≈ A[1,1]
        L∞11 = -X[1,end-1]
        L∞21 = -X[2,end-1]
        Q∞12 = -t1*ω1
        @test Qn.τ[1] ≈ (Q∞12 * L∞21 - A[1,1])/(L∞11*Q∞11) + 1

        Q∞, L∞ = ql(A)
        @test Q∞[1:10,1:11] * L∞[1:11,1:10] ≈ A[1:10,1:10]
    end
    @testset "3-periodic" begin
        ain =    [ -2.531640004434771-0.0im , 0.36995310821558014+2.5612894011525276im, -0.22944284364953327+0.39386202384951985im, -0.2700241133710857 + 0.8984628598798804im, 4.930380657631324e-32 + 0.553001215633963im ]
        l = 3
        A = _BandedMatrix(ain * Ones{ComplexF64}(1,∞), ∞, l, 1)
        n = 1000; Qn,Ln = ql(A[1:n,1:n])
        a = reverse(A.data.args[1])
        de = tail_de(a)
        X =  [transpose(a); 0 transpose(de)]
        F = ql_X!(copy(X))
        @test F.factors[1,1:end-1] ≈ de
        @test F.factors[2,:] ≈ Ln[l+2,1:l+2]
        @test F.factors[2,:] ≈ -Ln[l+3,2:l+3]

        (σ,τ,v) = householderparams(F)
        H = I - τ *[v,1]*[v,1]'
        @test H ≈ QLPackedQ([NaN v; NaN NaN], [0,τ])
        Qt = [σ 0; 0 1] * H
        @test Qt ≈ F.Q'
        @test  (householderiterate(Qt, 11)' * diagm(0 => [1; mortar(Fill([1,-1,-1],4))[1:10]]))[1:10,1:10] ≈ Qn[1:10,1:10]

        t1, ω1 = Qn.τ[2], Qn.factors[1,2]
        Q1 = QLPackedQ([NaN ω1; NaN NaN], [0, t1])
        t2, ω2 = Qn.τ[3], Qn.factors[2,3]
        Q2 = QLPackedQ([NaN ω2; NaN NaN], [0, t2])
        t3, ω3 = Qn.τ[4], Qn.factors[3,4]
        Q3 = QLPackedQ([NaN ω3; NaN NaN], [0, t3])

        @test abs.(householderiterate(Q3', Q2', Q1', 11)')[1:10,1:10] ≈ abs.(Qn[1:10,1:10])
        @test abs.(householderiterate(Q1', Q2', Q3', 11)')[1:10,1:10] ≈ abs.(Qn[1:10,1:10])
        @test abs.(householderiterate(Q2', Q3', Q1', 11)')[1:10,1:10] ≈ abs.(Qn[1:10,1:10])

        # remove last
        Q̃n = QLPackedQ(Qn.factors, [0; Qn.τ[2:end]])
        @test Q̃n[1:10,1:10] ≈ (householderiterate(Q2', Q3', Q1', 15)')[1:10,1:10]

        S1 = [1 0; 0 -1] * Qt * [1 0; 0 -1]
        @test S1 ≈ [σ 0; 0 1] * (I - τ *[-v,1]*[-v,1]')
        S2 = [-1 0; 0 -1] * Qt * [1 0; 0 -1]
        @test S2 ≈ [-σ 0; 0 1] * (I - τ *[-v,1]*[-v,1]')
        S3 = [-1 0; 0 1] * Qt
        @test S3 ≈ [-σ 0; 0 1] * (I - τ *[v,1]*[v,1]')
        @test (diagm(0 => [-1; Ones(10)]) * householderiterate(S1, S2, S3, 11))'[1:10,1:10] ≈ Qn[1:10,1:10]

        s1 = sign(S1[1,1]/Q2[1,1]')
        s2 = sign(Q2[2,2]')
        s3 = sign(Q3[2,2]')
        @test [s1 0; 0 1] * Q2' * [1 0; 0 -1/s2] ≈ S1
        @test [-s2 0; 0 1] * Q3' * [1 0 ; 0 -1/s3] ≈ S2
        @test [-s3 0; 0 1] * Q1' * [1 0 ; 0 1/s1] ≈ S3

        @test (1-τ)*(-s2) ≈ 1-t2'  # [2,2] entry (1)
        @test (1-τ)*s1 ≈ 1-t1'
        @test (1-τ)*(-s3) ≈ 1-t3'
        @test -1/s1 * (-σ) * (1-τ*abs2(v)) ≈ 1-t2' * abs2(ω2) # [1,1] entry (4)
        @test 1/s2 * (σ) * (1-τ*abs2(v)) ≈ 1-t3' * abs2(ω3) # [1,1] entry
        @test 1/s3 * (σ) * (1-τ*abs2(v)) ≈ 1-t1' * abs2(ω1) # [1,1] entry
        @test 1/s1 * (-s2) * (σ*τ*v) ≈ -t2' * ω2 # (7)
        @test 1/s2 * (-s3) * (σ*τ*v) ≈ -t3' * ω3
        @test 1/s3 * (-s1) * (σ*τ*v) ≈ -t1' * ω1
        @test τ*v ≈ -t2'*ω2' # (10)
        @test τ*v ≈ t1'*ω1'
        @test τ*v ≈ -t3'*ω3'

        @test σ * (1-τ*abs2(v)) ≈ s1 - s2 * (σ*τ*v) * ω2'  # (13) = (4) & (7)
        @test σ * (1-τ*abs2(v)) ≈ s2 - s3 * (σ*τ*v) * ω3'  # (14) = (5) & (8)
        @test σ * (1-τ*abs2(v)) ≈ s3 - s1 * (σ*τ*v)  * ω1' # (15) = (6) & (9)

        @test (τ*v)/(-(1-(1-τ)*(-s2) )) ≈ ω2' # (16) = (10) & (1)
        @test (τ*v)/((1-(1-τ)*(s1) )) ≈ ω1'
        @test (τ*v)/(-(1-(1-τ)*(-s3) )) ≈ ω3'

        @test σ * (1-τ*abs2(v)) ≈ s1 - s2 * (σ*τ*v) * (τ*v)/(-(1-(1-τ)*(-s2) )) #  (19) = (16) & (13)
        @test σ * (1-τ*abs2(v)) ≈ s2 - s3 * (σ*τ*v) * (τ*v)/(-(1-(1-τ)*(-s3) ))
        @test σ * (1-τ*abs2(v)) ≈ s3 - s1 * (σ*τ*v) * (τ*v)/((1-(1-τ)*(s1) ))

        α = σ * (1-τ*abs2(v))
        β = α - τ*σ
        @test α  + β*s2  ≈ s1 + (1-τ)*s1*s2
        @test  α + β*s3  ≈ s2 + (1-τ)*s2*s3
        @test  α - β*s1  ≈ s3 - (1-τ)*s1*s3

        # from Mathematica
        @test s1 ≈ -((1 + 3α + 3*α*β + β^3 - 3*α*τ - 3*α*β*τ - sqrt((1 + β^3 - 3*α*(1 + β)*(-1 + τ))^2 + 4*α*(-1 + τ)*(1 + α + β + β^2 - α*τ)^2))/(2*(-1 + τ)*(1 + α + β + β^2 - α*τ)))
        @test s2 ≈ -((1 - α - α*β + β^3 + α*τ + α*β*τ - sqrt((1 + β^3 - 3*α*(1 + β)*(-1 + τ))^2 + 4*α*(-1 + τ)*(1 + α + β + β^2 - α*τ)^2))/(2*(-1 + β + β^2 + α*(-1 + τ))*(-1 + τ)))
        @test s3 ≈ (1 - α - α*β + β^3 + α*τ + α*β*τ - sqrt((1 + β^3 - 3*α*(1 + β)*(-1 + τ))^2 + 4*α*(-1 + τ)*(1 + α + β + β^2 - α*τ)^2))/(2*(1 + β - β^2 + α*(-1 + τ))*(-1 + τ))

        @test t2' ≈ 1-(1-τ)*(-s2) # (1)
        @test t1' ≈ 1-(1-τ)*s1 # (2)
        @test t3' ≈ 1+(1-τ)*s3 # (3)

        @test ω1' ≈ τ*v/t1'  # (11)
        @test ω2' ≈ -τ*v / t2'
        @test ω3' ≈ -τ*v / t3'

        X = F.factors
        m = size(X,2)
        μ = mortar(Fill([1 -1 -1], 1, ∞))
        μ0 = mortar(Fill([1 0 0], 1, ∞))
        μ1 = mortar(Fill([0 1 0], 1, ∞))
        μ2 = mortar(Fill([0 0 1], 1, ∞))

        data = Hcat([X[1,end-1]; μ[1,1:m-1] .* X[2,end-1:-1:1]], ApplyArray(*, [μ[1,1:m].*X[2,end:-1:1] μ[1,2:m+1].*X[2,end:-1:1] μ[1,3:m+2].*X[2,end:-1:1]], Vcat(μ0,μ1,μ2)))
        L∞ = _BandedMatrix(data,∞,m-1,0)
        @test L∞[1:10,1:10] ≈ Ln[1:10,1:10]

        Q∞11 = 1 - ω1*t1*conj(ω1)  # Q[1,1] without the callowing correction
        Q̃n = QLPackedQ(Qn.factors, [0; Qn.τ[2:end]])
        @test Q∞11 ≈ Q̃n[1,1]

        data2 = Vcat(Hcat(0.0+0im, mortar(Fill([ω1 ω2 ω3],1,∞))), data)
        factors = _BandedMatrix(data2,∞,m-1,1)
        Q̃∞ = QLPackedQ(factors, Vcat(0.0+0im,mortar(Fill([t1,t2,t3],∞))))
        @test Q̃∞[1:10,1:10] ≈ Q̃n[1:10,1:10]

        @test Q̃∞[1:10,1:11]*L∞[1:11,2:10] ≈ A[1:10,2:10]

        @test (1-Qn.τ[1])*Q̃∞[1,1] * L∞[1,1] + Q̃∞[1,2] * L∞[2,1] ≈ A[1,1]
        L∞11 = X[1,end-1]
        @test L∞11 ≈ Ln[1,1]
        L∞21 = X[2,end-1]
        @test L∞21 ≈ Ln[2,1]
        Q∞12 = -t1*ω1
        @test Q∞12 ≈ Qn[1,2]
        @test Qn.τ[1] ≈ (Q∞12 * L∞21 - A[1,1])/(L∞11*Q∞11) + 1

        Q∞, L∞ = ql(A)
        @test Q∞[1:10,1:11] * L∞[1:11,1:10] ≈ A[1:10,1:10]
    end
end

@testset "Toeplitz" begin
    for (Z,A,B) in ((2,2.1+0.01im,0.5),) # (2,-0.1+0.1im,0.5))
        n = 1_000; T = Tridiagonal(Fill(ComplexF64(Z),∞), Fill(ComplexF64(A),∞), Fill(ComplexF64(B),∞)); Q,L = ql(BandedMatrix(T)[1:n,1:n]);
        d,e = tail_de([Z,A,B])
        @test L[1,1] ≈ e
        @test (Q'*[Z; zeros(n-1)])[1] ≈ d
        @test ql([Z A B; 0 d e]).L[1,1:2] ≈ [d;e]
        @test ql([Z A B; 0 d e]).L[2,:] ≈ -L[3,1:3]

        F = ql_X!([Z A B; 0 d e])
        @test F.Q'*[Z A B; 0 d e] ≈ F.L

        σ,τ,v = householderparams(F)
        @test [σ 0; 0 1] * (I - τ*[v,1]*[v,1]') ≈  F.Q'

        ql!([Z A B; 0 d e]).Q'
        householderiterate(F.Q', 100)[1:10,1:10]
        [σ/sqrt(σ) 0; 0 sqrt(σ)] * (I - τ*[v,1]*[v,1]') * [sqrt(σ) 0 ; 0 1/sqrt(σ)]

        s,t,ω = combine_two_Q(σ,τ,v)
        @test Q.τ[2] ≈ t
        @test Q.factors[1,2] ≈ ω

        Q∞,L∞ = F = ql(T)

        @test Q.τ[2] ≈ F.τ[2] ≈ t
        @test Q.factors[1,2] ≈ F.factors[1,2] ≈ ω
        @test Q∞[1:10,1:10] ≈ Q[1:10,1:10]
        @test L∞[1:10,1:10] ≈ L[1:10,1:10]
        @test Q∞[1:10,1:12] * L∞[1:12,1:10] ≈ T[1:10,1:10]
    end
end

@testset "QL Construction" begin
    c,a,b =  (2,0.0+0.0im,0.5)
    J = Tridiagonal(Vcat(ComplexF64[], Fill(ComplexF64(c),∞)),
                        Vcat(ComplexF64[2], Fill(ComplexF64(a),∞)),
                        Vcat(ComplexF64[], Fill(ComplexF64(b),∞)))
    A = J + (2.1+0.01im)I
    B = BandedMatrix(A,(2,1))
    T = toeptail(B)
    F∞ = ql(T)
    σ = 1- F∞.τ[1]
    F∞.factors[1,1] *= σ
    F∞ = QL(F∞.factors, Vcat(zero(ComplexF64),F∞.τ.args[2]))
    Q∞, L∞ = F∞
    @test Q∞[1:10,1:12] * L∞[1:12,1:10] ≈ T[1:10,1:10]

    n = 100_000; Q,L = ql(B[1:n,1:n]);
    @test Q.τ[3] ≈ F∞.τ[2]

    data = bandeddata(B).args[1]
    e = F∞.factors[1,1]
    d = conj(F∞.Q[1,1]*c)
    @test d == (F∞.Q'*Vcat(c,Zeros(∞)))[1]
    B̃ = _BandedMatrix(data, size(data,2), 2,1)
    B̃[end,end-1:end] .= (d,e)
    F = ql!(B̃)
    @test Q.τ[1] ≈ F.τ[1]
    B̃.data[3:end,end] .= (L∞[2,1], L∞[3,1]) # fill in L
    B̃.data[4,end-1] = L∞[3,1] # fill in L
    H = Hcat(B̃.data, F∞.factors[1:4,2] * Ones{eltype(F∞)}(1,∞))
    @test Q.τ[1:10] ≈ Vcat(F.τ, F∞.τ.args[2])[1:10]

    Q2,L2 = ql(A)

    @test Q.τ[1:10] ≈ Q2.τ[1:10]
    @test Q2[1:10,1:12]*L2[1:12,1:10] ≈ A[1:10,1:10]
end

@testset "Hessenberg Toeplitz" begin
    a = [1,2,3,0.5]
    T = _BandedMatrix(reverse(a) * Ones(1,∞), ∞, 2, 1)
    n = 10_000; Q, L = ql(T[1:n,1:n])

    de = tail_de(a)
    @test L[1,1] ≈ de[end]
    @test vec((Q')[1:1,1:2]*T[3:4,1:2]) ≈ de[1:end-1]

    X = [transpose(a); [0 transpose(de)]]
    @test ql(X).L[1,1:3] ≈ de
    @test ql(X).τ[1] == 0
    @test ql(X).τ[2] ≈ Q.τ[2]

    @test T isa InfToeplitz

    T = BandedMatrix(-2 => Fill(1,∞), 0 => Fill(0.5+eps()im,∞), 1 => Fill(0.25,∞))
    Q,L = QL(ql(T))
    Qn,Ln = ql(T[1:1000,1:1000])
    @test Qn.τ[1:10] ≈ Q.τ[1:10]
    @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]
end

function Toep_L11(T)
    l,u = bandwidths(T)
    @assert u == 2
    # shift by one
    H = _BandedMatrix(T.data, ∞, l+1, 1)
    Q1,L1 = ql(H)

    d = Q1[1:3,1]'T[1:1+l,1]
    ℓ = Q1.factors.data.args[2].args[1][2:end] # new L
    T2 = _BandedMatrix(Hcat([[zero(d); d; ℓ[3:end]] L1[1:5,1]], ℓ*Ones{eltype(T)}(1,∞)), ∞, 3, 1)
    Q2,L2 = ql(T2)
    D = (Q2')[1:5,1:4] * (Q1')[1:4,1:3] * T[3:5,1:3]
    X = [Matrix(T[3:4,1:6]); [zeros(2,2) [-1 0; 0 1]*D[1:2,:] [-1 0; 0 1]*L2[1:2,2]]]
    ql(X).L[1,3]
end

@testset "Pentadiagonal Toeplitz" begin
    a = [1,2,5+0im,0.5,0.2]
    T = _BandedMatrix(reverse(a) * Ones{eltype(a)}(1,∞), ∞, 2, 2)

    n = 1_000; Q, L = ql(T[1:n,1:n])
    H = _BandedMatrix(T.data, ∞, 3, 1)
    F1 = ql(H)
    Q1,L1 = F1
    @test Q1[1:10,1:11]*L1[1:11,1:10] ≈ H[1:10,1:10]

    d = Q1[1:3,1]'T[1:1+T.l,1]
    ℓ = F1.factors.data.args[2].args[1][2:end]
    @test Vcat(zero(d), d, ℓ[3:end])[2:end] ≈ (Q1')[1:4,1:3]*T[1:3,1]

    T2 = _BandedMatrix(Hcat([[zero(d); d; ℓ[3:end]] L1[1:5,1]], ℓ*Ones{eltype(a)}(1,∞)), ∞, 3, 1)
    @test T2 isa PertToeplitz
    Q2,L2 = ql(T2)
    @test Q2[1:10,1:11]*L2[1:11,1:10] ≈ T2[1:10,1:10]
    @test Q1[1:10,1:11]*Q2[1:11,1:12]*L2[1:12,1:10] ≈ T[1:10,1:10]

    D = (Q2')[1:5,1:4] * (Q1')[1:4,1:3] * T[3:5,1:3]
    X = real([Matrix(T[3:4,1:6]); [zeros(2,2) [-1 0; 0 1]*D[1:2,:] [-1 0; 0 1]*L2[1:2,2]]])
    @test ql(X).L[1:2,1:4] ≈ X[3:4,3:end]

    @test L[1,1] ≈ -ql(X).L[1,3]

    @test abs(L[1,1]) ≈ abs(Toep_L11(T))

    a = [1,2,0,0.5,0.2]
    T = _BandedMatrix(reverse(a) * Ones{eltype(a)}(1,∞), ∞, 2, 2)
    A = T+(5+0im)I
    @test abs(ql((A)[1:1000,1:100]).L[1,1]) ≈ abs(Toep_L11(A))
    A = T+(-5+0im)I
    @test abs(ql((A)[1:1000,1:100]).L[1,1]) ≈ abs(Toep_L11(A))

    A = T+(-5+1im)I
    @test abs(ql((A)[1:1000,1:100]).L[1,1]) ≈ abs(Toep_L11(A))
end

@testset "3-diagonals" begin
    @testset "No periodic" begin
        A = BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(5,∞), -1 => Fill(2im,∞))
        l,u = bandwidths(A)
        H = _BandedMatrix(A.data, ∞, l+u-1, 1)
        Q1,L1 = ql(H)
        @test Q1[1:10,1:11] * L1[1:11,1:10] ≈ H[1:10,1:10]
        @test L1[1:10,1:10] ≈ Q1[1:13,1:10]'H[1:13,1:10]

        @test (Q1[1:13,1:10]'A[1:13,1:12])[1:10,u:10+u-1] ≈ L1[1:10,1:10]
        # to the left of H
        D1, Q1, L1 = reduceband(A)
        T2 = _BandedMatrix(rightasymptotics(parent(L1).data).args[1][2:end] * Ones{ComplexF64}(1,∞), ∞, l, u)
        l1 = L1[1,1]

        A2 = [[D1 l1 zeros(1,10-size(D1,2)-1)]; T2[1:10-1,1:10]]
        @test Q1[1:13,1:10]'A[1:13,1:10] ≈ A2

        B2 = _BandedMatrix(T2.data, ∞, l+u-2, 2)
        D2, Q2, L2 = reduceband(B2)
        l2 = L2[1,1]
        T3 = _BandedMatrix(rightasymptotics(parent(L2).data).args[1][2:end] * Ones{ComplexF64}(1,∞), ∞, l+1, u-1)
        A3 = [[D2 l2 zeros(1,10-size(D2,2)-1)]; T3[1:10-1,1:10]]
        @test Q2[1:13,1:10]'B2[1:13,1:10] ≈ A3

        n,m = 10,10
        Q̃2 = [1 zeros(1,m-1);  zeros(n-1,1) Q2[1:n-1,1:m-1]]
        @test norm((Q̃2'A2[:,1:end-2])[band(2)][2:end]) ≤ 10eps() # banded apart from (1,1) entry
        @test (Q̃2'A2[:,1:end-2])[2:end,2:end] ≈ A3[1:9,1:7]
        @test (Q̃2'A2[:,1:end-2])[3:end,1] ≈ A3[3:10,1]
        @test (Q̃2'A2[:,1:end-2])[1:5,1:2] ≈  Q̃2[1:4,1:5]' * [D1; T2[1:size(D1,2)+1,1:2]]
        @test (Q̃2'A2[:,1:end-2]) ≈ [A2[1,1] A2[1:1,2:8]; [Q2[1:3,1:3]' * T2[1:3,1]; Zeros(10-4)]  A3[1:end-1,1:7] ]

        # fix last entry
        @test (Q̃2'A2[:,1:3])[1:2,1:3] ≈ [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]]
        Q3,L3 = ql( [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]])
        Q̃3 = [Q3 zeros(2, n-2); zeros(n-2,2) I]
        @test norm((Q̃3'Q̃2'A2[:,1:end-2])[band(2)]) ≤ 10eps()

        @test Q̃3'Q̃2'A2[:,1:end-2] ≈ [L3 zeros(2,8-3); [[Q2[1:3,2:3]' * T2[1:3,1]; Zeros(10-4)] A3[2:end-1,1:7] ] ]

        fd_data = hcat([0; L3[:,1]; Q2[1:3,2:3]' * T2[1:3,1]], [L3[:,2]; T3[1:3,1]], [L3[2,3]; T3[1:4,2]])
        B3 = _BandedMatrix(Hcat(fd_data, T3.data), ∞, l+u-1, 1)
        @test B3[1:10,1:8] ≈ Q̃3'Q̃2'A2[:,1:end-2]

        @test ql(B3).L[1,1] ≈ ql(A[1:1000,1:1000]).L[1,1]
    end
    @testset "periodic" begin
        A = BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(0.5-0.1im,∞), -1 => Fill(-2im,∞))
        l,u = bandwidths(A)
        H = _BandedMatrix(A.data, ∞, l+u-1, 1)
        Q1,L1 = ql(H)
        @test Q1[1:10,1:11] * L1[1:11,1:10] ≈ H[1:10,1:10]
        @test L1[1:10,1:10] ≈ Q1[1:13,1:10]'H[1:13,1:10]

        @test (Q1[1:13,1:10]'A[1:13,1:12])[1:10,u:10+u-1] ≈ L1[1:10,1:10]
        # to the left of H
        D1, Q1, L1 = reduceband(A)
        T2 = _BandedMatrix(rightasymptotics(parent(L1).data).args[1][2:end] * Ones{ComplexF64}(1,∞), ∞, l, u)
        l1 = L1[1,1]

        A2 = [[D1 l1 zeros(1,10-size(D1,2)-1)]; T2[1:10-1,1:10]]
        @test Q1[1:13,1:10]'A[1:13,1:10] ≈ A2

        B2 = _BandedMatrix(T2.data, ∞, l+u-2, 2)
        D2, Q2, L2 = reduceband(B2)
        l2 = L2[1,1]

        # peroidic tail
        T3 = _BandedMatrix(rightasymptotics(parent(L2).data).args[2], ∞, l+1, u-1)
        A3 = [[D2 l2 zeros(1,10-size(D2,2)-1)]; T3[1:10-1,1:10]]
        @test Q2[1:13,1:10]'B2[1:13,1:10] ≈ A3

        n,m = 10,10
        Q̃2 = [1 zeros(1,m-1);  zeros(n-1,1) Q2[1:n-1,1:m-1]]
        @test norm((Q̃2'A2[:,1:end-2])[band(2)][2:end]) ≤ 20eps() # banded apart from (1,1) entry
        @test (Q̃2'A2[:,1:end-2])[2:end,2:end] ≈ A3[1:9,1:7]
        @test (Q̃2'A2[:,1:end-2])[3:end,1] ≈ -A3[3:10,1]
        @test (Q̃2'A2[:,1:end-2])[1:5,1:2] ≈  Q̃2[1:4,1:5]' * [D1; T2[1:size(D1,2)+1,1:2]]
        @test (Q̃2'A2[:,1:end-2]) ≈ [A2[1,1] A2[1:1,2:8]; [Q2[1:3,1:3]' * T2[1:3,1]; Zeros(10-4)]  A3[1:end-1,1:7] ]

        # fix last entry
        @test (Q̃2'A2[:,1:3])[1:2,1:3] ≈ [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]]
        Q3,L3 = ql( [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]])
        Q̃3 = [Q3 zeros(2, n-2); zeros(n-2,2) I]
        @test norm((Q̃3'Q̃2'A2[:,1:end-2])[band(2)]) ≤ 50eps()

        @test Q̃3'Q̃2'A2[:,1:end-2] ≈ [L3 zeros(2,8-3); [[Q2[1:3,2:3]' * T2[1:3,1]; Zeros(10-4)] A3[2:end-1,1:7] ] ]

        fd_data = hcat([0; L3[:,1]; Q2[1:3,2:3]' * T2[1:3,1]], [L3[:,2]; T3[1:3,1]], [L3[2,3]; T3[1:4,2]])
        B3 = _BandedMatrix(Hcat(fd_data, T3.data), ∞, l+u-1, 1)
        @test B3[1:10,1:8] ≈ Q̃3'Q̃2'A2[:,1:end-2]

        # remove oscillation

        fd_data_s = diagm(0 => (-1).^(0:size(fd_data,1)-1)) * (fd_data * diagm(0 => (-1).^(1:size(fd_data,2))))
        T3_data_s = (-1)^size(fd_data,2) * (-1).^(1:u+l+1) .* T3.data[:,1]
        B3_s = _BandedMatrix(Hcat(fd_data_s, T3_data_s*Ones{ComplexF64}(1,∞)), ∞, l+u-1, 1)

        @test diagm(0 => (-1).^(0:9)) * B3[1:10,1:10] ≈ B3_s[1:10,1:10]

        @test ql(B3_s).L[1,1] ≈ ql(A[1:1000,1:1000]).L[1,1]
    end

    @testset "3rd case" begin
        A = BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(3-0.1im,∞), -1 => Fill(-2im,∞))
        l,u = bandwidths(A)
        H = _BandedMatrix(A.data, ∞, l+u-1, 1)
        Q1,L1 = ql(H)
        @test Q1[1:10,1:11] * L1[1:11,1:10] ≈ H[1:10,1:10]
        @test L1[1:10,1:10] ≈ Q1[1:13,1:10]'H[1:13,1:10]

        @test (Q1[1:13,1:10]'A[1:13,1:12])[1:10,u:10+u-1] ≈ L1[1:10,1:10]
        # to the left of H
        D1, Q1, L1 = reduceband(A)
        T2 = _BandedMatrix(rightasymptotics(parent(L1).data).args[1][2:end] * Ones{ComplexF64}(1,∞), ∞, l, u)
        l1 = L1[1,1]

        A2 = [[D1 l1 zeros(1,10-size(D1,2)-1)]; T2[1:10-1,1:10]]
        @test Q1[1:13,1:10]'A[1:13,1:10] ≈ A2

        B2 = _BandedMatrix(T2.data, ∞, l+u-2, 2)
        D2, Q2, L2 = reduceband(B2)
        l2 = L2[1,1]

        # peroidic tail
        T3 = _BandedMatrix(rightasymptotics(parent(L2).data).args[1][2:end] * Ones{ComplexF64}(1,∞), ∞, l+1, u-1)
        A3 = [[D2 l2 zeros(1,10-size(D2,2)-1)]; T3[1:10-1,1:10]]
        @test Q2[1:13,1:10]'B2[1:13,1:10] ≈ A3

        n,m = 10,10
        Q̃2 = [1 zeros(1,m-1);  zeros(n-1,1) Q2[1:n-1,1:m-1]]
        @test norm((Q̃2'A2[:,1:end-2])[band(2)][2:end]) ≤ 20eps() # banded apart from (1,1) entry
        @test (Q̃2'A2[:,1:end-2])[2:end,2:end] ≈ A3[1:9,1:7]
        @test (Q̃2'A2[:,1:end-2])[3:end,1] ≈ A3[3:10,1] # swapped sign
        @test (Q̃2'A2[:,1:end-2])[1:5,1:2] ≈  Q̃2[1:4,1:5]' * [D1; T2[1:size(D1,2)+1,1:2]]
        @test (Q̃2'A2[:,1:end-2]) ≈ [A2[1,1] A2[1:1,2:8]; [Q2[1:3,1:3]' * T2[1:3,1]; Zeros(10-4)]  A3[1:end-1,1:7] ]

        # fix last entry
        @test (Q̃2'A2[:,1:3])[1:2,1:3] ≈ [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]]
        Q3,L3 = ql( [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]])
        Q̃3 = [Q3 zeros(2, n-2); zeros(n-2,2) I]
        @test norm((Q̃3'Q̃2'A2[:,1:end-2])[band(2)]) ≤ 50eps()

        @test Q̃3'Q̃2'A2[:,1:end-2] ≈ [L3 zeros(2,8-3); [[Q2[1:3,2:3]' * T2[1:3,1]; Zeros(10-4)] A3[2:end-1,1:7] ] ]

        fd_data = hcat([0; L3[:,1]; Q2[1:3,2:3]' * T2[1:3,1]], [L3[:,2]; T3[1:3,1]], [L3[2,3]; T3[1:4,2]])
        B3 = _BandedMatrix(Hcat(fd_data, T3.data), ∞, l+u-1, 1)
        @test B3[1:10,1:8] ≈ Q̃3'Q̃2'A2[:,1:end-2]

        # test tail
        T4 = _BandedMatrix(T3.data, ∞, l+u-1,1)
        Q4, L4 = ql(T4)
        @test Q4[1:10,1:11]*L4[1:11,1:10] ≈ T4[1:10,1:10]

        @test ql(B3).L[1,1] ≈ ql(B3[1:10000,1:10000]).L[1,1]
        @test abs(ql(B3).L[1,1]) ≈ abs(ql(A[1:10000,1:10000]).L[1,1])
    end
end

_Lrightasymptotics(D::Vcat) = D.args[2]
_Lrightasymptotics(D::ApplyArray) = D.args[1][2:end] * Ones{ComplexF64}(1,∞)
Lrightasymptotics(L) = _Lrightasymptotics(rightasymptotics(parent(L).data))

function qdL(A)
    l,u = bandwidths(A)
    H = _BandedMatrix(A.data, ∞, l+u-1, 1)
    Q1,L1 = ql(H)
    D1, Q1, L1 = reduceband(A)
    T2 = _BandedMatrix(Lrightasymptotics(L1), ∞, l, u)
    l1 = L1[1,1]
    A2 = [[D1 l1 zeros(1,10-size(D1,2)-1)]; T2[1:10-1,1:10]] # TODO: remove
    B2 = _BandedMatrix(T2.data, ∞, l+u-2, 2)
    D2, Q2, L2 = reduceband(B2)
    l2 = L2[1,1]
    # peroidic tail
    T3 = _BandedMatrix(Lrightasymptotics(L2), ∞, l+1, u-1)
    A3 = [[D2 l2 zeros(1,10-size(D2,2)-1)]; T3[1:10-1,1:10]] # TODO: remove

    Q3,L3 = ql( [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]])

    fd_data = hcat([0; L3[:,1]; Q2[1:3,2:3]' * T2[1:3,1]], [L3[:,2]; T3[1:3,1]], [L3[2,3]; T3[1:4,2]])
    B3 = _BandedMatrix(Hcat(fd_data, T3.data), ∞, l+u-1, 1)

    if T3 isa InfToeplitz
        ql(B3).L
    else
        # remove periodicity
        fd_data_s = diagm(0 => (-1).^(0:size(fd_data,1)-1)) * (fd_data * diagm(0 => (-1).^(1:size(fd_data,2))))
        T3_data_s = (-1)^size(fd_data,2) * (-1).^(1:u+l+1) .* T3.data[:,1]
        B3_s = _BandedMatrix(Hcat(fd_data_s, T3_data_s*Ones{ComplexF64}(1,∞)), ∞, l+u-1, 1)
        ql(B3_s).L
    end
end

@testset "quick-and-dirty L" begin
    for λ in (5,1,0.1+0.1im,-0.5-0.1im), A in (BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(-λ,∞), -1 => Fill(2im,∞)),
                                        BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(-conj(λ),∞), -1 => Fill(-2im,∞)))
        L∞ = qdL(A)[1:10,1:10]
        Ln = ql(A[1:1000,1:1000]).L[1:10,1:10]
        @test L∞ .* sign.(diag(L∞)) ≈ Matrix(Ln) .* sign.(diag(Ln))
    end
    for λ in (-3-0.1im, 0.0, -1im)
        @show λ
        A = BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(-conj(λ),∞), -1 => Fill(-2im,∞))
        @test abs(qdL(A)[1,1]) ≈ abs(ql(A[1:10000,1:10000]).L[1,1])
    end
    for λ in (1+2im,)
        A = BandedMatrix(3 => Fill(7/10,∞), 2 => Fill(1,∞), 0 => Fill(-λ,∞), -1 => Fill(2im,∞))
        @test_throws DomainError qdL(A)
    end
end

@testset "bi-infinite" begin
    Δ = BandedMatrix(-2 => Vcat(Float64[],Fill(1.0,∞)), -1 =>Vcat([1.0], Fill(0.0,∞)), 1 => Vcat([1.0], Fill(0.0,∞)), 2 => Vcat(Float64[],Fill(1.0,∞)))
    A = (Δ - 4I)
    B = BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2)))
    T = toeptail(B)
    H = _BandedMatrix(T.data, ∞, bandwidth(T,1)+bandwidth(T,2)-1, 1)
    Q,L = ql(H)

    A2 = Q'A
    L
end


text:
#=
---------------------------------------------------------------------
This function computes the Metric Terms for General 2D Quad Grids.

The strategy follows Algorithm 12.3 in F.X. Giraldo's Introduction to Element-based
Galerkin Methods using Tensor-Product Bases: Analysis, Algorithms, and Applications.

Written by F.X. Giraldo on 5/2019
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

include("map_deriv.jl")

function compute_metrics(coord,intma,ψ,dψ,Ne,Np,Nq,DFloat)

    #Initialize Global Arrays
    ξ_x=zeros(DFloat,Nq,Nq,Ne)
    ξ_y=zeros(DFloat,Nq,Nq,Ne)
    η_x=zeros(DFloat,Nq,Nq,Ne)
    η_y=zeros(DFloat,Nq,Nq,Ne)
    jac=zeros(DFloat,Nq,Nq,Ne)

    #Initialize Local Arrays
    x=zeros(DFloat,Np,Np)
    y=zeros(DFloat,Np,Np)

    #loop thru the elements
    for e=1:Ne

        #Store Element Variables
        for j=1:Np
            for i=1:Np
                I=intma[i,j,e]
                x[i,j]=coord[1,I]
                y[i,j]=coord[2,I]
            end #i
        end #j

        #Construct Mapping Derivatives: dx/dξ; dx/dη; dy/dξ; dy/dη
        (x_ξ,x_η)=map_deriv(ψ,dψ,x,Np,Nq,DFloat)
        (y_ξ,y_η)=map_deriv(ψ,dψ,y,Np,Nq,DFloat)

        #Construct Inverse Mapping: dξ/dx; dξ/dy; dη/dx; dη/dy
        for j=1:Nq
            for i=1:Nq
                xjac=x_ξ[i,j]*y_η[i,j] - x_η[i,j]*y_ξ[i,j]
                ξ_x[i,j,e]=+1.0/xjac*y_η[i,j]
                ξ_y[i,j,e]=-1.0/xjac*x_η[i,j]
                η_x[i,j,e]=-1.0/xjac*y_ξ[i,j]
                η_y[i,j,e]=+1.0/xjac*x_ξ[i,j]
                jac[i,j,e]=abs(xjac)
            end #i
        end #j
    end #e

    return (ξ_x,ξ_y,η_x,η_y,jac)
end


text:
using DFTK
using JLD2
import DFTK: ScfDefaultCallback, ScfSaveCheckpoints
using Test
include("testcases.jl")

function test_scfres_agreement(tested, ref)
    @test tested.basis.model.lattice           == ref.basis.model.lattice
    @test tested.basis.model.smearing          == ref.basis.model.smearing
    @test tested.basis.model.symmetries        == ref.basis.model.symmetries
    @test tested.basis.model.spin_polarization == ref.basis.model.spin_polarization

    @test length(tested.basis.model.atoms) == length(ref.basis.model.atoms)
    @test tested.basis.model.atoms[1][2]   == ref.basis.model.atoms[1][2]

    @test tested.basis.Ecut      == ref.basis.Ecut
    @test tested.basis.kweights  == ref.basis.kweights
    @test tested.basis.ksymops   == ref.basis.ksymops
    @test tested.basis.fft_size  == ref.basis.fft_size

    kcoords_test = getproperty.(tested.basis.kpoints, :coordinate)
    kcoords_ref  = getproperty.(ref.basis.kpoints, :coordinate)
    @test kcoords_test == kcoords_ref

    @test tested.n_iter         == ref.n_iter
    @test tested.energies.total ≈  ref.energies.total atol=1e-13
    @test tested.eigenvalues    == ref.eigenvalues
    @test tested.occupation     == ref.occupation
    @test tested.ψ              == ref.ψ
    @test tested.ρ.real         == ref.ρ.real
    @test tested.ρspin.real     == ref.ρspin.real
end

if mpi_nprocs() == 1  # Not yet implemented for MPI-parallel runs
    @testset "Test checkpointing" begin
        O = ElementPsp(o2molecule.atnum, psp=load_psp("hgh/pbe/O-q6.hgh"))
        magnetic_moments = [O => [1., 1.]]
        model = model_PBE(o2molecule.lattice, [O => o2molecule.positions],
                          temperature=0.02, smearing=smearing=Smearing.Gaussian(),
                          magnetic_moments=magnetic_moments)
        basis  = PlaneWaveBasis(model, 4; kgrid=[1, 1, 1])
        ρspin  = guess_spin_density(basis, magnetic_moments)

        # Run SCF and do checkpointing along the way
        scfres = mktempdir() do tmpdir
            checkpointfile = joinpath(tmpdir, "scfres.jld2")
            callback = ScfDefaultCallback() ∘ ScfSaveCheckpoints(checkpointfile; keep=true)
            scfres = self_consistent_field(basis, tol=5e-2, ρspin=ρspin, callback=callback)
            test_scfres_agreement(scfres, load_scfres(checkpointfile))
            scfres
        end

        # Save & load scfres
        mktempdir() do tmpdir
            checkpointfile = joinpath(tmpdir, "scfres.jld2")
            save_scfres(checkpointfile, scfres)
            test_scfres_agreement(scfres, load_scfres(checkpointfile))
        end
    end
end


text:
__precompile__()
using BenchmarkTools, LinearAlgebra

function AddMatRowWise!(MatrixA::Matrix{Float64})
    NumRows,NumCols = size(MatrixA)
    for  i = 1:NumRows
        for j = 1:NumCols
            MatrixA[i,j] = 10*i + j
        end
    end
end

function AddMatColWise!(MatrixA::Matrix{Float64})
    NumRows,NumCols = size(MatrixA)
    for  j = 1:NumCols
        for i = 1:NumRows
            MatrixA[i,j] = 10*i + j
        end
    end
end

function askdimensions()
	print("Enter rows and column for matrix: ");
    Num = parse.(Int,split(readline()))
    return Num[1],Num[2]
end

NumRows, NumCols = askdimensions()

MatrixA = Matrix{Float64}(undef,NumRows,NumCols)

cpu_time_used_row = @belapsed AddMatRowWise!($MatrixA)
cpu_time_used_col = @belapsed AddMatColWise!($MatrixA)

println("Time for allocate row-wise: $cpu_time_used_row seconds")

println("Time for allocate column-wise: $cpu_time_used_col seconds")

max_stars_repo_name: lucasmoschen/tutorial-julia-opt

max_issues_repo_name: lucasmoschen/tutorial-julia-opt

max_forks_repo_name: lucasmoschen/tutorial-julia-opt


text:
"""
`module JAC.TableStrings`
    ... a submodel of JAC that contains various methods to quickly printout the headers of various tables,
        especially if these headers require several lines.
"""
module TableStrings

    using Printf, ..Basics,  ..Defaults

    """
    `TableStrings.center(n::Int64, sa::String)`
        ... a string sa is centered within a string of length n, and optionally, na blanks are appended; a string is returned.
    """
    function  center(n::Int64, sa::String; na::Int64=0)
        sb = sa;   mb = n - length(sa)
        while mb > 2    sb = " "*sb*" ";            mb = mb - 2    end
        if    mb == 2   sb = " "*sb*" "     else    sb = " "*sb    end
        return( sb * repeat(" ", na) )
    end

    """
    `TableStrings.flushleft(n::Int64, sa::String)`
        ... a string sa is left-placed within a string of length n, and optionally, na blanks are appended; a string is returned.
    """
    function  flushleft(n::Int64, sa::String; na::Int64=0)
        sb = sa;   sb = sb * repeat(" ", n+na)
        return( sb[1:n+na] )
    end

    """
    `TableStrings.flushright(n::Int64, sa::String)`
        ... a string sa is right-placed within a string of length n, and optionally, na blanks are appended; a string is returned.
    """
    function  flushright(n::Int64, sa::String; na::Int64=0)
        sb = sa;   sb = repeat(" ", n) * sb;   sb = sb[end-n+1:end]
        return( sb * repeat(" ", na) )
    end

    """
    `TableStrings.hBlank(n::Int64)`  ... a string of n 'blanks' is returned.
    """
    function  hBlank(n::Int64)
        return( repeat(" ", n) )
    end

    """
    `TableStrings.hLine(n::Int64)`  ... a string of n * "-" is returned.
    """
    function  hLine(n::Int64)
        return( repeat("-", n) )
    end

    """
    `TableStrings.inUnits(sa::String)`  ... a string is returned that displays the requested unit in the form: '[unit]'
    """
    function  inUnits(sa::String)
        if       sa == "energy"               return( "[" * Defaults.GBL_ENERGY_UNIT * "]" )
        elseif   sa == "rate"                 return( "[" * Defaults.GBL_RATE_UNIT   * "]" )
        elseif   sa == "time"                 return( "[" * Defaults.GBL_TIME_UNIT   * "]" )
        elseif   sa == "cross section"        return( "[" * Defaults.GBL_CROSS_SECTION_UNIT * "]" )
        elseif   sa == "strength"             return( "[" * Defaults.GBL_STRENGTH_UNIT * "]" )
        elseif   sa == "reduced strength"     return( "[" * Defaults.GBL_STRENGTH_UNIT * "^2" * "]" )
        else     error("stop a")
        end
    end

    """
    `TableStrings.level(i::Int64)`  ... a string "....i" is returned.
    """
    function  level(i::Int64)
        si = string(i);   ni = length(si);    sa = repeat(" ", 5)
        return( sa[1:5-ni] * si )
    end

    """
    `TableStrings.levels_if(i::Int64, f::Int64)`  ... a string "....i -- ....f" is returned
    """
    function  levels_if(i::Int64, f::Int64)
        si = string(i);   ni = length(si);    sf = string(f);   nf = length(sf);    sa = repeat(" ", 5)
        return( sa[1:5-ni] * si * " -- " * sa[1:4-nf] * sf )
    end

    """
    `TableStrings.levels_imf(i::Int64, m::Int64, f::Int64)`  ... a string "....i -- ....m -- ....f" is returned
    """
    function  levels_imf(i::Int64, m::Int64, f::Int64)
        si = string(i);   ni = length(si);    sm = string(m);   nm = length(sm);    sf = string(f);   nf = length(sf);
        sa = repeat(" ", 5)
        return( sa[1:5-ni] * si * " -- " * sa[1:5-nm] * sm * " -- " * sa[1:5-nf] * sf )
    end

    """
    `TableStrings.gaugeMultipolesTupels(n::Int64, mpList::Array{Tuple{Basics.UseGauge,Array{Basics.EmMultipole,1}},1})`
        ... a string of tupels Used gauge[E1,M1,E2,E1, ...] is returned.
    """
    function  gaugeMultipolesTupels(n::Int64, mpList::Array{Tuple{Basics.UseGauge,Array{Basics.EmMultipole,1}},1})
        sa = "";   wa = String[]
        for tt in mpList
           sa = sa * string(tt[1])[1:3] * "["
           for  mp in tt[2]
              sa = sa * string(mp) * ","
           end
           sa = sa[1:end-1] * "], "
           if  length(sa) + 15 > n    push!(wa, sa[1:end-2]);    sa = ""    end
        end
        if  sa != ""    push!(wa, sa[1:end-2])    end
        return( wa )
    end

    """
    `TableStrings.kappaKappaSymmetryTupels(n::Int64, kappaList::Array{Tuple{Int64,Int64,LevelSymmetry},1})`
        ... a list of Strings with maximal length n is returned; each string in this list comprises a number of
            'shell -> symmetry -> shell' descriptors.
    """
    function kappaKappaSymmetryTupels(n::Int64, kappaList::Array{Tuple{Int64,Int64,LevelSymmetry},1})
        sa = "";   wa = String[]
        for k in kappaList
           shin = Subshell(9, k[1]);    shout = Subshell(9, k[2])
           sa = sa * string(shin)[2:end] * "(" * string(k[3]) * ")" * string(shout)[2:end]  * ",  "
           if  length(sa) + 15 > n    push!(wa, sa[1:end-3]);    sa = ""    end
        end
        if  sa != ""    push!(wa, sa[1:end-3])    end
        return( wa )
    end

    """
    `TableStrings.kappaKappaKappaSymmetryTupels(n::Int64, kappaList::Array{Tuple{Int64,Int64,Int64,LevelSymmetry},1})`
        ... a list of Strings with maximal length n is returned; each string in this list comprises a number of
            '(shell -> shell) symmetry [shell]' descriptors.
    """
    function kappaKappaKappaSymmetryTupels(n::Int64, kappaList::Array{Tuple{Int64,Int64,Int64,LevelSymmetry},1})
        sa = "";   wa = String[]
        for k in kappaList
           shin = Subshell(9, k[1]);    shout = Subshell(9, k[2]);    shauger = Subshell(9, k[3])
           sa = sa * "(" * string(shin)[2:end] * "->" * string(shout)[2:end] * ") " * string(k[4]) * "[" * string(shauger)[2:end]  * "],  "
           if  length(sa) + 18 > n    push!(wa, sa[1:end-3]);    sa = ""    end
        end
        if  sa != ""    push!(wa, sa[1:end-3])    end
        return( wa )
    end

    """
    `TableStrings.kappaMultipoleSymmetryTupels(n::Int64, kappaList::Array{Tuple{Int64,EmMultipole,EmGauge,LevelSymmetry},1})`
        ... a list of Strings with maximal length n is returned; each string in this list comprises a number of
            'shell (multipole, gauge, symmetry)' descriptors.
    """
    function  kappaMultipoleSymmetryTupels(n::Int64, kappaList::Array{Tuple{Int64,EmMultipole,EmGauge,LevelSymmetry},1})
        sa = "";   wa = String[]
        for k in kappaList
           sh = Subshell(9, k[1])
           sa = sa * string(sh)[2:end] * " (" * string(k[2]) * ", " * string(k[3])[1:3] * "; " * string(k[4])* "),  "
           if  length(sa) + 15 > n    push!(wa, sa[1:end-3]);    sa = ""    end
        end
        if  sa != ""    push!(wa, sa[1:end-3])    end
        return( wa )
    end

    """
    `TableStrings.processSymmetryEnergyTupels(n::Int64, pList::Array{Tuple{Basics.AbstractProcess, Int64, LevelSymmetry, Float64},1}, sc::String)`
         ... a list of Strings with maximal length n is returned; each string in this list comprises a number of
            'sc(process: NoElectrons, symmetry, energy),  ' descriptors.
    """
    function  processSymmetryEnergyTupels(n::Int64, pList::Array{Tuple{Basics.AbstractProcess, Int64, LevelSymmetry, Float64},1}, sc::String)
        sa = "";   wa = String[]
        for p in pList
           sa = sa * sc * "[" * string(p[1])[1:1] * ": " * string(p[2]) * ", " * string(p[3]) * ", "
           sa = sa * @sprintf("%.4e", p[4]) * "],  "
           if  length(sa) + 20 > n    push!(wa, sa[1:end-3]);    sa = ""    end
        end
        if  sa != ""    push!(wa, sa[1:end-3])    end
        return( wa )
    end

    """
    `TableStrings.kappaSymmetryTupels(n::Int64, kappaList::Array{Tuple{Int64,LevelSymmetry},1})`
        ... a string of shell (Symmetry) is returned.
    """
    function  kappaSymmetryTupels(n::Int64, kappaList::Array{Tuple{Int64,LevelSymmetry},1})
        sa = ""
        for k in kappaList
           sh = Subshell(9, k[1])
           sa = sa * string(sh)[2:end] * " (" * string(k[2]) * "),  "
        end
        return( sa[1:end-3] )
    end

    """
    `TableStrings.multipoleGaugeTupels(n::Int64, mpList::Array{Tuple{Basics.EmMultipole,Basics.EmGauge},1})`
        ... a string of tupels E1(Gauge), M1(Gauge) is returned.
    """
    function  multipoleGaugeTupels(n::Int64, mpList::Array{Tuple{Basics.EmMultipole,Basics.EmGauge},1})
        sa = ""
        for mp in mpList
           sa = sa * string(mp[1]) * "(" * string(mp[2]) * "), "
        end
        return( sa[1:end-2] )
    end

    """
    `TableStrings.multipoleList(mpList::Array{EmMultipole,1})`  ... a string of multipoles is returned.
    """
    function  multipoleList(mpList::Array{EmMultipole,1})
        sa = ""
        for mp in mpList
           sa = sa * string(mp) * ", "
        end
        return( sa[1:end-2] )
    end

    """
    `TableStrings.symmetries_if()`  ... a string 'J^P --> J^P' is returned
    """
    function  symmetries_if(isym::LevelSymmetry,  fsym::LevelSymmetry)
        sa = string(isym) * " --> " * string(fsym)
        return( sa )
    end

    """
    `TableStrings.symmetries_imf()`  ... a string 'J^P --> J^P --> J^P' is returned
    """
    function  symmetries_imf(isym::LevelSymmetry, msym::LevelSymmetry, fsym::LevelSymmetry)
        sa = string(isym) * " --> " * string(msym) * " --> " * string(fsym)
        return( sa )
    end

    """
    `TableStrings.twoPhotonGaugeTupels(n::Int64, mpList::Array{Tuple{AngularJ64, Basics.EmMultipole, LevelSymmetry, Basics.EmMultipole, Basics.EmGauge},1})
        ... a string of tupels (K-rank, multipole_1, Jsym, multipole_2, gauge) is returned.
    """
    function  twoPhotonGaugeTupels(n::Int64, mpList::Array{Tuple{AngularJ64, Basics.EmMultipole, LevelSymmetry, Basics.EmMultipole, Basics.EmGauge},1})
        sa = "";   wa = String[];   mpReduced = String[]
        for tt in mpList
           sb = "(" * string(tt[1]) * ", " * string(tt[2]) * ", " * string(tt[3]) * ", " * string(tt[4]) * ", " * string(tt[5])[1:3] * ")"
           if    sb in mpReduced    else    push!(mpReduced, sb)    end
        end
        for sb in mpReduced
           sa = sa * sb * ", "
           if  length(sa) + 20 > n    push!(wa, sa[1:end-2]);    sa = ""    end
        end
        if  sa != ""    push!(wa, sa[1:end-2])    end
        return( wa )
    end

end # module


text:
module Model

using Turing
using ..MangroveBase: ReviewsSummary, Reviews, kids, subs, normalize
import ..MangroveBase.subs, ..MangroveBase.kids

"""
Construct a `Beta` distribution with given `mean` and spikiness given by `certainty`.
`certainty` is uninformative at 2.
"""
MeanBeta(mean::Real, certainty::Real = 3) = Beta(certainty * mean, certainty * (1 - mean))

# Generative Bayesian model representing the Mangrove Reviews data.
@model mangrove_model(data, typical_rating) = begin
    # Assume that each subject has a fundamental quality associated with it.
    # Extreme qualities are not as likely as moderate ones.
    qualities = Dict()
    for sub in subs(data)
        qualities[sub] ~ MeanBeta(typical_rating)
    end
    # Assume that each reviewer is either neutral (0) or biased (1).
    biases = Dict()
    for kid in kids(data)
        biases[kid] ~ Bernoulli(0.3)
    end
    # Observe the ratings from the dataset, which depend on quality and bias.
    for info in keys(data)
      if biases[info.kid] == 0
            # Neutral reviewer leaves a review which reflects the quality.
            data[info][1] ~ MeanBeta(qualities[info.sub])
            # Neutral reviewer is fairly likely to leave some text opinion.
            data[info][2] ~ Bernoulli(0.8)
        else
            # Biased reviewers are mostly negative and pay not attention to actual quality.
            data[info][1] ~ Beta(0.04, 0.08)
            # Biased reviewers are less likely to leave a text opinion.
            data[info][2] ~ Bernoulli(0.1)
        end
    end
end

"Normalize the data and run inference using the model."
function get_chains(raw_data::Reviews)::Chains
    data = convert(ReviewsSummary, raw_data)
    typical_rating = mean(values(data))
    sampler = PG(1000)
    sample(mangrove_model(data, typical_rating), sampler, 100)
end

# Extract useful data from the resulting MarkovChains: `mean(chn::Chains)`

subs(qualities_mean::ChainDataFrame)::Vector{String} =
  [match(r"qualities\[(.+)\]", string(p))[1] for p in qualities_mean[:, :parameters]]
qualities(qualities_mean::ChainDataFrame)::Vector{Int16} =
  [round(Int16, m*100) for m in qualities_mean[:, :mean]]

kids(biases_mean::ChainDataFrame)::Vector{String} =
  [match(r"biases\[(.+)\]", string(p))[1] for p in biases_mean[:, :parameters]]
neutralities(biases_mean::ChainDataFrame)::Vector{Float32} =
  [1 - m for m in biases_mean[:, :mean]]

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
module DuckDB
using DataFrames, Dates, DuckDB_jll
using DBInterface, Tables

export DBInterface

include("api.jl")
include("consts.jl")

const DBHandle = Ref{Ptr{Cvoid}}    # DuckDB DB connection handle (DB.Interface)

"""
    toDataFrame(connection::Ref{Ptr{Cvoid}},query::String)::DataFrame

Creates a DataFrame from a SQL query within a connection.

* `connection`: The connection to perform the query in.
* `query`: The SQL query to run.
* returns: the abstract dataframe

"""
function toDataFrame(connection::Ref{Ptr{Cvoid}}, query::String)::DataFrame
    res = execute(connection, query)::Ref{duckdb_result}
    return toDataFrame(res)
end
"""
    toDataFrame(result::Ref{duckdb_result})::DataFrame

Creates a DataFrame from the full result
* `result`: the full result from `execute`
* returns: the abstract dataframe

"""
function toDataFrame(result::Ref{duckdb_result})::DataFrame
    columns = unsafe_wrap(Array{duckdb_column}, result[].columns, Int64(result[].column_count))
    df = DataFrame()
    for i = 1:Int64(result[].column_count)
        rows = Int64(result[].row_count)
        name = (unsafe_string(columns[i].name))
        type = DUCKDB_TYPE(Int64(columns[i].type))
        if type == DUCKDB_TYPE_INVALID
            print("invalid type for column - \"" * name * "\"")
        else
            mask = unsafe_wrap(Array, columns[i].nullmask, rows)
            data = unsafe_wrap(Array, Ptr{DUCKDB_TYPES[type]}(columns[i].data), rows)
            bmask = reinterpret(Bool, mask)

            if 0 != sum(mask)
                data = data[.!bmask]
            end

            # Format the different datatypes for DataFrame
            if type == DUCKDB_TYPE_BOOLEAN
                data = Bool.(data)
            elseif type == DUCKDB_TYPE_DATE
                data = Dates.epochdays2date.(data .+ 719528)
            elseif type == DUCKDB_TYPE_TIME
                data = Dates.Time.(Dates.Nanosecond.(data .* 1000))
            elseif type == DUCKDB_TYPE_TIMESTAMP
                data = Dates.epochms2datetime.((data ./ 1000) .+ 62167219200000)
            elseif type == DUCKDB_TYPE_INTERVAL
                data = map(
                    x -> Dates.CompoundPeriod(
                        Dates.Month(x.months),
                        Dates.Day(x.days),
                        Dates.Microsecond(x.micros),
                    ),
                    data,
                )
            elseif type == DUCKDB_TYPE_HUGEINT
                data = map(x -> x.upper < 1 ? (x.lower::UInt64)%Int64 : x, data)
            elseif type == DUCKDB_TYPE_VARCHAR
                data = unsafe_string.(data)
            end

            if 0 != sum(mask)
                fulldata = Array{Union{Missing,eltype(data)}}(missing, rows)
                fulldata[.!bmask] = data
                data = fulldata
            end

            df[!, name] = data
        end
    end
    return df
end

"""
	open(path, out_database)

Creates a new database or opens an existing database file stored at the the given path.
If no path is given a new in-memory database is created instead.

* `path`: Path to the database file on disk, or `nullptr` or `:memory:` to open an in-memory database.
* returns: `database`: The result database object.
"""
function open(path::String)::Ref{Ptr{Cvoid}}
    database = Ref{Ptr{Cvoid}}()
    duckdb_open(path,database)
    return database
end

"""
    close(database)
Closes the specified database and de-allocates all memory allocated for that database.\n
This should be called after you are done with any database allocated through duckdb_open.\n
Note that failing to call duckdb_close (in case of e.g. a program crash) will not cause data corruption. Still it is recommended to always correctly close a database object after you are done with it.
*`database`: the database object to shut down.

"""
function close(database)
    duckdb_close(database)
end

open() = open(":memory:")

"""
    connect(path)
Creates a new database or opens an existing database file stored at the the given path. If no path is given a new in-memory database is created instead.
* `path`: Path to the database file on disk or `:memory:` to open an in-memory database.
* returns: a connection handle

"""
function connect(database)::Ref{Ptr{Cvoid}}
    connection = Ref{Ptr{Cvoid}}()
    duckdb_connect(database, connection)
    return connection
end

"""
    disconnect(connection)
Closes the specified connection and de-allocates all memory allocated for that connection.
* `connection`: The connection to close.

"""
function disconnect(connection)
    duckdb_disconnect(connection)
end

"""
    execute(connection, query)

Executes a SQL query within a connection and returns the full (materialized) result. If the query fails to execute, `DuckDBError` is returned and the error message can be retrieved by calling `duckdb_result_error`.

Note that after running duckdb_query, duckdb_destroy_result must be called on the result object even if the query fails, otherwise the error stored within the result will not be freed correctly.
* `connection`: The connection to perform the query in.
* `query`: The SQL query to run.
* returns: the full result pointer

"""
function execute(connection::Ref{Ptr{Cvoid}}, query::String)::Ref{duckdb_result}
    result = Ref{duckdb_result}()
    duckdb_query(connection, query, result)
    if result[].error_message != Ptr{UInt8}(0)
        print(unsafe_string(result[].error_message))
    end
    return result
end

struct DuckDBException <: Exception
    msg::AbstractString
end

duckdb_errmsg(handle) = unsafe_string(handle[].error_message)
duckdbexception(handle::DBHandle) = DuckDBException(duckdb_errmsg(handle))
duckdberror(handle::DBHandle) = throw(duckdbexception(handle))

"""
    appender_create(connection, table, schema)

Creates appender object used to append new entries to an existing table

* `connection`: The connection to create the appender for
* `table`: Name of the table to append to
* `schema`: Name of the schema of the table to append to, optional argument that defaults to "main"
* returns: The appender object
"""
function appender_create(connection::Ref{Ptr{Cvoid}}, table::String, schema::String="main")
    appender = Ref{Ptr{Cvoid}}()
    duckdb_appender_create(connection, schema, table, appender)
    return appender
end

"""
    appendDataFrame(df, connection, table, schema)

Appends DataFrame df to a table in database.
Note that the datatype and order of the columns have to be consistent with the layout of the table.
The column names of the dataframe are irrelevant as the data in the dataframe will be stored with column names of the table in the database.
The datatypes currently supported by the appender are: Int, Float32, Float64, String, Dates.Time, Dates.Date, Missing

* `df`: Dataframe that will be appended to the table
* `connection`: Connection to the database the dataframe will be appended to
* `table`: Name of the table the data is to be appended to
* `schema`: Name of the schema of the table to append to, optional argument that defaults to "main"
"""
function appendDataFrame(df::DataFrame, connection::Ref{Ptr{Cvoid}}, table::String, schema::String="main")
    # Determine the type of each column
    dftypes = eltype.(eachcol(df))

    # Create a vector holding the right appender functions
    appender_function_vector= Vector{Function}(undef, ncol(df))
    for (pos, dtype) in enumerate(dftypes)
        appender_function_vector[pos] = _appender_function(dtype)
    end

    # Create the appender
    appender = appender_create(connection, table, schema)
    # Append the data to the table
    for row in eachrow(df)

        # For each value in the row apply the right appender function
        for (index, element) in enumerate(row)
            appender_function_vector[index](appender, element)
        end

        # End the appender row
        duckdb_appender_end_row(appender)
    end

    # After all rows have been added to the appender flush them to the table and destroy the appender
    duckdb_appender_destroy(appender)
end

"""
    _appender_function(dtype)

Returns the correct appender function for the respective data type to append the data to a table
* `dtpye`: Datatype of data to be appended
"""
function _appender_function(dtype::DataType)
    if dtype == Bool
        return duckdb_append_bool
    elseif dtype == Int8
        return duckdb_append_int8
    elseif dtype == Int16
        return duckdb_append_int16
    elseif dtype == Int32
        return duckdb_append_int32
    elseif dtype == Int64
        return duckdb_append_int64
    elseif dtype == UInt8
        return duckdb_append_uint8
    elseif dtype == UInt16
        return duckdb_append_uint16
    elseif dtype == UInt32
        return duckdb_append_uint32
    elseif dtype == UInt64
        return duckdb_append_uint64
    elseif dtype == Float32
        return duckdb_append_float
    elseif dtype == Float64
        return duckdb_append_double
    elseif dtype == Dates.Date
        return duckdb_append_date
    elseif dtype == Dates.Time
        return duckdb_append_time
    elseif dtype == String
        return duckdb_append_varchar
    else
        throw("Unsupported datatype $dtype datatype encountered. Supported types are: Int, Float32, Float64, String, Dates.Time, Dates.Date, Missing")
    end
end

_appender_function(dtype::Union) = return _union_appender

"""
    _union_appender

Appender function that can be used to append dataframe columns with union datatypes Union{T, Missing} to the database table
* `appender`: Appender object used to append the values to the table
* `value`: Value that is to be appended to the table
"""
function _union_appender(appender, value)
    value_type = typeof(value)
    if value_type == Missing
        return duckdb_append_null(appender)
    else
        return _appender_function(value_type)(appender, value)
    end
end

"""
    duckdb_append_date(appender, value)

Appender function that converts the Dates.Date format back to the stored date format in DuckDB and appends it to the table
* `appender`: Appender object used to append the values to the table
* `value`: Value that is to be appended to the table
"""
function duckdb_append_date(appender::Ref{Ptr{Cvoid}}, value::Dates.Date)
    converted_value = Dates.date2epochdays(value) - 719528
    duckdb_append_date(appender, converted_value)
end

"""
    duckdb_append_time(appender, value)

Appender function that converts the Dates.Time format back to the stored time format in DuckDB and appends it to the table
* `appender`: Appender object used to append the values to the table
* `value`: Value that is to be appended to the table
"""
function duckdb_append_time(appender::Ref{Ptr{Cvoid}}, value::Dates.Time)
    converted_value = hour(value) * 60 * 60 * 1000 * 1000 + minute(value) * 60 * 1000 * 1000 + second(value) * 1000 * 1000 + millisecond(value) * 1000 + microsecond(value)
    duckdb_append_time(appender, converted_value)
end

"""
    `DuckDB.DB()` => in-memory DuckDB database
    `DuckDB.DB(file)` => file-based DuckDB database

Constructors for a representation of a DuckDB database, either backed by an on-disk file or in-memory.

`DuckDB.DB` requires the `file` string argument in the 2nd definition
as the name of either a pre-defined DuckDB database to be opened,
or if the file doesn't exist, a database will be created.

The `DuckDB.DB` object represents a single connection to a DuckDB database.
All other DuckDB.jl functions take an `DuckDB.DB` as the first argument as context.

To create an in-memory temporary database, call `DuckDB.DB()`.

The `DuckDB.DB` will be automatically closed/shutdown when it goes out of scope
(i.e. the end of the Julia session, end of a function call wherein it was created, etc.)

NOTE: This borrows heavily from SQLite.jl, [here](https://github.com/JuliaDatabases/SQLite.jl/blob/9724a175ae30b22cb01775250ce2e87317706215/src/SQLite.jl)
"""
mutable struct DB <: DBInterface.Connection
    file::String
    handle::DBHandle

    # TODO: Support prepared statements

    function DB(f::AbstractString)
        f = String(isempty(f) ? f : expanduser(f))

        try
            dbhandle = open(f)
            handle = connect(dbhandle)

            db = new(f, handle)
            finalizer(_close, db)
            return db
        catch
            duckdberror(handle)
        end
    end
end

function _close(db::DB)
    db.handle == C_NULL || duckdb_disconnect(db.handle)
    db.handle = C_NULL
    return
end

DB() = DB(":memory:")
DBInterface.close!(db::DB) = _close(db)
Base.close(db::DB) = _close(db)

Base.show(io::IO, db::DB) = print(io, string("DB(", "\"$(db.file)\"", ")"))

"""
    DBInterface.execute(db::DB, sql::String, [params])

TODO: Support DBInterface statements

Take inputs given by `db` and `sql` execute the query and return an iterator of result rows.

Note that the returned result row iterator only supports a single-pass, forward-only iteration of the result rows.

TODO: Support Tables.jl
"""
function DBInterface.execute(db::DB, sql::String)

    result = execute(db.handle, sql)

    return result
end

end # module

max_issues_repo_name: kimmolinna/julia-duckdb


text:
using OptimNewton
using Test

@testset "OptimNewton.jl" begin
    # get_root.
    # x0 = pi/10, f(x) = sin(x), fd(x) = cos(x).
    @test isapprox(get_root(pi/10, x -> sin(x), x -> cos(x)), 0.0);
    # x0 = 1, f(x) = 2 - x^2 + x^3, fd(x) = -2x + 3x^2.
    @test isapprox(get_root(1, x -> 2 - x^2 + x^3, x -> -2x + 3x^2),-1);

    # optimize
    # x0 = pi/3, f(x) = sin(x), fd(x) = cos(x), sd(x) = -sin(x)
    @test isapprox(optimize(pi/3, x -> sin(x), x -> cos(x), x -> -sin(x)), pi/2);
    # x0 = 1, f(x) = 2 - x^2 + x^3, fd(x) = -2x + 3x^2, sd(x) = -2 + 6x.
    @test isapprox(optimize(1, x -> 2 - x^2 + x^3, x -> -2x + 3x^2, x -> -2 + 6x),2/3);
 end