text:
##############################################################################
#
# Text output
#
##############################################################################

function escapedprint(io::IO, x::Any, escapes::AbstractString)
    ourshowcompact(io, x)
end

function escapedprint(io::IO, x::AbstractString, escapes::AbstractString)
    escape_string(io, x, escapes)
end

function printtable(io::IO,
                    dt::AbstractDataTable;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    nastring::AbstractString = "NULL")
    n, p = size(dt)
    etypes = eltypes(dt)
    if header
        cnames = _names(dt)
        for j in 1:p
            print(io, quotemark)
            print(io, cnames[j])
            print(io, quotemark)
            if j < p
                print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    quotestr = string(quotemark)
    for i in 1:n
        for j in 1:p
            if !isnull(dt[j][i])
                if ! (etypes[j] <: Real)
                    print(io, quotemark)
                    x = unsafe_get(dt[i, j])
                    escapedprint(io, x, quotestr)
                    print(io, quotemark)
                else
                    print(io, dt[i, j])
                end
            else
                print(io, nastring)
            end
            if j < p
                print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    return
end

function printtable(dt::AbstractDataTable;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    nastring::AbstractString = "NULL")
    printtable(STDOUT,
               dt,
               header = header,
               separator = separator,
               quotemark = quotemark,
               nastring = nastring)
    return
end
##############################################################################
#
# HTML output
#
##############################################################################

function html_escape(cell::AbstractString)
    cell = replace(cell, "&", "&")
    cell = replace(cell, "<", "<")
    cell = replace(cell, ">", ">")
    return cell
end

@compat function Base.show(io::IO, ::MIME"text/html", dt::AbstractDataTable)
    cnames = _names(dt)
    write(io, "
")
    write(io, "
")
    write(io, "	")
    for column_name in cnames
        write(io, "$column_name	")
    end
    write(io, "
")
    haslimit = get(io, :limit, true)
    n = size(dt, 1)
    if haslimit
        tty_rows, tty_cols = _displaysize(io)
        mxrow = min(n,tty_rows)
    else
        mxrow = n
    end
    for row in 1:mxrow
        write(io, "
")
        write(io, "$row	")
        for column_name in cnames
            cell = sprint(ourshowcompact, dt[row, column_name])
            write(io, "$(html_escape(cell))	")
        end
        write(io, "
")
    end
    if n > mxrow
        write(io, "
")
        write(io, "⋮	")
        for column_name in cnames
            write(io, "⋮	")
        end
        write(io, "
")
    end
    write(io, "
")
end

##############################################################################
#
# LaTeX output
#
##############################################################################

function latex_char_escape(char::AbstractString)
    if char == "\\"
        return "\\textbackslash{}"
    elseif char == "~"
        return "\\textasciitilde{}"
    else
        return string("\\", char)
    end
end

function latex_escape(cell::AbstractString)
    cell = replace(cell, ['\\','~','#','$','%','&','_','^','{','}'], latex_char_escape)
    return cell
end

function Base.show(io::IO, ::MIME"text/latex", dt::AbstractDataTable)
    nrows = size(dt, 1)
    ncols = size(dt, 2)
    cnames = _names(dt)
    alignment = repeat("c", ncols)
    write(io, "\\begin{tabular}{r|")
    write(io, alignment)
    write(io, "}\n")
    write(io, "\t& ")
    header = join(map(c -> latex_escape(string(c)), cnames), " & ")
    write(io, header)
    write(io, "\\\\\n")
    write(io, "\t\\hline\n")
    for row in 1:nrows
        write(io, "\t")
        write(io, @sprintf("%d", row))
        for col in 1:ncols
            write(io, " & ")
            cell = dt[row,col]
            if !isnull(cell)
                content = unsafe_get(cell)
                if mimewritable(MIME("text/latex"), content)
                    show(io, MIME("text/latex"), content)
                else
                    print(io, latex_escape(sprint(ourshowcompact, content)))
                end
            end
        end
        write(io, " \\\\\n")
    end
    write(io, "\\end{tabular}\n")
end

##############################################################################
#
# MIME
#
##############################################################################

@compat function Base.show(io::IO, ::MIME"text/csv", dt::AbstractDataTable)
    printtable(io, dt, true, ',')
end

@compat function Base.show(io::IO, ::MIME"text/tab-separated-values", dt::AbstractDataTable)
    printtable(io, dt, true, '\t')
end

##############################################################################
#
# DataStreams-based IO
#
##############################################################################

importall DataStreams
using WeakRefStrings

# DataTables DataStreams implementation
function Data.schema(df::DataTable, ::Type{Data.Column})
    return Data.Schema(map(string, names(df)),
                       DataType[typeof(A) for A in df.columns], size(df, 1))
end

# DataTable as a Data.Source
function Data.isdone(source::DataTable, row, col)
    rows, cols = size(source)
    return row > rows || col > cols
end

Data.streamtype(::Type{DataTable}, ::Type{Data.Column}) = true
Data.streamtype(::Type{DataTable}, ::Type{Data.Field}) = true

Data.streamfrom{T <: AbstractVector}(source::DataTable, ::Type{Data.Column}, ::Type{T}, col) =
    (@inbounds A = source.columns[col]::T; return A)
Data.streamfrom{T}(source::DataTable, ::Type{Data.Column}, ::Type{T}, col) =
    (@inbounds A = source.columns[col]; return A)
Data.streamfrom{T}(source::DataTable, ::Type{Data.Field}, ::Type{T}, row, col) =
    (@inbounds A = Data.streamfrom(source, Data.Column, T, col); return A[row]::T)

# DataTable as a Data.Sink
allocate{T}(::Type{T}, rows, ref) = Array{T}(rows)
allocate{T}(::Type{Vector{T}}, rows, ref) = Array{T}(rows)

allocate{T}(::Type{Nullable{T}}, rows, ref) =
    NullableArray{T, 1}(Array{T}(rows), fill(true, rows), isempty(ref) ? UInt8[] : ref)
allocate{T}(::Type{NullableVector{T}}, rows, ref) =
    NullableArray{T, 1}(Array{T}(rows), fill(true, rows), isempty(ref) ? UInt8[] : ref)

allocate{S,R}(::Type{CategoricalArrays.CategoricalValue{S,R}}, rows, ref) =
    CategoricalArray{S,1,R}(rows)
allocate{S,R}(::Type{CategoricalVector{S,R}}, rows, ref) =
    CategoricalArray{S,1,R}(rows)

allocate{S,R}(::Type{Nullable{CategoricalArrays.CategoricalValue{S,R}}}, rows, ref) =
    NullableCategoricalArray{S,1,R}(rows)
allocate{S,R}(::Type{NullableCategoricalVector{S,R}}, rows, ref) =
    NullableCategoricalArray{S,1,R}(rows)

function DataTable{T <: Data.StreamType}(sch::Data.Schema,
                                         ::Type{T}=Data.Field,
                                         append::Bool=false,
                                         ref::Vector{UInt8}=UInt8[], args...)
    rows, cols = size(sch)
    rows = max(0, T <: Data.Column ? 0 : rows) # don't pre-allocate for Column streaming
    columns = Vector{Any}(cols)
    types = Data.types(sch)
    for i = 1:cols
        columns[i] = allocate(types[i], rows, ref)
    end
    return DataTable(columns, map(Symbol, Data.header(sch)))
end

# given an existing DataTable (`sink`), make any necessary changes for streaming source
# with Data.Schema `sch` to it, given we know if we'll be `appending` or not
function DataTable(sink, sch::Data.Schema, ::Type{Data.Field}, append::Bool,
                   ref::Vector{UInt8})
    rows, cols = size(sch)
    newsize = max(0, rows) + (append ? size(sink, 1) : 0)
    # need to make sure we don't break a NullableVector{WeakRefString{UInt8}} when appending
    if append
        for (i, T) in enumerate(Data.types(sch))
            if T <: Nullable{WeakRefString{UInt8}}
                sink.columns[i] = NullableArray(String[string(get(x, "")) for x in sink.columns[i]])
                sch.types[i] = Nullable{String}
            end
        end
    end
    newsize != size(sink, 1) && foreach(x->resize!(x, newsize), sink.columns)
    sch.rows = newsize
    return sink
end
function DataTable(sink, sch::Data.Schema, ::Type{Data.Column}, append::Bool, ref::Vector{UInt8})
    rows, cols = size(sch)
    append ? (sch.rows += size(sink, 1)) : foreach(empty!, sink.columns)
    return sink
end

Data.streamtypes(::Type{DataTable}) = [Data.Column, Data.Field]

Data.streamto!{T}(sink::DataTable, ::Type{Data.Field}, val::T, row, col, sch::Data.Schema{false}) =
    push!(sink.columns[col]::Vector{T}, val)
Data.streamto!{T}(sink::DataTable, ::Type{Data.Field}, val::Nullable{T}, row, col, sch::Data.Schema{false}) =
    push!(sink.columns[col]::NullableVector{T}, val)
Data.streamto!{T, R}(sink::DataTable, ::Type{Data.Field}, val::CategoricalValue{T, R}, row, col, sch::Data.Schema{false}) =
    push!(sink.columns[col]::CategoricalVector{T, R}, val)
Data.streamto!{T, R}(sink::DataTable, ::Type{Data.Field}, val::Nullable{CategoricalValue{T, R}}, row, col, sch::Data.Schema{false}) =
    push!(sink.columns[col]::NullableCategoricalVector{T, R}, val)
Data.streamto!{T}(sink::DataTable, ::Type{Data.Field}, val::T, row, col, sch::Data.Schema{true}) =
    (sink.columns[col]::Vector{T})[row] = val
Data.streamto!{T}(sink::DataTable, ::Type{Data.Field}, val::Nullable{T}, row, col, sch::Data.Schema{true}) =
    (sink.columns[col]::NullableVector{T})[row] = val
Data.streamto!(sink::DataTable, ::Type{Data.Field}, val::Nullable{WeakRefString{UInt8}}, row, col, sch::Data.Schema{true}) =
    sink.columns[col][row] = val
Data.streamto!{T, R}(sink::DataTable, ::Type{Data.Field}, val::CategoricalValue{T, R}, row, col, sch::Data.Schema{true}) =
    (sink.columns[col]::CategoricalVector{T, R})[row] = val
Data.streamto!{T, R}(sink::DataTable, ::Type{Data.Field}, val::Nullable{CategoricalValue{T, R}}, row, col, sch::Data.Schema{true}) =
    (sink.columns[col]::NullableCategoricalVector{T, R})[row] = val

function Data.streamto!{T}(sink::DataTable, ::Type{Data.Column}, column::T, row, col, sch::Data.Schema)
    if row == 0
        sink.columns[col] = column
    else
        append!(sink.columns[col]::T, column)
    end
    return length(column)
end

function Base.append!{T}(dest::NullableVector{WeakRefString{T}}, column::NullableVector{WeakRefString{T}})
    offset = length(dest.values)
    parentoffset = length(dest.parent)
    append!(dest.isnull, column.isnull)
    append!(dest.parent, column.parent)
    # appending new data to `dest` would invalid all existing WeakRefString pointers
    resize!(dest.values, length(dest) + length(column))
    for i = 1:offset
        old = dest.values[i]
        dest.values[i] = WeakRefString{T}(pointer(dest.parent, old.ind), old.len, old.ind)
    end
    for i = 1:length(column)
        old = column.values[i]
        dest.values[offset + i] = WeakRefString{T}(pointer(dest.parent, parentoffset + old.ind), old.len, parentoffset + old.ind)
    end
    return length(dest)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule dSFMT_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("dSFMT")
JLLWrappers.@generate_main_file("dSFMT", UUID("05ff407c-b0c1-5878-9df8-858cc2e60c36"))
end  # module dSFMT_jll

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Twice-precision arithmetic.

# Necessary for creating nicely-behaved ranges like r = 0.1:0.1:0.3
# that return r[3] == 0.3.  Otherwise, we have roundoff error due to
#     0.1 + 2*0.1 = 0.30000000000000004

"""
    hi, lo = splitprec(F::Type{<:AbstractFloat}, i::Integer)

Represent an integer `i` as a pair of floating-point numbers `hi` and
`lo` (of type `F`) such that:
- `widen(hi) + widen(lo) ≈ i`. It is exact if 1.5 * (number of precision bits in `F`) is greater than the number of bits in `i`.
- all bits in `hi` are more significant than any of the bits in `lo`
- `hi` can be exactly multiplied by the `hi` component of another call to `splitprec`.

In particular, while `convert(Float64, i)` can be lossy since Float64
has only 53 bits of precision, `splitprec(Float64, i)` is exact for
any Int64/UInt64.
"""
function splitprec(::Type{F}, i::Integer) where {F<:AbstractFloat}
    hi = truncbits(F(i), cld(precision(F), 2))
    ihi = oftype(i, hi)
    hi, F(i - ihi)
end

function truncmask(x::F, mask) where {F<:IEEEFloat}
    reinterpret(F, mask & reinterpret(uinttype(F), x))
end
truncmask(x, mask) = x

function truncbits(x::F, nb) where {F<:IEEEFloat}
    truncmask(x, typemax(uinttype(F)) << nb)
end
truncbits(x, nb) = x

## Dekker arithmetic

"""
    hi, lo = canonicalize2(big, little)

Generate a representation where all the nonzero bits in `hi` are more
significant than any of the nonzero bits in `lo`. `big` must be larger
in absolute value than `little`.
"""
function canonicalize2(big, little)
    h = big+little
    h, (big - h) + little
end

"""
    zhi, zlo = add12(x, y)

A high-precision representation of `x + y` for floating-point
numbers. Mathematically, `zhi + zlo = x + y`, where `zhi` contains the
most significant bits and `zlo` the least significant.

Because of the way floating-point numbers are printed, `lo` may not
look the way you might expect from the standpoint of decimal
representation, even though it is exact from the standpoint of binary
representation.

Example:
```julia
julia> 1.0 + 1.0001e-15
1.000000000000001

julia> big(1.0) + big(1.0001e-15)
1.000000000000001000100000000000020165767380775934141445417482375879192346701529

julia> hi, lo = Base.add12(1.0, 1.0001e-15)
(1.000000000000001, -1.1012302462515652e-16)

julia> big(hi) + big(lo)
1.000000000000001000100000000000020165767380775934141445417482375879192346701529
```

`lo` differs from 1.0e-19 because `hi` is not exactly equal to
the first 16 decimal digits of the answer.
"""
function add12(x::T, y::T) where {T}
    x, y = ifelse(abs(y) > abs(x), (y, x), (x, y))
    canonicalize2(x, y)
end
add12(x, y) = add12(promote(x, y)...)

"""
    zhi, zlo = mul12(x, y)

A high-precision representation of `x * y` for floating-point
numbers. Mathematically, `zhi + zlo = x * y`, where `zhi` contains the
most significant bits and `zlo` the least significant.

Example:
```julia
julia> x = Float32(π)
3.1415927f0

julia> x * x
9.869605f0

julia> Float64(x) * Float64(x)
9.869604950382893

julia> hi, lo = Base.mul12(x, x)
(9.869605f0, -1.140092f-7)

julia> Float64(hi) + Float64(lo)
9.869604950382893
```
"""
function mul12(x::T, y::T) where {T<:AbstractFloat}
    h = x * y
    ifelse(iszero(h) | !isfinite(h), (h, h), canonicalize2(h, fma(x, y, -h)))
end
mul12(x::T, y::T) where {T} = (p = x * y; (p, zero(p)))
mul12(x, y) = mul12(promote(x, y)...)

"""
    zhi, zlo = div12(x, y)

A high-precision representation of `x / y` for floating-point
numbers. Mathematically, `zhi + zlo ≈ x / y`, where `zhi` contains the
most significant bits and `zlo` the least significant.

Example:
```julia
julia> x, y = Float32(π), 3.1f0
(3.1415927f0, 3.1f0)

julia> x / y
1.013417f0

julia> Float64(x) / Float64(y)
1.0134170444063078

julia> hi, lo = Base.div12(x, y)
(1.013417f0, 3.8867366f-8)

julia> Float64(hi) + Float64(lo)
1.0134170444063066
"""
function div12(x::T, y::T) where {T<:AbstractFloat}
    # We lose precision if any intermediate calculation results in a subnormal.
    # To prevent this from happening, standardize the values.
    xs, xe = frexp(x)
    ys, ye = frexp(y)
    r = xs / ys
    rh, rl = canonicalize2(r, -fma(r, ys, -xs)/ys)
    ifelse(iszero(r) | !isfinite(r), (r, r), (ldexp(rh, xe-ye), ldexp(rl, xe-ye)))
end
div12(x::T, y::T) where {T} = (p = x / y; (p, zero(p)))
div12(x, y) = div12(promote(x, y)...)

## TwicePrecision

"""
    TwicePrecision{T}(hi::T, lo::T)
    TwicePrecision{T}((num, denom))

A number with twice the precision of `T`, e.g., quad-precision if `T =
Float64`. `hi` represents the high bits (most significant bits) and
`lo` the low bits (least significant bits). Rational values
`num//denom` can be approximated conveniently using the syntax
`TwicePrecision{T}((num, denom))`.

When used with `T<:Union{Float16,Float32,Float64}` to construct an "exact"
`StepRangeLen`, `ref` should be the range element with smallest
magnitude and `offset` set to the corresponding index.  For
efficiency, multiplication of `step` by the index is not performed at
twice precision: `step.hi` should have enough trailing zeros in its
`bits` representation that `(0:len-1)*step.hi` is exact (has no
roundoff error).  If `step` has an exact rational representation
`num//denom`, then you can construct `step` using

    step = TwicePrecision{T}((num, denom), nb)

where `nb` is the number of trailing zero bits of `step.hi`.  For
ranges, you can set `nb = ceil(Int, log2(len-1))`.
"""
struct TwicePrecision{T}
    hi::T    # most significant bits
    lo::T    # least significant bits
end

TwicePrecision{T}(x::T) where {T} = TwicePrecision{T}(x, zero(T))

function TwicePrecision{T}(x) where {T}
    xT = convert(T, x)
    Δx = x - xT
    TwicePrecision{T}(xT, T(Δx))
end

TwicePrecision{T}(i::Integer) where {T<:AbstractFloat} =
    TwicePrecision{T}(canonicalize2(splitprec(T, i)...)...)

TwicePrecision(x) = TwicePrecision{typeof(x)}(x)

# Numerator/Denominator constructors
function TwicePrecision{T}(nd::Tuple{Integer,Integer}) where {T<:Union{Float16,Float32}}
    n, d = nd
    TwicePrecision{T}(n/d)
end

function TwicePrecision{T}(nd::Tuple{Any,Any}) where {T}
    n, d = nd
    TwicePrecision{T}(n) / d
end

function TwicePrecision{T}(nd::Tuple{I,I}, nb::Integer) where {T,I}
    twiceprecision(TwicePrecision{T}(nd), nb)
end

# Truncating constructors. Useful for generating values that can be
# exactly multiplied by small integers.
function twiceprecision(val::T, nb::Integer) where {T<:IEEEFloat}
    hi = truncbits(val, nb)
    TwicePrecision{T}(hi, val - hi)
end

function twiceprecision(val::TwicePrecision{T}, nb::Integer) where {T<:IEEEFloat}
    hi = truncbits(val.hi, nb)
    TwicePrecision{T}(hi, (val.hi - hi) + val.lo)
end

nbitslen(r::StepRangeLen) = nbitslen(eltype(r), length(r), r.offset)
nbitslen(::Type{T}, len, offset) where {T<:IEEEFloat} =
    min(cld(precision(T), 2), nbitslen(len, offset))
# The +1 here is for safety, because the precision of the significand
# is 1 bit higher than the number that are explicitly stored.
nbitslen(len, offset) = len < 2 ? 0 : ceil(Int, log2(max(offset-1, len-offset))) + 1

eltype(::Type{TwicePrecision{T}}) where {T} = T

promote_rule(::Type{TwicePrecision{R}}, ::Type{TwicePrecision{S}}) where {R,S} =
    TwicePrecision{promote_type(R,S)}
promote_rule(::Type{TwicePrecision{R}}, ::Type{S}) where {R,S<:Number} =
    TwicePrecision{promote_type(R,S)}

(::Type{T})(x::TwicePrecision) where {T<:Number} = T(x.hi + x.lo)::T
TwicePrecision{T}(x::Number) where {T} = TwicePrecision{T}(T(x), zero(T))

convert(::Type{TwicePrecision{T}}, x::TwicePrecision{T}) where {T} = x
convert(::Type{TwicePrecision{T}}, x::TwicePrecision) where {T} =
    TwicePrecision{T}(convert(T, x.hi), convert(T, x.lo))

convert(::Type{T}, x::TwicePrecision) where {T<:Number} = T(x)
convert(::Type{TwicePrecision{T}}, x::Number) where {T} = TwicePrecision{T}(x)

float(x::TwicePrecision{<:AbstractFloat}) = x
float(x::TwicePrecision) = TwicePrecision(float(x.hi), float(x.lo))

big(x::TwicePrecision) = big(x.hi) + big(x.lo)

-(x::TwicePrecision) = TwicePrecision(-x.hi, -x.lo)

function zero(::Type{TwicePrecision{T}}) where {T}
    z = zero(T)
    TwicePrecision{T}(z, z)
end

# Arithmetic

function +(x::TwicePrecision, y::Number)
    s_hi, s_lo = add12(x.hi, y)
    TwicePrecision(canonicalize2(s_hi, s_lo+x.lo)...)
end
+(x::Number, y::TwicePrecision) = y+x

function +(x::TwicePrecision{T}, y::TwicePrecision{T}) where T
    r = x.hi + y.hi
    s = abs(x.hi) > abs(y.hi) ? (((x.hi - r) + y.hi) + y.lo) + x.lo : (((y.hi - r) + x.hi) + x.lo) + y.lo
    TwicePrecision(canonicalize2(r, s)...)
end
+(x::TwicePrecision, y::TwicePrecision) = +(promote(x, y)...)

-(x::TwicePrecision, y::TwicePrecision) = x + (-y)
-(x::TwicePrecision, y::Number) = x + (-y)
-(x::Number, y::TwicePrecision) = x + (-y)

function *(x::TwicePrecision, v::Number)
    v == 0 && return TwicePrecision(x.hi*v, x.lo*v)
    x * TwicePrecision(oftype(x.hi*v, v))
end
function *(x::TwicePrecision{<:IEEEFloat}, v::Integer)
    v == 0 && return TwicePrecision(x.hi*v, x.lo*v)
    nb = ceil(Int, log2(abs(v)))
    u = truncbits(x.hi, nb)
    TwicePrecision(canonicalize2(u*v, ((x.hi-u) + x.lo)*v)...)
end
*(v::Number, x::TwicePrecision) = x*v

function *(x::TwicePrecision{T}, y::TwicePrecision{T}) where {T}
    zh, zl = mul12(x.hi, y.hi)
    ret = TwicePrecision{T}(canonicalize2(zh, (x.hi * y.lo + x.lo * y.hi) + zl)...)
    ifelse(iszero(zh) | !isfinite(zh), TwicePrecision{T}(zh, zh), ret)
end
*(x::TwicePrecision, y::TwicePrecision) = *(promote(x, y)...)

function /(x::TwicePrecision, v::Number)
    x / TwicePrecision(oftype(x.hi/v, v))
end

function /(x::TwicePrecision, y::TwicePrecision)
    hi = x.hi / y.hi
    uh, ul = mul12(hi, y.hi)
    lo = ((((x.hi - uh) - ul) + x.lo) - hi*y.lo)/y.hi
    ret = TwicePrecision(canonicalize2(hi, lo)...)
    ifelse(iszero(hi) | !isfinite(hi), TwicePrecision(hi, hi), ret)
end

## StepRangeLen

# Use TwicePrecision only for Float64; use Float64 for T<:Union{Float16,Float32}
# See also _linspace1
# Ratio-of-integers constructors
function steprangelen_hp(::Type{Float64}, ref::Tuple{Integer,Integer},
                         step::Tuple{Integer,Integer}, nb::Integer,
                         len::Integer, offset::Integer)
    StepRangeLen(TwicePrecision{Float64}(ref),
                 TwicePrecision{Float64}(step, nb), Int(len), offset)
end

function steprangelen_hp(::Type{T}, ref::Tuple{Integer,Integer},
                         step::Tuple{Integer,Integer}, nb::Integer,
                         len::Integer, offset::Integer) where {T<:IEEEFloat}
    StepRangeLen{T}(ref[1]/ref[2], step[1]/step[2], Int(len), offset)
end

# AbstractFloat constructors (can supply a single number or a 2-tuple
const F_or_FF = Union{AbstractFloat, Tuple{AbstractFloat,AbstractFloat}}
asF64(x::AbstractFloat) = Float64(x)
asF64(x::Tuple{AbstractFloat,AbstractFloat}) = Float64(x[1]) + Float64(x[2])

function steprangelen_hp(::Type{Float64}, ref::F_or_FF,
                         step::F_or_FF, nb::Integer,
                         len::Integer, offset::Integer)
    StepRangeLen(TwicePrecision{Float64}(ref...),
                 twiceprecision(TwicePrecision{Float64}(step...), nb), Int(len), offset)
end

function steprangelen_hp(::Type{T}, ref::F_or_FF,
                         step::F_or_FF, nb::Integer,
                         len::Integer, offset::Integer) where {T<:IEEEFloat}
    StepRangeLen{T}(asF64(ref),
                    asF64(step), Int(len), offset)
end

StepRangeLen(ref::TwicePrecision{T}, step::TwicePrecision{T},
             len::Integer, offset::Integer=1) where {T} =
    StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}(ref, step, len, offset)

# Construct range for rational start=start_n/den, step=step_n/den
function floatrange(::Type{T}, start_n::Integer, step_n::Integer, len::Integer, den::Integer) where T
    if len < 2 || step_n == 0
        return steprangelen_hp(T, (start_n, den), (step_n, den), 0, Int(len), 1)
    end
    # index of smallest-magnitude value
    imin = clamp(round(Int, -start_n/step_n+1), 1, Int(len))
    # Compute smallest-magnitude element to 2x precision
    ref_n = start_n+(imin-1)*step_n  # this shouldn't overflow, so don't check
    nb = nbitslen(T, len, imin)
    steprangelen_hp(T, (ref_n, den), (step_n, den), nb, Int(len), imin)
end

function floatrange(a::AbstractFloat, st::AbstractFloat, len::Real, divisor::AbstractFloat)
    T = promote_type(typeof(a), typeof(st), typeof(divisor))
    m = maxintfloat(T, Int)
    if abs(a) <= m && abs(st) <= m && abs(divisor) <= m
        ia, ist, idivisor = round(Int, a), round(Int, st), round(Int, divisor)
        if ia == a && ist == st && idivisor == divisor
            # We can return the high-precision range
            return floatrange(T, ia, ist, Int(len), idivisor)
        end
    end
    # Fallback (misses the opportunity to set offset different from 1,
    # but otherwise this is still high-precision)
    steprangelen_hp(T, (a,divisor), (st,divisor), nbitslen(T, len, 1), Int(len), 1)
end

function (:)(start::T, step::T, stop::T) where T<:Union{Float16,Float32,Float64}
    step == 0 && throw(ArgumentError("range step cannot be zero"))
    # see if the inputs have exact rational approximations (and if so,
    # perform all computations in terms of the rationals)
    step_n, step_d = rat(step)
    if step_d != 0 && T(step_n/step_d) == step
        start_n, start_d = rat(start)
        stop_n, stop_d = rat(stop)
        if start_d != 0 && stop_d != 0 &&
                T(start_n/start_d) == start && T(stop_n/stop_d) == stop
            den = lcm_unchecked(start_d, step_d) # use same denominator for start and step
            m = maxintfloat(T, Int)
            if den != 0 && abs(start*den) <= m && abs(step*den) <= m &&  # will round succeed?
                    rem(den, start_d) == 0 && rem(den, step_d) == 0      # check lcm overflow
                start_n = round(Int, start*den)
                step_n = round(Int, step*den)
                len = max(0, div(den*stop_n - stop_d*start_n + step_n*stop_d, step_n*stop_d))
                # Integer ops could overflow, so check that this makes sense
                if isbetween(start, start + (len-1)*step, stop + step/2) &&
                        !isbetween(start, start + len*step, stop)
                    # Return a 2x precision range
                    return floatrange(T, start_n, step_n, len, den)
                end
            end
        end
    end
    # Fallback, taking start and step literally
    lf = (stop-start)/step
    if lf < 0
        len = 0
    elseif lf == 0
        len = 1
    else
        len = round(Int, lf) + 1
        stop′ = start + (len-1)*step
        # if we've overshot the end, subtract one:
        len -= (start < stop < stop′) + (start > stop > stop′)
    end
    steprangelen_hp(T, start, step, 0, len, 1)
end

step(r::StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}) where {T<:AbstractFloat} = T(r.step)
step(r::StepRangeLen{T,TwicePrecision{T},TwicePrecision{T}}) where {T} = T(r.step)

function range_start_step_length(a::T, st::T, len::Integer) where T<:Union{Float16,Float32,Float64}
    start_n, start_d = rat(a)
    step_n, step_d = rat(st)
    if start_d != 0 && step_d != 0 &&
            T(start_n/start_d) == a && T(step_n/step_d) == st
        den = lcm_unchecked(start_d, step_d)
        m = maxintfloat(T, Int)
        if abs(den*a) <= m && abs(den*st) <= m &&
                rem(den, start_d) == 0 && rem(den, step_d) == 0
            start_n = round(Int, den*a)
            step_n = round(Int, den*st)
            return floatrange(T, start_n, step_n, len, den)
        end
    end
    steprangelen_hp(T, a, st, 0, len, 1)
end

# This assumes that r.step has already been split so that (0:len-1)*r.step.hi is exact
function unsafe_getindex(r::StepRangeLen{T,<:TwicePrecision,<:TwicePrecision}, i::Integer) where T
    # Very similar to _getindex_hiprec, but optimized to avoid a 2nd call to add12
    @_inline_meta
    u = i - r.offset
    shift_hi, shift_lo = u*r.step.hi, u*r.step.lo
    x_hi, x_lo = add12(r.ref.hi, shift_hi)
    T(x_hi + (x_lo + (shift_lo + r.ref.lo)))
end

function _getindex_hiprec(r::StepRangeLen{<:Any,<:TwicePrecision,<:TwicePrecision}, i::Integer)
    u = i - r.offset
    shift_hi, shift_lo = u*r.step.hi, u*r.step.lo
    x_hi, x_lo = add12(r.ref.hi, shift_hi)
    x_hi, x_lo = add12(x_hi, x_lo + (shift_lo + r.ref.lo))
    TwicePrecision(x_hi, x_lo)
end

function getindex(r::StepRangeLen{T,<:TwicePrecision,<:TwicePrecision}, s::OrdinalRange{<:Integer}) where T
    @boundscheck checkbounds(r, s)
    soffset = 1 + round(Int, (r.offset - first(s))/step(s))
    soffset = clamp(soffset, 1, length(s))
    ioffset = first(s) + (soffset-1)*step(s)
    if step(s) == 1 || length(s) < 2
        newstep = r.step
    else
        newstep = twiceprecision(r.step*step(s), nbitslen(T, length(s), soffset))
    end
    if ioffset == r.offset
        StepRangeLen(r.ref, newstep, length(s), max(1,soffset))
    else
        StepRangeLen(r.ref + (ioffset-r.offset)*r.step, newstep, length(s), max(1,soffset))
    end
end

*(x::Real, r::StepRangeLen{<:Real,<:TwicePrecision}) =
    StepRangeLen(x*r.ref, twiceprecision(x*r.step, nbitslen(r)), length(r), r.offset)
*(r::StepRangeLen{<:Real,<:TwicePrecision}, x::Real) = x*r
/(r::StepRangeLen{<:Real,<:TwicePrecision}, x::Real) =
    StepRangeLen(r.ref/x, twiceprecision(r.step/x, nbitslen(r)), length(r), r.offset)

StepRangeLen{T,R,S}(r::StepRangeLen{T,R,S}) where {T<:AbstractFloat,R<:TwicePrecision,S<:TwicePrecision} = r

StepRangeLen{T,R,S}(r::StepRangeLen) where {T<:AbstractFloat,R<:TwicePrecision,S<:TwicePrecision} =
    _convertSRL(StepRangeLen{T,R,S}, r)

StepRangeLen{Float64}(r::StepRangeLen) =
    _convertSRL(StepRangeLen{Float64,TwicePrecision{Float64},TwicePrecision{Float64}}, r)
StepRangeLen{T}(r::StepRangeLen) where {T<:IEEEFloat} =
    _convertSRL(StepRangeLen{T,Float64,Float64}, r)

StepRangeLen{Float64}(r::AbstractRange) =
    _convertSRL(StepRangeLen{Float64,TwicePrecision{Float64},TwicePrecision{Float64}}, r)
StepRangeLen{T}(r::AbstractRange) where {T<:IEEEFloat} =
    _convertSRL(StepRangeLen{T,Float64,Float64}, r)

function _convertSRL(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen{<:Integer}) where {T,R,S}
    StepRangeLen{T,R,S}(R(r.ref), S(r.step), length(r), r.offset)
end

function _convertSRL(::Type{StepRangeLen{T,R,S}}, r::AbstractRange{<:Integer}) where {T,R,S}
    StepRangeLen{T,R,S}(R(first(r)), S(step(r)), length(r))
end

function _convertSRL(::Type{StepRangeLen{T,R,S}}, r::AbstractRange{U}) where {T,R,S,U}
    # if start and step have a rational approximation in the old type,
    # then we transfer that rational approximation to the new type
    f, s = first(r), step(r)
    start_n, start_d = rat(f)
    step_n, step_d = rat(s)
    if start_d != 0 && step_d != 0 &&
            U(start_n/start_d) == f && U(step_n/step_d) == s
        den = lcm_unchecked(start_d, step_d)
        m = maxintfloat(T, Int)
        if den != 0 && abs(f*den) <= m && abs(s*den) <= m &&
                rem(den, start_d) == 0 && rem(den, step_d) == 0
            start_n = round(Int, f*den)
            step_n = round(Int, s*den)
            return floatrange(T, start_n, step_n, length(r), den)
        end
    end
    __convertSRL(StepRangeLen{T,R,S}, r)
end

function __convertSRL(::Type{StepRangeLen{T,R,S}}, r::StepRangeLen{U}) where {T,R,S,U}
    StepRangeLen{T,R,S}(R(r.ref), S(r.step), length(r), r.offset)
end
function __convertSRL(::Type{StepRangeLen{T,R,S}}, r::AbstractRange{U}) where {T,R,S,U}
    StepRangeLen{T,R,S}(R(first(r)), S(step(r)), length(r))
end

function sum(r::StepRangeLen)
    l = length(r)
    # Compute the contribution of step over all indices.
    # Indexes on opposite side of r.offset contribute with opposite sign,
    #    r.step * (sum(1:np) - sum(1:nn))
    np, nn = l - r.offset, r.offset - 1  # positive, negative
    # To prevent overflow in sum(1:n), multiply its factors by the step
    sp, sn = sumpair(np), sumpair(nn)
    W = widen(Int)
    Δn = W(sp[1]) * W(sp[2]) - W(sn[1]) * W(sn[2])
    s = r.step * Δn
    # Add in contributions of ref
    ref = r.ref * l
    convert(eltype(r), s + ref)
end
function sum(r::StepRangeLen{<:Any,<:TwicePrecision,<:TwicePrecision})
    l = length(r)
    # Compute the contribution of step over all indices.
    # Indexes on opposite side of r.offset contribute with opposite sign,
    #    r.step * (sum(1:np) - sum(1:nn))
    np, nn = l - r.offset, r.offset - 1  # positive, negative
    # To prevent overflow in sum(1:n), multiply its factors by the step
    sp, sn = sumpair(np), sumpair(nn)
    tp = _tp_prod(r.step, sp[1], sp[2])
    tn = _tp_prod(r.step, sn[1], sn[2])
    s_hi, s_lo = add12(tp.hi, -tn.hi)
    s_lo += tp.lo - tn.lo
    # Add in contributions of ref
    ref = r.ref * l
    sm_hi, sm_lo = add12(s_hi, ref.hi)
    add12(sm_hi, sm_lo + ref.lo)[1]
end

# sum(1:n) as a product of two integers
sumpair(n::Integer) = iseven(n) ? (n+1, n>>1) : (n, (n+1)>>1)

function +(r1::StepRangeLen{T,R}, r2::StepRangeLen{T,R}) where T where R<:TwicePrecision
    len = length(r1)
    (len == length(r2) ||
        throw(DimensionMismatch("argument dimensions must match")))
    if r1.offset == r2.offset
        imid = r1.offset
        ref = r1.ref + r2.ref
    else
        imid = round(Int, (r1.offset+r2.offset)/2)
        ref1mid = _getindex_hiprec(r1, imid)
        ref2mid = _getindex_hiprec(r2, imid)
        ref = ref1mid + ref2mid
    end
    step = twiceprecision(r1.step + r2.step, nbitslen(T, len, imid))
    StepRangeLen{T,typeof(ref),typeof(step)}(ref, step, len, imid)
end

## LinRange

# For Float16, Float32, and Float64, this returns a StepRangeLen
function range_start_stop_length(start::T, stop::T, len::Integer) where {T<:IEEEFloat}
    len < 2 && return _linspace1(T, start, stop, len)
    if start == stop
        return steprangelen_hp(T, start, zero(T), 0, len, 1)
    end
    # Attempt to find exact rational approximations
    start_n, start_d = rat(start)
    stop_n, stop_d = rat(stop)
    if start_d != 0 && stop_d != 0
        den = lcm_unchecked(start_d, stop_d)
        m = maxintfloat(T, Int)
        if den != 0 && abs(den*start) <= m && abs(den*stop) <= m
            start_n = round(Int, den*start)
            stop_n = round(Int, den*stop)
            if T(start_n/den) == start && T(stop_n/den) == stop
                return _linspace(T, start_n, stop_n, len, den)
            end
        end
    end
    _linspace(start, stop, len)
end

function _linspace(start::T, stop::T, len::Integer) where {T<:IEEEFloat}
    (isfinite(start) && isfinite(stop)) || throw(ArgumentError("start and stop must be finite, got $start and $stop"))
    # Find the index that returns the smallest-magnitude element
    Δ, Δfac = stop-start, 1
    if !isfinite(Δ)   # handle overflow for large endpoints
        Δ, Δfac = stop/len - start/len, Int(len)
    end
    tmin = -(start/Δ)/Δfac            # t such that (1-t)*start + t*stop == 0
    imin = round(Int, tmin*(len-1)+1) # index approximately corresponding to t
    if 1 < imin < len
        # The smallest-magnitude element is in the interior
        t = (imin-1)/(len-1)
        ref = T((1-t)*start + t*stop)
        step = imin-1 < len-imin ? (ref-start)/(imin-1) : (stop-ref)/(len-imin)
    elseif imin <= 1
        imin = 1
        ref = start
        step = (Δ/(len-1))*Δfac
    else
        imin = Int(len)
        ref = stop
        step = (Δ/(len-1))*Δfac
    end
    if len == 2 && !isfinite(step)
        # For very large endpoints where step overflows, exploit the
        # split-representation to handle the overflow
        return steprangelen_hp(T, start, (-start, stop), 0, 2, 1)
    end
    # 2x calculations to get high precision endpoint matching while also
    # preventing overflow in ref_hi+(i-offset)*step_hi
    m, k = prevfloat(floatmax(T)), max(imin-1, len-imin)
    step_hi_pre = clamp(step, max(-(m+ref)/k, (-m+ref)/k), min((m-ref)/k, (m+ref)/k))
    nb = nbitslen(T, len, imin)
    step_hi = truncbits(step_hi_pre, nb)
    x1_hi, x1_lo = add12((1-imin)*step_hi, ref)
    x2_hi, x2_lo = add12((len-imin)*step_hi, ref)
    a, b = (start - x1_hi) - x1_lo, (stop - x2_hi) - x2_lo
    step_lo = (b - a)/(len - 1)
    ref_lo = a - (1 - imin)*step_lo
    steprangelen_hp(T, (ref, ref_lo), (step_hi, step_lo), 0, Int(len), imin)
end

# range for rational numbers, start = start_n/den, stop = stop_n/den
# Note this returns a StepRangeLen
_linspace(::Type{T}, start::Integer, stop::Integer, len::Integer) where {T<:IEEEFloat} = _linspace(T, start, stop, len, 1)
function _linspace(::Type{T}, start_n::Integer, stop_n::Integer, len::Integer, den::Integer) where T<:IEEEFloat
    len < 2 && return _linspace1(T, start_n/den, stop_n/den, len)
    start_n == stop_n && return steprangelen_hp(T, (start_n, den), (zero(start_n), den), 0, len, 1)
    tmin = -start_n/(Float64(stop_n) - Float64(start_n))
    imin = round(Int, tmin*(len-1)+1)
    imin = clamp(imin, 1, Int(len))
    ref_num = Int128(len-imin) * start_n + Int128(imin-1) * stop_n
    ref_denom = Int128(len-1) * den
    ref = (ref_num, ref_denom)
    step_full = (Int128(stop_n) - Int128(start_n), ref_denom)
    steprangelen_hp(T, ref, step_full,  nbitslen(T, len, imin), Int(len), imin)
end

# For len < 2
function _linspace1(::Type{T}, start, stop, len::Integer) where T<:IEEEFloat
    len >= 0 || throw(ArgumentError("range($start, stop=$stop, length=$len): negative length"))
    if len <= 1
        len == 1 && (start == stop || throw(ArgumentError("range($start, stop=$stop, length=$len): endpoints differ")))
        # Ensure that first(r)==start and last(r)==stop even for len==0
        # The output type must be consistent with steprangelen_hp
        if T<:Union{Float32,Float16}
            return StepRangeLen{T}(Float64(start), Float64(start) - Float64(stop), len, 1)
        else
            return StepRangeLen(TwicePrecision(start, zero(T)), TwicePrecision(start, -stop), len, 1)
        end
    end
    throw(ArgumentError("should only be called for len < 2, got $len"))
end

### Numeric utilities

# Approximate x with a rational representation. Guaranteed to return,
# but not guaranteed to return a precise answer.
# https://en.wikipedia.org/wiki/Continued_fraction#Best_rational_approximations
function rat(x)
    y = x
    a = d = 1
    b = c = 0
    m = maxintfloat(narrow(typeof(x)), Int)
    while abs(y) <= m
        f = trunc(Int,y)
        y -= f
        a, c = f*a + c, a
        b, d = f*b + d, b
        max(abs(a), abs(b)) <= convert(Int,m) || return c, d
        oftype(x,a)/oftype(x,b) == x && break
        y = inv(y)
    end
    return a, b
end

# This version of lcm does not check for overflows
lcm_unchecked(a::T, b::T) where T<:Integer = a * div(b, gcd(a, b))

narrow(::Type{T}) where {T<:AbstractFloat} = Float64
narrow(::Type{Float64}) = Float32
narrow(::Type{Float32}) = Float16
narrow(::Type{Float16}) = Float16

function _tp_prod(t::TwicePrecision, x, y...)
    @_inline_meta
    _tp_prod(t * x, y...)
end
_tp_prod(t::TwicePrecision) = t
<(x::TwicePrecision{T}, y::TwicePrecision{T}) where {T} =
    x.hi < y.hi || ((x.hi == y.hi) & (x.lo < y.lo))

isbetween(a, x, b) = a <= x <= b || b <= x <= a

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
module Quadrature
using RandomNumbers

# Functions we will integrate
f1(x) = (1.0 / sqrt(2 * pi)) * exp((-x^2.0) / 2.0)
f1lim = (0.0, 1.0)

f2(x) = ((sin(2.0 * x))^2.0)/(x^2.0)
f2lim = (pi, 2 * pi)

f3(x) = x^(-.5)
f3lim = (0.0, 2.0)

# Change of variable for f1 -> Inf
f4(x) = f1(-log(x)) / x
f4lim = (0.0, e^-f1lim[1])

# Infinite limits of integration for printing
f1inflim = (0.0, Inf)

# PDF of the normal distribution
function normpdf(x, mu=0, sigma=1)
	return (sigma * sqrt(2 * pi))^-1.0 * e^(-(x - mu)^2 / (2 * sigma^2))
end

# Limits to integrate normpdf over
normlim = (-50.0, 50.0)

# PDF of the Cauchy distribution
cauchypdf(x) = (pi * (1.0 + x^2.0))^-1.0

# PDF of the t-distribution
function tpdf(x, dof)
	term1 = gamma(.5 * (dof + 1.0)) / (sqrt(dof * pi) * gamma(dof / 2.0))
	term2 = (1.0 + (x^2.0)/dof)^(-.5*(dof+1.0))
	return term1 * term2
end

# Generate a single random t variable
function rn_t(dof)
	z = rn_standard_normal()
	v = rn_chi_sq(dof)
	return z / sqrt(v / dof)
end

# Utility function to print each entry in an array on its own line
mapln(x) = map(println, x)

# Utility function to truncate
function truncate(x::Float64, minvalue::Float64, maxvalue::Float64)
	if x > maxvalue
		return maxvalue
	elseif x < minvalue
		return minvalue
	else
		return x
	end
end

# Truncated Cauchy
rn_trunc_cauchy(a, b) = truncate(rn_standard_cauchy(), a, b)
rn_trunc_cauchy(a, b, n) = [rn_trunc_cauchy(a, b) for x in 1:n]

# Truncated t
rn_trunc_t(dof, a, b) = truncate(rn_t(dof), a, b)
rn_trunc_t(dof, a, b, n) = [rn_trunc_t(dof, a, b) for x in 1:n]

# Extended trapezoidal rule
function exttrap(f, a, b; nmax=20, epsilon=1e-6)
	# Initialize I
	n = 2
	I = [(b - a) * .5 * ( f(a) + f(b))]

	while n < nmax
		# Determine new points to evaluate at
		irange = 1:2^(n-2)
		evalpoints = map(x -> (a + (2*x - 1) * (b - a) / (2^(n -1))), irange)

		# Sum function evals at new points
		sumterm = sum(map(f, evalpoints))

		# Add new I onto results vector
		push!(I, .5 * I[n-1] + (sumterm * (b - a) / (2^(n-1))))

		# Check for convergence
		if abs(I[n] - I[n-1]) <= epsilon * abs(I[n-1])
			return I
		end

		n += 1
	end

	return I
end

# Romberg integration
function romberg(f, a, b; nmax=20, epsilon=1e-6)
	# Get the results of the trapezoidal rule
	traps = exttrap(f, a, b, nmax=nmax, epsilon=epsilon)

	# Find a size for our results
	resultsize = min(nmax, length(traps))

	# Create a zero matrix and assign traps to the first column
	r = zeros(Float64, (resultsize, resultsize))
	r[:, 1] = traps[1:resultsize]

	# For each row past the first
	for i in 2:resultsize
		# For each column past the first
		for j in 2:i
			r[i,j] = (1 / (4^(j - 1) - 1)) * (4^(j-1) * r[i,j-1] - r[i-1,j-1])
		end

		# Check stopping criteria
		if abs(r[i,i] - r[i,i-1]) <= epsilon * abs(r[i,i-1])
			return r[1:i,1:i]
		end
	end
end

# Extended midpoint rule
function midpoint(f, a, b; nmax=20, epsilon=1e-6)
	# Initialize I
	n = 2
	I = [f(.5 * (a + b))]

	while n < nmax
		# Determine new points to evaluate at:
		# Keep odd numbers from 1:2*3^(n-1) that aren't divisible by 3
		newodds = filter(x -> x % 3 != 0, [1:2:2*3^(n-1)])

		# Divide by (2*3^(n-1)) and scale to (a,b)
		newpoints = a + (b - a) * newodds / (2*3^(n - 1))

		# Sum function evals at this point
		sumterm = sum(map(f, newpoints))

		# Add new I onto results vector
		push!(I, I[n-1] / 3.0 + sumterm * (b - a) / 3.0^(n-1))

		# Check for convergence
		if abs(I[n] - I[n-1]) <= epsilon * abs(I[n-1])
			return I
		end

		n += 1
	end

	return I
end

# Monte Carlo Integration with antithetic variables
function mci(f, a, b; nmax=1e4, epsilon=1e-6)

	# Store our random variables in a vector
	x = [rn_unif(a, b)]

	# Calculate our initial result
	I = [(b - a) * .5 * (f(x[1]) + f(a + b - x[1]))]

	# Initialize n
	n = 2
	while n < nmax
		# Generate a new random variable and add it to the vector
		push!(x, rn_unif(a, b))

		# Calculate new result and push it onto the results vector
		push!(I, (b - a) * .5 * mean(map(f, x) + map(f, a + b - x)))

		# Check for convergence
		if abs(I[n] - I[n-1]) <= epsilon * abs(I[n-1])
			return I
		end

		n += 1
	end

	return I
end

# Monte Carlo Integration with Importance Sampling (Uniform)
function mciunif(f, a, b; n=1000)
	# Store our random variables in a vector
	x = rn_unif(a, b, n)

	# Calculate our estimate
	I = (b - a) * mean(map(f, x))

	return I
end

# Monte Carlo Integration with Importance Sampling (Truncated Standard Cauchy)
function mcicauchy(f, a, b; n=1000)
	# Store our random variables in a vector
	x = rn_trunc_cauchy(a, b, n)

	# Calculate our estimate
	I = mean(map(f, x) ./ map(cauchypdf, x))

	return I
end

# Monte Carlo Integration with Importance Sampling (Truncated t-distribution)
function mcit(f, a, b; n=1000)
	# Use 30 degrees of freedom - in the future this might be worth parameterizing
	dof = 30

	# Store our random variables in a vector
	x = rn_trunc_t(dof, a, b, n)

	# Calculate our estimate
	I = mean(map(f, x) ./ map(x -> tpdf(x, dof), x))

	return I
end

# Run integration function f n times and write final results to a csv file
function mcincsv(f, n, filename)
	println("Writing $n runs to $filename")
	results = [f(x) for x in 1:n]
	writecsv(filename, results)
end

end

max_stars_repo_name: andrewjlm/StatComputing

max_issues_repo_name: andrewjlm/StatComputing

max_forks_repo_name: andrewjlm/StatComputing


text:
module TestModel

using Mimi

include("region_parameters.jl")
include("gross_economy.jl")
include("emissions.jl")

export create_model

function create_model()

	m = Model()

	set_dimension!(m, :time, collect(2015:5:2110))
	set_dimension!(m, :regions, [:Region1, :Region2, :Region3])

	add_comp!(m, grosseconomy)
	add_comp!(m, emissions)

	# set parameters for grosseconomy component
	update_param!(m, :grosseconomy, :l, l)
	update_param!(m, :grosseconomy, :tfp, tfp)
	update_param!(m, :grosseconomy, :s, s)
	update_param!(m, :grosseconomy, :depk,depk)
	update_param!(m, :grosseconomy, :k0, k0)
	update_param!(m, :grosseconomy, :share, 0.3)
	update_param!(m, :grosseconomy, :tester, zeros(Mimi.dim_count(m.md, :time),
												Mimi.dim_count(m.md, :regions)))

	# set and connect parameters for emissions component
	update_param!(m, :emissions, :sigma, sigma)
	connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    return m
end

end # module


text:
function zad36(A)
    a,b = size(A)
    s = 1
    if( a!=b)
        return NaN
    end
    for i in 1:a
        for j in 1:a
            if i == j
                s *= A[i,j]
            end
        end
    end
    return s
end


text:
include("runtests-prerequisite.jl")

mutable struct MyMovie <: PostgresORM.IEntity

  title::Union{Missing,String}
  year::Union{Missing,Int32}

  MyMovie(args::NamedTuple) = MyMovie(;args...)
  MyMovie(;title = missing,
           year = missing,) = (
            x = new(missing, missing);
            x.title = title;
            x.year = year;
            return x
            )
end

mutable struct MyGreatStruct <: PostgresORM.IEntity

  id::Union{Missing,Int32}
  first_name::Union{Missing,String}
  favorite_movie::Union{Missing, MyMovie}

  MyGreatStruct(args::NamedTuple) = MyGreatStruct(;args...)
  MyGreatStruct(;
        id = missing,
        first_name = missing,
        favorite_movie = missing,
 ) = (
            x = new(missing, missing, missing);
            x.id = id;
            x.first_name = first_name;
            x.favorite_movie = favorite_movie;
            return x
            )

end

module MyGreatStructORM
    data_type = Main.MyGreatStruct
    Main.PostgresORM.get_orm(x::Main.MyGreatStruct) = return(MyGreatStructORM)
    get_table_name() = "public.my_new_struct"
    const columns_selection_and_mapping =
      Dict(
         :id => "id",
         :first_name => "first_name",
         :favorite_movie => ["favorite_movie_title","favorite_movie_year"]
      )
    get_id_props() = return [:id]

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict()

    const track_changes = false
end

module MyMovieORM
    data_type = Main.MyMovie
    Main.PostgresORM.get_orm(x::Main.MyMovie) = return(MyMovieORM)
    get_table_name() = "public.my_movie"
    const columns_selection_and_mapping =
      Dict(
         :title => "title",
         :year => "year"
      )
    get_id_props() = return [:title,:year]

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict()

    const track_changes = false
end

@testset "Test `util_convert_flatdictfromdb_to_structuredrenameddict`" begin

    flatdict = Dict(:id => 1,
                    :first_name=> "Vincent",
                    :favorite_movie_title => "The Godfather",
                    :favorite_movie_year => 1972)

    PostgresORM.Controller.util_convert_flatdictfromdb_to_structuredrenameddict(flatdict,
                                                         MyGreatStruct)
end


text:
using Test, IonSim, IonSim.PhysicalConstants, Suppressor

@testset "constants -- PhysicalConstant" begin
    # test that algebraic operations on physical constants return a <:Number
    @test typeof(c - c) <: Number
    @test typeof(c - 1) <: Number
    @test typeof(1 - c) <: Number
    @test typeof(c + c) <: Number
    @test typeof(1 + c) <: Number
    @test typeof(c + 1) <: Number
    @test typeof(c * c) <: Number
    @test typeof(1 * c) <: Number
    @test typeof(c * 1) <: Number
    @test typeof(c / c) <: Number
    @test typeof(1 / c) <: Number
    @test typeof(c / 1) <: Number
    @test typeof(c^3) <: Number
    @test typeof(c^α) <: Number
    @test typeof(2^α) <: Number
    @test typeof(sqrt(c)) <: Number
    @test typeof(α^α) <: Number

    # Test comparisons for physical constants
    @test c.x == c
    @test c == c.x
    @test c == c
    @test c.x > c-1
    @test c > c.x-1
    @test c > c-1
    @test c.x < c+1
    @test c < c.x+1
    @test c < c+1
    @test c.x ≥ c
    @test c ≥ c.x
    @test c ≥ c
    @test c.x ≤ c
    @test c ≤ c.x
    @test c ≤ c
    @test α < c
    @test α > ħ

    # test print/show for physicalconstants
    out = @capture_out print(c)
    @test out == "$(c.x) [$(c.units)]"

    @test IonSim._print_axis(x̂) == "x̂"
    @test IonSim._print_axis(ŷ) == "ŷ"
    @test IonSim._print_axis(ẑ) == "ẑ"
    @test IonSim._print_axis((x=1,y=1,z=1)) == string((x=1,y=1,z=1))
end

@testset "constants -- other constants" begin
    @test x̂/2 == (x=1/2, y=0, z=0)
    @test 2x̂ == (x=2, y=0, z=0)
    @test x̂*2 == (x=2, y=0, z=0)
    @test x̂ + x̂ == (x=2, y=0, z=0)
    @test x̂ - x̂ == (x=0, y=0, z=0)
    @test ndot(x̂, x̂) == 1
    @test ndot(x̂, ŷ) == 0
end


text:
module DataDeps

import JSON
include("./segregate.jl")
export add, list, traindata, testdata, setseries

function add(url::AbstractString, name::AbstractString, dir::AbstractString="")
    if (dir=="")
        dir=joinpath(Pkg.dir("DataDeps"), "datasets/$(name)")
    end
    getdata(dir, name, url)
    return 0
end

function list(dir::AbstractString="")
    dirhere=joinpath(Pkg.dir("DataDeps"), "datasets")
    if (dir=="")
        dir=dirhere
    end
    dictall=Dict()
    try
        f=open(joinpath(dirhere, "downloads.json"), "r")
        dicttxt = readstring(f)
        close(f)
        dictall=JSON.parse(dicttxt)
    end
    textout=""
    for i in dictall
        if (dir==i[2][1:length(dir)])
            textout*="$(i[1])\n"
        end
    end
    println(textout)
    return textout
end

end #module


text:
const NATIVE_LABEL_ENCODING = LabelEnc.NativeLabels(["inlier", "outlier"])
const LABEL_ENCODING = LabelEnc.NativeLabels([:inlier,:outlier])
const LEARNING_LABEL_ENCODING = LabelEnc.NativeLabels([:Lin,:Lout])

function load_data(file_path; header=false, native_label_encoding=NATIVE_LABEL_ENCODING)
    raw_data, _ = header ? readdlm(file_path, ',', header=header) : (readdlm(file_path, ','), nothing)
    data = copy(transpose(float.(raw_data[:, 1:end-1])))
    labels = convert_labels_from_raw(raw_data[:,end])
    return data, labels
end

function convert_labels_from_raw(raw_labels, native_label_encoding=NATIVE_LABEL_ENCODING)
    return convertlabel(LABEL_ENCODING, raw_labels, native_label_encoding)
end

function convert_labels_to_learning(labels, label_encoding=LABEL_ENCODING)
    return convertlabel(LEARNING_LABEL_ENCODING, labels, label_encoding)
end

abstract type SplitStrategy end

struct FullSplitStrat <: SplitStrategy end
struct UnlabeledSplitStrat <: SplitStrategy end
struct UnlabeledAndLabeledInlierSplitStrat <: SplitStrategy end
struct LabeledSplitStrat <: SplitStrategy end
struct LabeledInlierSplitStrat <: SplitStrategy end
struct LabeledOutlierSplitStrat <: SplitStrategy end

struct DataSplits
    train::BitArray
    test::BitArray
    train_strat::SplitStrategy
    test_strat::SplitStrategy
    query_strat::SplitStrategy
    DataSplits(train::BitArray) = DataSplits(train, FullSplitStrat())
    DataSplits(train::BitArray, strat::SplitStrategy) = DataSplits(train, train, strat)
    DataSplits(train::BitArray, train_strat::SplitStrategy, test_strat::SplitStrategy) = new(train, train, train_strat, test_strat, train_strat)
    DataSplits(train::BitArray, train_strat::SplitStrategy, test_strat::SplitStrategy, query_strat::SplitStrategy) = new(train, train, train_strat, test_strat, query_strat)
    DataSplits(train::BitArray, test::BitArray, strat::SplitStrategy) = new(train, test, strat, strat, strat)
    DataSplits(train::BitArray, test::BitArray, train_strat::SplitStrategy, test_strat::SplitStrategy) = new(train, test, train_strat, test_strat, train_strat)
end

get_train(ds::DataSplits, data::Array{T, 2}, pools::Vector{Symbol}) where T <: Real = select_subset(ds.train_strat, ds.train, data, pools)
get_test(ds::DataSplits, data::Array{T, 2}, pools::Vector{Symbol}) where T <: Real = select_subset(ds.test_strat, ds.test, data, pools)
get_query(ds::DataSplits, data::Array{T, 2}, pools::Vector{Symbol}) where T <: Real = select_subset(ds.query_strat, ds.train, data, pools)

function select_subset(strat::SplitStrategy, init_mask::BitArray, data::Array{T, 2}, pools::Vector{Symbol}) where T <: Real
    mask = calc_mask(strat, init_mask, pools)
    return (data[:, mask], pools[mask], findall(mask))
end

calc_mask(strat::FullSplitStrat, init_mask::BitArray, pools::Vector{Symbol}) = init_mask
calc_mask(strat::UnlabeledSplitStrat, init_mask::BitArray, pools::Vector{Symbol}) = init_mask .& (pools .== :U)
calc_mask(strat::UnlabeledAndLabeledInlierSplitStrat, init_mask::BitArray, pools::Vector{Symbol}) = init_mask .& ((pools .== :U) .| (pools .== :Lin))
calc_mask(strat::LabeledSplitStrat, init_mask::BitArray, pools::Vector{Symbol}) = init_mask .& ((pools .== :Lin) .| (pools .== :Lout))
calc_mask(strat::LabeledInlierSplitStrat, init_mask::BitArray, pools::Vector{Symbol}) = init_mask .& (pools .== :Lin)
calc_mask(strat::LabeledOutlierSplitStrat, init_mask::BitArray, pools::Vector{Symbol}) = init_mask .& (pools .== :Lout)

function get_initial_pools(data, labels, data_splits, initial_pool_strategy; n=20, p=0.1, x=10)
    if initial_pool_strategy ∉ ["Pu", "Pp", "Pn", "Pa"]
        throw(ArgumentError("Unknown initial pools strategy '$(initial_pool_strategy)'."))
    end
    l = fill(:U, size(data, 2))
    if initial_pool_strategy == "Pa"
        n = min(size(data, 1) + x, size(data, 2))
        l = fill(:U, size(data, 2))
        label_candidates = findall(data_splits.train .& (labels .== :inlier))
        l[label_candidates[1:n]] .= convert_labels_to_learning(labels[label_candidates[1:n]])
    elseif initial_pool_strategy ∈ ["Pp", "Pn"]
        p = initial_pool_strategy == "Pp" ? p : min(1.0, n / sum(data_splits.train))
        label_candidates = findall(data_splits.train)
        (label_indices, _), _ = stratifiedobs((label_candidates, labels[data_splits.train]), p=p)
        l[label_indices] .= convert_labels_to_learning(labels[label_indices])
    end
    return l
end

function get_splits_and_init_pools(data, labels, split_strategy, initial_pool_strategy; n=20, p=0.1, x=10)
    train = trues(size(data, 2))
    if split_strategy == "Sf"
        data_splits = DataSplits(train, FullSplitStrat())
    elseif split_strategy == "Sh"
        (train_indices, train_labels), (test_indices, test_labels) = stratifiedobs((collect(1:length(labels)), labels), p=0.8)
        train[test_indices] .= false
        data_splits = DataSplits(train, .~train, FullSplitStrat())
    elseif split_strategy == "Si"
        data_splits = DataSplits(train, LabeledInlierSplitStrat(), FullSplitStrat(), FullSplitStrat())
    else
        throw(ArgumentError("Unknown split strategy '$(split_strategy)'."))
    end
    initial_pools = get_initial_pools(data, labels, data_splits, initial_pool_strategy, x=x)
    return data_splits, initial_pools
end


text:
# Tests for the code in the API.jl file.
# Author: Noah C. Benson <n@nben.net>

@testset "API" begin

    @testset "setindex" begin
        d1 = Dict()
        d2 = setindex(d1, 10, :a)
        @test d1 != d2
        @test d2[:a] == 10
        @test length(d1) == 0
        d1 = IdDict()
        d2 = setindex(d1, 10, :a)
        @test d1 != d2
        @test d2[:a] == 10
        @test length(d1) == 0
        u1 = [:x, :y, :z]
        u2 = setindex(u1, :a, 1)
        @test u1 != u2
        @test u2[1] == :a
        @test u1[1] == :x
        u1 = view(u1, 1:2)
        u2 = setindex(u1, :a, 1)
        @test u2[1] == :a
        @test u1[1] == :x
        u1 = reshape(u1, (1,2))
        u2 = setindex(u1, :a, 1, 1)
        @test u2[1,1] == :a
        @test u1[1,1] == :x
        t1 = (:x, :y, :z)
        t2 = setindex(t1, :a, 1)
        @test t2[1] == :a
        @test length(t2) == 3
        b1 = BitArray([1,1,1,1,1])
        b2 = setindex(b1, 0, 3)
        @test b1[3] == true
        @test b2[3] == false
    end

    @testset "PArray" begin
        a = collect(reshape(1:100, (10,10)))
        b = convert(Array{Float64,2}, collect(reshape(1:100, (10,10))))
        p = PArray(a)
        q = PArray(b)
        @test a == b
        @test p == q
        @test a == p
        @test a == q
        @test b == p
        @test b == q
        @test !(a === b)
        @test !(a === p)
        @test !(a === q)
        @test !(b === p)
        @test !(b === q)
        @test !(p === q)
        @test hash(a) == hash(b)
        @test hash(p) == hash(q)
        @test hash(a) == hash(p)
        @test hash(a) == hash(q)
        @test hash(b) == hash(p)
        @test hash(b) == hash(q)
    end

    @testset "PDict" begin
        intels = [gensym() => x for x in 1:2000]
        fltels = [k => Float64(v) for (k,v) in intels]
        for (DT,PDT) in ((Dict,      PDict),
                         (IdDict,    PIdDict))
            a = DT(intels...)
            b = DT(fltels...)
            p = PDT(a)
            q = PDT(b)
            @test a == b
            @test a == p
            @test a == q
            @test b == p
            @test b == q
            @test !(a === b)
            @test !(a === p)
            @test !(a === q)
            @test !(b === p)
            @test !(b === q)
            @test !(p === q)
            @test hash(a) == hash(b)
            @test hash(p) == hash(q)
            if DT === Dict
                @test hash(a) == hash(p)
                @test hash(a) == hash(q)
                @test hash(b) == hash(p)
                @test hash(b) == hash(q)
            end
        end
    end

    make_idset(u::Vector{T}) where {T} = begin
        s = IdSet{T}()
        for x in u
            push!(s, x)
        end
        return s
    end
    @testset "PSet" begin
        intels = collect(1:1000)
        fltels = [Float64(v) for v in intels]
        for (ST,PST) in ((Set, PSet), (IdSet, PIdSet))
            if ST === IdSet
                a = make_idset(intels)
                b = make_idset(fltels)
            else
                a = ST(intels)
                b = ST(fltels)
            end
            p = PST(a)
            q = PST(b)
            @test a == p
            @test b == q
            @test (a == b) == (p == q)
            @test !(a === b)
            @test !(a === p)
            @test !(a === q)
            @test !(b === p)
            @test !(b === q)
            @test !(p === q)
            @test hash(a) == hash(b)
            @test hash(p) == hash(q)
            if ST === Set
                @test hash(a) == hash(p)
                @test hash(a) == hash(q)
                @test hash(b) == hash(p)
                @test hash(b) == hash(q)
            end
        end
    end

    @testset "PWDict" begin
        intels = [gensym() => (x,rand(Float64)) for x in 1:2000]
        fltels = [k => (Float64(v),w) for (k,(v,w)) in intels]
        p = PWDict{Symbol,Int,Float64}(intels...)
        q = PWDict{Symbol,Float64,Float64}(fltels...)
        @test p == q
        @test !(p === q)
        @test hash(p) == hash(q)
    end

    @testset "PWSet" begin
        intels = [gensym() => rand(Float64) for x in 1:2000]
        p = PWSet{Symbol,Float64}(intels)
        q = PWSet{Symbol,Float64}(reverse(intels))
        @test p == q
        @test !(p === q)
        @test hash(p) == hash(q)
    end
end


text:
using AcuteML, Test
@testset "xmlutils" begin

    nsc = createnode(AbsEmpty, "some")

    @testset "Node" begin
        n = createnode(AbsNormal, "a")
        ################################################################
        addnode!(n, "normalString", "val", AbsNormal)
        @test "val" == findcontent(String, "normalString", n, AbsNormal)
        @test ["val"] == findcontent("normalString", n, AbsNormal)

        updatecontent!("val2", "normalString", n, AbsNormal)
        @test "val2" == findcontent(String, "normalString", n, AbsNormal)

        addnode!(n, "attString", "val", AbsAttribute)
        @test "val" == findcontent(String, "attString", n, AbsAttribute)
        updatecontent!("val2", "attString", n, AbsAttribute)
        @test "val2" == findcontent(String, "attString", n, AbsAttribute)
        @test "val2" == get(n, "attString", nothing)
        @test nothing == get(n, "attString1", nothing)

        using Dates
        addnode!(n, "Date", Date(2013,7,1), AbsNormal)
        @test Date(2013,7,1) == findcontent(Date, "Date", n, AbsNormal)
        updatecontent!(Date(2013,7,2), "Date", n, AbsNormal)
        @test Date(2013,7,2) == findcontent(Date, "Date", n, AbsNormal)

        addnode!(n, "Time", Time(12,53,40), AbsNormal)
        @test Time(12,53,40) == findcontent(Time, "Time", n, AbsNormal)
        updatecontent!(Time(12,53,41), "Time", n, AbsNormal)
        @test Time(12,53,41) == findcontent(Time, "Time", n, AbsNormal)

        addnode!(n, "DateTime", DateTime(2013,5,1,12,53,40), AbsNormal)
        @test DateTime(2013,5,1,12,53,40) == findcontent(DateTime, "DateTime", n, AbsNormal)
        updatecontent!( DateTime(2013,5,1,12,53,41), "DateTime", n, AbsNormal)
        @test DateTime(2013,5,1,12,53,41) == findcontent(DateTime, "DateTime", n, AbsNormal)

        using DataFrames
        addnode!(n, "DataFrame",  DataFrame(course = ["Artificial Intelligence", "Robotics"], professor = ["Prof. A", "Prof. B"] ), AbsNormal)
        @test_skip  DataFrame(course = ["Artificial Intelligence", "Robotics"], professor = ["Prof. A", "Prof. B"] ) ==  findcontent(DataFrame, "DataFrame", n, AbsNormal)
        # updatecontent!( DataFrame(course = ["Artificial Intelligence", "Robotics2"], professor = ["Prof. A", "Prof. B"] ),"DataFrame", n, AbsNormal)
        # @test_skip  DataFrame(course = ["Artificial Intelligence", "Robotics2"], professor = ["Prof. A", "Prof. B"] ) ==  findcontent(DataFrame, "DataFrame", n, AbsNormal)

        addnode!(n, "Nothing", nothing, AbsAttribute)
        @test nothing == findcontent(Nothing, "Nothing", n, AbsAttribute)
        @test nothing == findcontent(Union{Nothing,String}, "Nothing", n, AbsAttribute)

        # updatecontent!(nothing, "Nothing", n, AbsAttribute)
        # @test nothing == findcontent(Nothing, "Nothing", n, AbsAttribute)

        ################################################################
        addnode!(n, "stringVect", ["aa", "bb"], AbsNormal)
        @test ["aa", "bb"] == findcontent(Vector{String}, "stringVect", n, AbsNormal)
        updatecontent!(["aa2", "bb"], "stringVect", n, AbsNormal)
        @test ["aa2", "bb"] == findcontent(Vector{String}, "stringVect", n, AbsNormal)

        addnode!(n, "floatVect", [5.6, 7.8], AbsNormal)
        @test [5.6, 7.8] == findcontent(Vector{Float64}, "floatVect", n, AbsNormal)
        updatecontent!([5.6, 7.9], "floatVect", n, AbsNormal)
        @test [5.6, 7.9] == findcontent(Vector{Float64}, "floatVect", n, AbsNormal)

        addnode!(n, "intVect", [5, 6], AbsNormal)
        @test [5, 6] == findcontent(Vector{Int64}, "intVect", n, AbsNormal)
        updatecontent!([5, 7], "intVect", n, AbsNormal)
        @test [5, 7] == findcontent(Vector{Int64}, "intVect", n, AbsNormal)

        addnode!(n, "DateVect", [Date(2013,7,1), Date(2014,7,1)], AbsNormal)
        @test [Date(2013,7,1), Date(2014,7,1)] == findcontent(Vector{Date}, "DateVect", n, AbsNormal)
        updatecontent!([Date(2013,7,2), Date(2014,7,1)], "DateVect", n, AbsNormal)
        @test [Date(2013,7,2), Date(2014,7,1)] == findcontent(Vector{Date}, "DateVect", n, AbsNormal)

        addnode!(n, "TimeVect", [Time(12,53,42), Time(12,53,40)], AbsNormal)
        @test [Time(12,53,42), Time(12,53,40)] == findcontent(Vector{Time}, "TimeVect", n, AbsNormal)
        updatecontent!([Time(12,53,43), Time(12,53,40)], "TimeVect", n, AbsNormal)
        @test [Time(12,53,43), Time(12,53,40)] == findcontent(Vector{Time}, "TimeVect", n, AbsNormal)

        addnode!(n, "AnyVect", ["aa", Time(12,53,40), 2, nothing], AbsNormal)
        @test string.(["aa", Time(12,53,40), 2]) == findcontent(typeof(["aa", Time(12,53,40), 2, nothing]), "AnyVect", n, AbsNormal)
        updatecontent!( ["aa", Time(12,53,40), 3, nothing], "AnyVect", n, AbsNormal)
        @test string.(["aa", Time(12,53,40), 3]) == findcontent(typeof(["aa", Time(12,53,40), 2, nothing]), "AnyVect", n, AbsNormal)

        addnode!(n, "Dict", Dict("class"=>"h1", "id"=>"h2"), AbsNormal)
        pprint(n)
        @test_throws MethodError findcontent(Dict("class"=>"h1", "id"=>"h2"), "Dict", n, AbsNormal)
        @test_broken updatecontent!(Dict("class"=>"h1", "id"=>"h2"), "Dict", n, AbsNormal)

        n = createnode(AbsNormal, "a")
        addnode!(n, "1", "val1", AbsText)
        @test "val1" == findcontent(String, "1", n, AbsText)
        updatecontent!("val2", "", n, AbsText)
        @test "val2" == findcontent(String, "", n, AbsText)

        n = createnode(AbsNormal, "a")
        addnode!(n, "1", 1, AbsText)
        @test 1 == findcontent(Int64, "1", n, AbsText)
        updatecontent!(2, "", n, AbsText)
        @test 2 == findcontent(Int64, "", n, AbsText)

        n = createnode(AbsNormal, "p")
        addnode!(n, "n1", "1", AbsNormal)
        addnode!(n, "n2", "2", AbsNormal)
        addnode!(n, "n3", "3", AbsNormal)
        addnode!(n, "n4", "4", AbsNormal)
        addnode!(n, "[2:3]", ["txt2", "txt3"], AbsText) # Adds text inside n1 and n2
        @test ["txt2", "txt3"] == findcontent(Vector{String}, "[2:3]", n, AbsText)
        updatecontent!(["txt2-u", "txt3-u"] , "[2:3]", n, AbsText)
        @test ["txt2-u", "txt3-u"] == findcontent(Vector{String}, "[2:3]", n, AbsText)
        # TODO will throw error if text nodes don't exist
        @test_broken updatecontent!(["txt1-uu", "txt2-uu"] , "[1:2]", n, AbsText)
        @test_broken ["txt1-uu", "txt2-uu"] == findcontent(Vector{String}, "[1:2]", n, AbsText)
        @test_broken ["txt1-uu", "txt2-uu", "txt3-u"] == findcontent(Vector{String}, "[1:2]", n, AbsText)

        n = createnode(AbsNormal, "p")
        addnode!(n, "n1", "1", AbsNormal)
        addnode!(n, "n2", "2", AbsNormal)
        addnode!(n, "n3", "3", AbsNormal)
        addnode!(n, "n4", "4", AbsNormal)
        # puts form 1 to end
        addnode!(n, "[:]", ["txt1", "txt2"], AbsText) # Adds text inside n1 and n2
        @test ["txt1", "txt2"] == findcontent(Vector{String}, "[:]", n, AbsText)
        updatecontent!(["txt1-u", "txt2"], "[:]", n, AbsText)
        @test ["txt1-u", "txt2"] == findcontent(Vector{String}, "[:]", n, AbsText)

        n = createnode(AbsNormal, "a")
        addnode!(n, "1", 1, AbsText)
        @test 1 == findcontent(Int64, "1", n, AbsText)
        updatecontent!(2, "", n, AbsText)
        @test 2 == findcontent(Int64, "", n, AbsText)

        @testset "Node link" begin
            # Linking two nodes
            n1 = createnode(AbsNormal, "n1")
            n2 = createnode(AbsNormal, "n2")
            pprint(n1)
            pprint(n2)
            addnode!(n1, "n2", n2, AbsNormal)
            pprint(n1)
            pprint(n2)

            # Linking a node and a type
            n3 = createnode(AbsNormal, "n3")
            @aml mutable struct n4 "~"
                a::UN{String} = nothing, "~"
            end
            n4i = n4(a="1")
            pprint(n3)
            pprint(n4i)
            addnode!(n3, "n2", n4i, AbsNormal)
            pprint(n3)
            pprint(n4i)

            # Linking a node and an empty type
            n5 = createnode(AbsNormal, "n5")
            @aml mutable struct n6 "~"
                a::UN{String} = nothing, "~"
            end
            n6i = n6()
            pprint(n5)
            pprint(n6i)
            addnode!(n5, "n2", n6i, AbsNormal)
            pprint(n5)
            pprint(n6i)
        end
    end

    @testset "Html Document" begin
        dhtml = createnode(AbsHtml, "html")

        addnode!(dhtml, "normalString", "val", AbsNormal)
        @test "val" == findcontent(String, "normalString", dhtml, AbsNormal)
        @test "val" == findfirst("normalString", dhtml, AbsNormal).content
        @test ["val"] == findcontent("normalString", dhtml, AbsNormal)

        updatecontent!("val2", "normalString", dhtml, AbsNormal)
        @test "val2" == findcontent(String, "normalString", dhtml, AbsNormal)

        addnode!(dhtml, "attString", "val", AbsAttribute)
        @test "val" == findcontent(String, "attString", dhtml, AbsAttribute)
        updatecontent!("val2", "attString", dhtml, AbsAttribute)
        @test "val2" == findcontent(String, "attString", dhtml, AbsAttribute)

        using Dates
        addnode!(dhtml, "Date", Date(2013,7,1), AbsNormal)
        @test Date(2013,7,1) == findcontent(Date, "Date", dhtml, AbsNormal)
        updatecontent!(Date(2013,7,2), "Date", dhtml, AbsNormal)
        @test Date(2013,7,2) == findcontent(Date, "Date", dhtml, AbsNormal)

        addnode!(dhtml, "Time", Time(12,53,40), AbsNormal)
        @test Time(12,53,40) == findcontent(Time, "Time", dhtml, AbsNormal)
        updatecontent!(Time(12,53,41), "Time", dhtml, AbsNormal)
        @test Time(12,53,41) == findcontent(Time, "Time", dhtml, AbsNormal)

        addnode!(dhtml, "DateTime", DateTime(2013,5,1,12,53,40), AbsNormal)
        @test DateTime(2013,5,1,12,53,40) == findcontent(DateTime, "DateTime", dhtml, AbsNormal)
        updatecontent!( DateTime(2013,5,1,12,53,41), "DateTime", dhtml, AbsNormal)
        @test DateTime(2013,5,1,12,53,41) == findcontent(DateTime, "DateTime", dhtml, AbsNormal)

        using DataFrames
        addnode!(dhtml, "DataFrame",  DataFrame(course = ["Artificial Intelligence", "Robotics"], professor = ["Prof. A", "Prof. B"] ), AbsNormal)
        @test_skip  DataFrame(course = ["Artificial Intelligence", "Robotics"], professor = ["Prof. A", "Prof. B"] ) ==  findcontent(DataFrame, "DataFrame", dhtml, AbsNormal)
        # updatecontent!( DataFrame(course = ["Artificial Intelligence", "Robotics2"], professor = ["Prof. A", "Prof. B"] ),"DataFrame", dhtml, AbsNormal)
        # @test_skip  DataFrame(course = ["Artificial Intelligence", "Robotics2"], professor = ["Prof. A", "Prof. B"] ) ==  findcontent(DataFrame, "DataFrame", dhtml, AbsNormal)

        addnode!(dhtml, "Nothing", nothing, AbsAttribute)
        @test nothing == findcontent(Nothing, "Nothing", dhtml, AbsAttribute)
        @test nothing == findcontent(Union{Nothing,String}, "Nothing", dhtml, AbsAttribute)

        # updatecontent!(nothing, "Nothing", dhtml, AbsAttribute)
        # @test nothing == findcontent(Nothing, "Nothing", dhtml, AbsAttribute)

        ################################################################
        addnode!(dhtml, "stringVect", ["aa", "bb"], AbsNormal)
        @test ["aa", "bb"] == findcontent(Vector{String}, "stringVect", dhtml, AbsNormal)
        updatecontent!(["aa2", "bb"], "stringVect", dhtml, AbsNormal)
        @test ["aa2", "bb"] == findcontent(Vector{String}, "stringVect", dhtml, AbsNormal)
        @test ["aa2", "bb"] == getproperty.(findall("stringVect", dhtml, AbsNormal), :content)

        addnode!(dhtml, "floatVect", [5.6, 7.8], AbsNormal)
        @test [5.6, 7.8] == findcontent(Vector{Float64}, "floatVect", dhtml, AbsNormal)
        updatecontent!([5.6, 7.9], "floatVect", dhtml, AbsNormal)
        @test [5.6, 7.9] == findcontent(Vector{Float64}, "floatVect", dhtml, AbsNormal)

        addnode!(dhtml, "intVect", [5, 6], AbsNormal)
        @test [5, 6] == findcontent(Vector{Int64}, "intVect", dhtml, AbsNormal)
        updatecontent!([5, 7], "intVect", dhtml, AbsNormal)
        @test [5, 7] == findcontent(Vector{Int64}, "intVect", dhtml, AbsNormal)

        addnode!(dhtml, "DateVect", [Date(2013,7,1), Date(2014,7,1)], AbsNormal)
        @test [Date(2013,7,1), Date(2014,7,1)] == findcontent(Vector{Date}, "DateVect", dhtml, AbsNormal)
        updatecontent!([Date(2013,7,2), Date(2014,7,1)], "DateVect", dhtml, AbsNormal)
        @test [Date(2013,7,2), Date(2014,7,1)] == findcontent(Vector{Date}, "DateVect", dhtml, AbsNormal)

        addnode!(dhtml, "TimeVect", [Time(12,53,42), Time(12,53,40)], AbsNormal)
        @test [Time(12,53,42), Time(12,53,40)] == findcontent(Vector{Time}, "TimeVect", dhtml, AbsNormal)
        updatecontent!([Time(12,53,43), Time(12,53,40)], "TimeVect", dhtml, AbsNormal)
        @test [Time(12,53,43), Time(12,53,40)] == findcontent(Vector{Time}, "TimeVect", dhtml, AbsNormal)

        addnode!(dhtml, "AnyVect", ["aa", Time(12,53,40), 2, nothing], AbsNormal)
        @test string.(["aa", Time(12,53,40), 2]) == findcontent(typeof(["aa", Time(12,53,40), 2, nothing]), "AnyVect", dhtml, AbsNormal)
        updatecontent!( ["aa", Time(12,53,40), 3, nothing], "AnyVect", dhtml, AbsNormal)
        @test string.(["aa", Time(12,53,40), 3]) == findcontent(typeof(["aa", Time(12,53,40), 2, nothing]), "AnyVect", dhtml, AbsNormal)

        @testset "HTML Document link" begin
            # Linking two nodes
            x1 = createnode(AbsHtml)
            nx1 = createnode(AbsNormal, "nx1")
            pprint(x1)
            pprint(nx1)
            addnode!(x1, "nx1", nx1, AbsNormal)
            pprint(x1)
            pprint(nx1)

            # Linking a node and a type
            x2 = createnode(AbsHtml)
            @aml mutable struct nx2 "~"
                a::UN{String} = nothing, "~"
            end
            nx2i = nx2(a="1")
            pprint(x2)
            pprint(nx2i)
            addnode!(x2, "nx2", nx2i, AbsNormal)
            pprint(x2)
            pprint(nx2i)

            # Linking a node and an empty type
            x3 = createnode(AbsHtml)
            @aml mutable struct nx3 "~"
                a::UN{String} = nothing, "~"
            end
            nx3i = nx3()
            pprint(x3)
            pprint(nx3i)
            addnode!(x3, "nx3", nx3i, AbsNormal)
            pprint(x3)
            pprint(nx3i)
        end

    end

    @testset "XML Document" begin
        dxml = createnode(AbsXml, "xml")

        import EzXML: setroot!
        setroot!(dxml, createnode(AbsNormal, "node"))

        addnode!(dxml, "normalString", "val", AbsNormal)
        @test "val" == findcontent(String, "normalString", dxml, AbsNormal)
        @test ["val"] == findcontent("normalString", dxml, AbsNormal)

        updatecontent!("val2", "normalString", dxml, AbsNormal)
        @test "val2" == findcontent(String, "normalString", dxml, AbsNormal)

        addnode!(dxml, "attString", "val", AbsAttribute)
        @test "val" == findcontent(String, "attString", dxml, AbsAttribute)
        updatecontent!("val2", "attString", dxml, AbsAttribute)
        @test "val2" == findcontent(String, "attString", dxml, AbsAttribute)

        using Dates
        addnode!(dxml, "Date", Date(2013,7,1), AbsNormal)
        @test Date(2013,7,1) == findcontent(Date, "Date", dxml, AbsNormal)
        updatecontent!(Date(2013,7,2), "Date", dxml, AbsNormal)
        @test Date(2013,7,2) == findcontent(Date, "Date", dxml, AbsNormal)

        addnode!(dxml, "Time", Time(12,53,40), AbsNormal)
        @test Time(12,53,40) == findcontent(Time, "Time", dxml, AbsNormal)
        updatecontent!(Time(12,53,41), "Time", dxml, AbsNormal)
        @test Time(12,53,41) == findcontent(Time, "Time", dxml, AbsNormal)

        addnode!(dxml, "DateTime", DateTime(2013,5,1,12,53,40), AbsNormal)
        @test DateTime(2013,5,1,12,53,40) == findcontent(DateTime, "DateTime", dxml, AbsNormal)
        updatecontent!( DateTime(2013,5,1,12,53,41), "DateTime", dxml, AbsNormal)
        @test DateTime(2013,5,1,12,53,41) == findcontent(DateTime, "DateTime", dxml, AbsNormal)

        using DataFrames
        addnode!(dxml, "DataFrame",  DataFrame(course = ["Artificial Intelligence", "Robotics"], professor = ["Prof. A", "Prof. B"] ), AbsNormal)
        @test_skip  DataFrame(course = ["Artificial Intelligence", "Robotics"], professor = ["Prof. A", "Prof. B"] ) ==  findcontent(DataFrame, "DataFrame", dxml, AbsNormal)
        # updatecontent!( DataFrame(course = ["Artificial Intelligence", "Robotics2"], professor = ["Prof. A", "Prof. B"] ),"DataFrame", dxml, AbsNormal)
        # @test_skip  DataFrame(course = ["Artificial Intelligence", "Robotics2"], professor = ["Prof. A", "Prof. B"] ) ==  findcontent(DataFrame, "DataFrame", dxml, AbsNormal)

        addnode!(dxml, "Nothing", nothing, AbsAttribute)
        @test nothing == findcontent(Nothing, "Nothing", dxml, AbsAttribute)
        @test nothing == findcontent(Union{Nothing,String}, "Nothing", dxml, AbsAttribute)

        # updatecontent!(nothing, "Nothing", dxml, AbsAttribute)
        # @test nothing == findcontent(Nothing, "Nothing", dxml, AbsAttribute)

        ################################################################
        addnode!(dxml, "stringVect", ["aa", "bb"], AbsNormal)
        @test ["aa", "bb"] == findcontent(Vector{String}, "stringVect", dxml, AbsNormal)
        updatecontent!(["aa2", "bb"], "stringVect", dxml, AbsNormal)
        @test ["aa2", "bb"] == findcontent(Vector{String}, "stringVect", dxml, AbsNormal)

        addnode!(dxml, "floatVect", [5.6, 7.8], AbsNormal)
        @test [5.6, 7.8] == findcontent(Vector{Float64}, "floatVect", dxml, AbsNormal)
        updatecontent!([5.6, 7.9], "floatVect", dxml, AbsNormal)
        @test [5.6, 7.9] == findcontent(Vector{Float64}, "floatVect", dxml, AbsNormal)

        addnode!(dxml, "intVect", [5, 6], AbsNormal)
        @test [5, 6] == findcontent(Vector{Int64}, "intVect", dxml, AbsNormal)
        updatecontent!([5, 7], "intVect", dxml, AbsNormal)
        @test [5, 7] == findcontent(Vector{Int64}, "intVect", dxml, AbsNormal)

        addnode!(dxml, "DateVect", [Date(2013,7,1), Date(2014,7,1)], AbsNormal)
        @test [Date(2013,7,1), Date(2014,7,1)] == findcontent(Vector{Date}, "DateVect", dxml, AbsNormal)
        updatecontent!([Date(2013,7,2), Date(2014,7,1)], "DateVect", dxml, AbsNormal)
        @test [Date(2013,7,2), Date(2014,7,1)] == findcontent(Vector{Date}, "DateVect", dxml, AbsNormal)

        addnode!(dxml, "TimeVect", [Time(12,53,42), Time(12,53,40)], AbsNormal)
        @test [Time(12,53,42), Time(12,53,40)] == findcontent(Vector{Time}, "TimeVect", dxml, AbsNormal)
        updatecontent!([Time(12,53,43), Time(12,53,40)], "TimeVect", dxml, AbsNormal)
        @test [Time(12,53,43), Time(12,53,40)] == findcontent(Vector{Time}, "TimeVect", dxml, AbsNormal)

        addnode!(dxml, "AnyVect", ["aa", Time(12,53,40), 2, nothing], AbsNormal)
        @test string.(["aa", Time(12,53,40), 2]) == findcontent(typeof(["aa", Time(12,53,40), 2, nothing]), "AnyVect", dxml, AbsNormal)
        updatecontent!( ["aa", Time(12,53,40), 3, nothing], "AnyVect", dxml, AbsNormal)
        @test string.(["aa", Time(12,53,40), 3]) == findcontent(typeof(["aa", Time(12,53,40), 2, nothing]), "AnyVect", dxml, AbsNormal)

        @testset "XML Document link" begin
            # Linking two nodes
            x1 = createnode(AbsXml)
            nx1 = createnode(AbsNormal, "nx1")
            pprint(x1)
            pprint(nx1)
            addnode!(x1, "nx1", nx1, AbsNormal)
            pprint(x1)
            pprint(nx1)

            # Linking a node and a type
            x2 = createnode(AbsXml)
            @aml mutable struct nx2 "~"
                a::UN{String} = nothing, "~"
            end
            nx2i = nx2(a="1")
            pprint(x2)
            pprint(nx2i)
            addnode!(x2, "nx2", nx2i, AbsNormal)
            pprint(x2)
            pprint(nx2i)

            # Linking a node and an empty type
            x3 = createnode(AbsXml)
            @aml mutable struct nx3 "~"
                a::UN{String} = nothing, "~"
            end
            nx3i = nx3()
            pprint(x3)
            pprint(nx3i)
            addnode!(x3, "nx3", nx3i, AbsNormal)
            pprint(x3)
            pprint(nx3i)
        end
    end

end


text:
################################################################################
#
#  Simplify
#
################################################################################

function simplify(K::NfRel; cached::Bool = true, prec::Int = 100)
  OK = maximal_order(K)
  B = lll_basis(OK)
  B1 = _sieve_primitive_elements(B)
  i = 6
  while isempty(B1)
    B1 = _sieve_primitive_elements(B, parameter = i)
    i += 3
  end
  a = B1[1]
  I = t2(a)
  for i = 2:min(50, length(B1))
    J = t2(B1[i])
    if J < I
      a = B1[i]
      I = J
    end
  end
  f = minpoly(a)
  @assert degree(f) == degree(K)
  Ks = number_field(f, cached = cached, check = false)[1]
  mKs = hom(Ks, K, a)
  return Ks, mKs
end

function _sieve_primitive_elements(B::Vector{T}) where T <: NumFieldElem
  K = parent(B[1])
  n = absolute_degree(K)
  B1 = typeof(B)()
  k = div(n, 2)
  for x in B
    c = conjugates_arb(x, 16)
    isprimitive = true
    for i = 2:k+1
      for j = 1:i-1
        if overlaps(c[i], c[j])
          isprimitive = false
          break
        end
      end
      if !isprimitive
        break
      end
    end
    if isprimitive
      push!(B1, x)
    end
  end
  return B1
end

function _is_primitive_via_block(a::NfRelElem{nf_elem}, rt::Dict{fq, Vector{fq}}, Fx, tmp::gfp_fmpz_poly)
  if iszero(a)
    return false
  end
  n = degree(parent(a))
  pol = data(a)
  conjs = Set{fq}()
  for (r, vr) in rt
    coeffs = Vector{fq}(undef, degree(pol)+1)
    for i = 0:degree(pol)
      nf_elem_to_gfp_fmpz_poly!(tmp, coeff(pol, i))
      coeffs[i+1] = evaluate(tmp, r)
    end
    g = Fx(coeffs)
    for i = 1:length(vr)
      ev = evaluate(g, vr[i])
      if ev in conjs
        return false
      end
      push!(conjs, ev)
    end
  end
  return true
end

function _find_prime(L::NfRel{nf_elem})
  p = 2^10
  K = base_field(L)
  OK = maximal_order(K)
  OL = maximal_order(L)

  n_attempts = max(5, min(degree(L), 10))
  candidates = Vector{Tuple{Int, Int}}(undef, n_attempts)
  i = 1
  f = L.pol
  threshold = degree(f)^2
  den = lcm(fmpz[denominator(coeff(f, i)) for i = 0:degree(f)])
  while i < n_attempts+1
    p = next_prime(p)
    if isindex_divisor(OK, p) || divisible(absolute_discriminant(OL), p) || divisible(den, p)
      continue
    end
    lp = prime_decomposition(OK, p)
    P = lp[1][1]
    F, mF = ResidueField(OK, P)
    mF1 = extend(mF, K)
    fF = map_coefficients(mF1, f)
    if degree(fF) != degree(f) || !issquarefree(fF)
      continue
    end
    FS = factor_shape(fF)
    d = lcm(Int[x for (x, v) in FS])*degree(P)
    acceptable = true
    for j = 2:length(lp)
      Q = lp[j][1]
      F2, mF2 = ResidueField(OK, Q)
      mF3 = extend(mF2, K)
      fF2 = map_coefficients(mF3, f)
      if degree(fF2) != degree(f) || !issquarefree(fF2)
        acceptable = false
        break
      end
      FS = factor_shape(fF2)
      d1 = lcm(Int[x for (x, v) in FS])
      d = lcm(d, d1*degree(Q))
    end
    if acceptable && d < threshold
      candidates[i] = (p, d)
      i += 1
    end
  end

  res = candidates[1]
  for j = 2:n_attempts
    if candidates[j][2] < res[2]
      res = candidates[j]
    end
  end
  return res[1], res[2]
end

function _setup_block_system(Lrel::NfRel{nf_elem})
  K = base_field(Lrel)
  OK = maximal_order(K)
  Zx = ZZ["x"][1]
  n = absolute_degree(Lrel)

  pint, d = _find_prime(Lrel)
  p = fmpz(pint)
  abs_deg = d
  #First, we search for elements that are primitive using block systems
  Fp = GF(p, cached = false)
  Fpx = PolynomialRing(Fp, cached = false)[1]
  F = FlintFiniteField(p, abs_deg, "w", cached = false)[1]
  Fx = PolynomialRing(F, cached = false)[1]
  rt_base_field = roots(Zx(K.pol), F)
  tmp = Fpx()
  g = Lrel.pol
  rt = Dict{fq, Vector{fq}}()
  nroots = 0
  roots_needed = div(n, 2)+1
  for r in rt_base_field
    coeff_gF = fq[]
    for i = 0:degree(g)
      nf_elem_to_gfp_fmpz_poly!(tmp, coeff(g, i))
      push!(coeff_gF, evaluate(tmp, r))
    end
    gF = Fx(coeff_gF)
    rt[r] = roots(gF)
    nroots += length(roots(gF))
    if nroots >= roots_needed
      break
    end
  end
  return rt, Fx, tmp
end

function _find_prime(L::NfRelNS{nf_elem})
  p = 2^10
  K = base_field(L)
  OK = maximal_order(K)
  OL = maximal_order(L)
  dL = numerator(discriminant(L, FlintQQ))

  n_attempts = min(degree(L), 10)
  candidates = Vector{Tuple{NfOrdIdl, Int}}(undef, n_attempts)
  i = 1
  pols = L.pol
  threshold = absolute_degree(L)^2
  polsR = Vector{fq_poly}(undef, length(pols))
  while i < n_attempts+1
    p = next_prime(p)
    if isindex_divisor(OK, p) || divisible(dL, p)
      continue
    end
    lp = prime_decomposition(OK, p)
    P = lp[1][1]
    @assert !isindex_divisor(OL, P)
    F, mF = ResidueField(OK, P)
    mF1 = extend(mF, K)
    Fx, _ = PolynomialRing(F, "x", cached = false)
    is_proj = true
    for j = 1:length(pols)
      fF = to_univariate(Fx, map_coefficients(mF1, pols[j]))
      if degree(fF) != total_degree(pols[j]) || !issquarefree(fF)
        is_proj = false
        break
      end
      polsR[j] = fF
    end
    if !is_proj
      continue
    end
    d = 1
    for j = 1:length(polsR)
      FS = factor_shape(polsR[j])
      d1 = lcm(Int[x for (x, v) in FS])
      d = lcm(d, d1)
    end
    acceptable = true
    for s = 2:length(lp)
      Q = lp[s][1]
      @assert !isindex_divisor(OL, Q)
      F, mF = ResidueField(OK, Q)
      Fx, _ = PolynomialRing(F, "x", cached = false)
      mF1 = extend(mF, K)
      is_proj = true
      for j = 1:length(pols)
        fF = to_univariate(Fx, map_coefficients(mF1, pols[j]))
        if degree(fF) != total_degree(pols[j]) || !issquarefree(fF)
          is_proj = false
          break
        end
        polsR[j] = fF
      end
      if !is_proj
        acceptable = false
        break
      end
      for j = 1:length(polsR)
        FS = factor_shape(polsR[j])
        d1 = lcm(Int[x for (x, v) in FS])*degree(Q)
        d = lcm(d, d1)
      end
    end
    if !acceptable
      continue
    end
    if d*degree(P) < threshold
      candidates[i] = (P, d)
      i += 1
    end
  end
  res = candidates[1]
  for j = 2:n_attempts
    if candidates[j][2]*degree(candidates[j][1]) < res[2]*degree(res[1])
      res = candidates[j]
    end
  end
  return res[1], res[2]
end

function _setup_block_system(Lrel::NfRelNS{nf_elem})
  K = base_field(Lrel)
  OK = maximal_order(K)
  Zx = ZZ["x"][1]

  n = absolute_degree(Lrel)

  P, d = _find_prime(Lrel)
  p = minimum(P, copy = false)
  abs_deg = d*degree(P)
  #First, we search for elements that are primitive using block systems
  Fp = GF(p, cached = false)
  Fpx = PolynomialRing(Fp, cached = false)[1]
  F = FlintFiniteField(p, abs_deg, "w", cached = false)[1]
  Fx = PolynomialRing(F, cached = false)[1]
  rt_base_field = roots(Zx(K.pol), F)
  rt = Dict{fq, Vector{Vector{fq}}}()
  Rxy = PolynomialRing(F, ngens(Lrel), cached = false)[1]
  tmp = Fpx()
  Kx, _ = PolynomialRing(K, "x", cached = false)
  for r in rt_base_field
    vr = Vector{Vector{fq}}()
    for f in Lrel.pol
      g = to_univariate(Kx, f)
      coeff_gF = fq[]
      for i = 0:degree(g)
        nf_elem_to_gfp_fmpz_poly!(tmp, coeff(g, i))
        push!(coeff_gF, evaluate(tmp, r))
      end
      gF = Fx(coeff_gF)
      push!(vr, roots(gF))
    end
    rt[r] = vr
  end
  rt1 = Dict{fq, Vector{Vector{fq}}}()
  ind = 1
  nconjs_needed = div(n, 2)+1
  for (r, v) in rt
    rtv = Vector{Vector{fq}}()
    it = cartesian_product_iterator([1:length(v[i]) for i in 1:length(v)], inplace = true)
    for i in it
      push!(rtv, [v[j][i[j]] for j = 1:length(v)])
      ind += 1
      if ind > nconjs_needed
        break
      end
    end
    rt1[r] = rtv
    if ind > nconjs_needed
      break
    end
  end
  @assert ind > nconjs_needed
  return rt1, Rxy, tmp
end

function _sieve_primitive_elements(B::Vector{T}; parameter::Int = div(absolute_degree(parent(B[1])), 2)) where T <: Union{NfRelNSElem{nf_elem}, NfRelElem{nf_elem}}
  Lrel = parent(B[1])
  #First, we choose the candidates
  ape = absolute_primitive_element(Lrel)
  B_test = vcat(B, T[ape])
  Bnew = typeof(B)()
  nrep = min(parameter, absolute_degree(Lrel))
  for i = 1:length(B_test)
    push!(Bnew, B_test[i])
    for j = 1:nrep
      if i != j
        push!(Bnew, B_test[i]+B_test[j])
        push!(Bnew, B_test[i]-B_test[j])
      end
    end
  end
  #Now, we test for primitiveness.
  rt1, Rxy, tmp = _setup_block_system(Lrel)
  indices = Int[]
  for i = 1:length(Bnew)
    if _is_primitive_via_block(Bnew[i], rt1, Rxy, tmp)
      push!(indices, i)
    end
  end
  return Bnew[indices]
end

function _is_primitive_via_block(a::NfRelNSElem{nf_elem}, rt::Dict{fq, Vector{Vector{fq}}}, Rxy, tmp)
  if length(vars(a.data)) < ngens(parent(a))
    return false
  end
  n = degree(parent(a))
  pol = data(a)
  conjs = Set{fq}()
  for (r, vr) in rt
    ctx = MPolyBuildCtx(Rxy)
    for (c, v) in zip(coefficients(pol), exponent_vectors(pol))
      nf_elem_to_gfp_fmpz_poly!(tmp, c)
      push_term!(ctx, evaluate(tmp, r), v)
    end
    g = finish(ctx)
    for i = 1:length(vr)
      ev = evaluate(g, vr[i])
      if ev in conjs
        return false
      end
      push!(conjs, ev)
    end
  end
  return true
end

function _find_short_primitive_element(L::NfRelNS)
  B = lll_basis(maximal_order(L))
  parameter = div(absolute_degree(L), 2)
  B1 = _sieve_primitive_elements(B, parameter = parameter)
  while isempty(B1)
    parameter += 1
    B1 = _sieve_primitive_elements(B, parameter = parameter)
  end
  a = B1[1]
  I = t2(a)
  for i = 2:length(B1)
    J = t2(B1[i])
    if J < I
      a = B1[i]
      I = J
    end
  end
  return a
end

function simplified_absolute_field(L::NfRelNS; cached = false)
  a = _find_short_primitive_element(L)
  f = absolute_minpoly(a)
  @assert degree(f) == absolute_degree(L)
  K = number_field(f, check = false, cached = cached)[1]
  mp = hom(K, L, a)
  return K, mp
end

function simplified_absolute_field(L::NfRel; cached::Bool = false)
  OL = maximal_order(L)
  B = lll_basis(OL)
  B1 = _sieve_primitive_elements(B)
  nrep = 4
  while isempty(B1)
    nrep += 1
    B1 = _sieve_primitive_elements(B, parameter = nrep)
  end
  a = B1[1]
  I = t2(a)
  for i = 2:min(50, length(B1))
    J = t2(B1[i])
    if J < I
      a = B1[i]
      I = J
    end
  end
  f = absolute_minpoly(a)
  @assert degree(f) == absolute_degree(L)
  K = number_field(f, check = false, cached = cached)[1]
  mp = hom(K, L, a)
  imp = inv(mp)
  return K, mp, hom(base_field(L), K,  imp(L(gen(base_field(L)))))
end


text:
using Beamforming
using Test

@testset "All tests" begin
    include("stack.jl")
    include("beamforming.jl")
    include("vespa.jl")
    include("synth.jl")
    include("arf.jl")
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming