text:
# fast robust geometric predicates
# Converted into Julia from C
# Original functions written by
#  Jonathan Richard Shewchuk
#  School of Computer Science
#  Carnegie Mellon University
#  5000 Forbes Avenue
#  Pittsburgh, Pennsylvania  15213-3891
#  jrs@cs.cmu.edu

using GeometryTypes

"""
original comment:
/* Which of the following two methods of finding the absolute values is      */
/*   fastest is compiler-dependent.  A few compilers can inline and optimize */
/*   the fabs() call; but most will incur the overhead of a function call,   */
/*   which is disastrously slow.  A faster way on IEEE machines might be to  */
/*   mask the appropriate bit, but that's difficult to do in C.              */
This function definition attempts to do the same thing as a C-Define
"""
@inline absolute(a) = ((a) >= 0.0 ? (a) : -(a))
# TODO: should also try the built in absolute
# @inline absolute(a) = abs(a)

# Skipping tails, condensing to one function (?)
#@inline function fast_two_sum_tail(a::Float64, b::Float64, x::Float64)::Tuple{Float64,Float64}
#    bvirt::Float64 = x - a
#    x, b - bvirt
#end

@inline function fast_two_sum(a::Float64, b::Float64)::Tuple{Float64,Float64}
    x::Float64 = a + b
    bvirt::Float64 = x - a
    x, b - bvirt
end

#@inline function fast_two_diff_tail(a::Float64, b::Float64, x::Float64)

const resulterrbound = Ref{Float64}()
const ccwerrboundA = Ref{Float64}()
const ccwerrboundB = Ref{Float64}()
const ccwerrboundC = Ref{Float64}()
const o3derrboundA = Ref{Float64}()
const o3derrboundB = Ref{Float64}()
const o3derrboundC = Ref{Float64}()
const iccerrboundA = Ref{Float64}()
const iccerrboundB = Ref{Float64}()
const iccerrboundC = Ref{Float64}()
const isperrboundA = Ref{Float64}()
const isperrboundB = Ref{Float64}()
const isperrboundC = Ref{Float64}()
"""
/*****************************************************************************/
/*                                                                           */
/*  exactinit()   Initialize the variables used for exact arithmetic.        */
/*                                                                           */
/*  `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in   */
/*  floating-point arithmetic.  `epsilon' bounds the relative roundoff       */
/*  error.  It is used for floating-point error analysis.                    */
/*                                                                           */
/*  `splitter' is used to split floating-point numbers into two half-        */
/*  length significands for exact multiplication.                            */
/*                                                                           */
/*  I imagine that a highly optimizing compiler might be too smart for its   */
/*  own good, and somehow cause this routine to fail, if it pretends that    */
/*  floating-point arithmetic is too much like real arithmetic.              */
/*                                                                           */
/*  Don't change this routine unless you fully understand it.                */
/*                                                                           */
/*****************************************************************************/
"""
function exactinit()::Nothing
    every_other::Bool = false
    epsilon::Float64 = 1.0
    splitter::Float64 = 1.0
    check::Float64 = 1.0
    #/* Repeatedly divide `epsilon' by two until it is too small to add to    */
    #/*   one without causing roundoff.  (Also check if the sum is equal to   */
    #/*   the previous sum, for machines that round up instead of using exact */
    #/*   rounding.  Not that this library will work on such machines anyway. */
    lastcheck = check;
    epsilon *= 0.5
    splitter *= 2.0
    check = 1.0 + epsilon
    while ((check != 1.0) && (check != lastcheck))
        lastcheck = check;
        epsilon *= 0.5
        if (every_other)
            splitter *= 2.0
        end
        every_other = !every_other
        check = 1.0 + epsilon
    end
    splitter += 1.0

    #/* Error bounds for orientation and incircle tests. */
    global resulterrbound[] = (3.0 + 8.0 * epsilon) * epsilon
    global ccwerrboundA[] = (3.0 + 16.0 * epsilon) * epsilon
    global ccwerrboundB[] = (2.0 + 12.0 * epsilon) * epsilon
    global ccwerrboundC[] = (9.0 + 64.0 * epsilon) * epsilon * epsilon
    global o3derrboundA[] = (7.0 + 56.0 * epsilon) * epsilon
    global o3derrboundB[] = (3.0 + 28.0 * epsilon) * epsilon
    global o3derrboundC[] = (26.0 + 288.0 * epsilon) * epsilon * epsilon
    global iccerrboundA[] = (10.0 + 96.0 * epsilon) * epsilon
    global iccerrboundB[] = (4.0 + 48.0 * epsilon) * epsilon
    global iccerrboundC[] = (44.0 + 576.0 * epsilon) * epsilon * epsilon
    global isperrboundA[] = (16.0 + 224.0 * epsilon) * epsilon
    global isperrboundB[] = (5.0 + 72.0 * epsilon) * epsilon
    global isperrboundC[] = (71.0 + 1408.0 * epsilon) * epsilon * epsilon
    nothing
end

"""
/*               Return a positive value if the point pd lies below the      */
/*               plane passing through pa, pb, and pc; "below" is defined so */
/*               that pa, pb, and pc appear in counterclockwise order when   */
/*               viewed from above the plane.  Returns a negative value if   */
/*               pd lies above the plane.  Returns zero if the points are    */
/*               coplanar.  The result is also a rough approximation of six  */
/*               times the signed volume of the tetrahedron defined by the   */
/*               four points.                                                */
"""
function orient3d(pa::Point3f0, pb::Point3f0, pc::Point3f0, pd::Point3f0)::Float64
    ad::Point3f0 = pa .- pd
    bd::Point3f0 = pb .- pd
    cd::Point3f0 = pc .- pd

    bdxcdy::Float64 = bd[1] * cd[2]
    cdxbdy::Float64 = cd[1] * bd[2]
    cdxady::Float64 = cd[1] * ad[2]
    adxcdy::Float64 = ad[1] * cd[2]
    adxbdy::Float64 = ad[1] * bd[2]
    bdxady::Float64 = bd[1] * ad[2]

    det::Float64 = adz * (bdxcdy - cdxbdy) +
                   bdz * (cdxady - adxcdy) +
                   cdz * (adxbdy - bdxady)
    permanent = (absolute(bdxcdy) + absolute(cdxbdy)) * absolute(adz) +
                (absolute(cdxady) + absolute(adxcdy)) * absolute(bdz) +
                (absolute(adxbdy) + absolute(bdxady)) * absolute(cdz)
    # TODO: define predetermined error bounds
    errbound = o3derrboundA * permanent
    if (det > errbound) || (-det > errbound)
        return det
    else
        println("WARNING: possible errors in floating point arithmetic!!!")
        return det
    end
    # TODO: implement robust predicates
    #return orient3dadapt(pa, pb, pc, pd, permanent)
end

"""
/*               Return a positive value if the point pe lies inside the     */
/*               sphere passing through pa, pb, pc, and pd; a negative value */
/*               if it lies outside; and zero if the five points are         */
/*               cospherical.  The points pa, pb, pc, and pd must be ordered */
/*               so that they have a positive orientation (as defined by     */
/*               orient3d()), or the sign of the result will be reversed.    */
"""
function insphere(pa::Point3f0, pb::Point3f0, pc::Point3f0, pd::Point3f0, pe::Point3f0)::Float64
    ae::Point3f0 = pa .- pe
    be::Point3f0 = pb .- pe
    ce::Point3f0 = pc .- pe
    de::Point3f0 = pd .- pe

    aexbey = ae[1] * be[2]
    bexaey = be[1] * ae[2]
    ab = aexbey - bexaey
    bexcey = be[1] * ce[2]
    cexbey = ce[1] * be[2]
    bc = bexcey - cexbey
    cexdey = ce[1] * de[2]
    dexcey = de[1] * ce[2]
    cd = cexdey - dexcey
    dexaey = de[1] * ae[2]
    aexdey = ae[1] * de[2]
    da = dexaey - aexdey

    aexcey = ae[1] * ce[2]
    cexaey = ce[1] * ae[2]
    ac = aexcey - cexaey
    bexdey = be[1] * de[2]
    dexbey = de[1] * be[2]
    bd = bexdey - dexbey

    abc = ae[3] * bc - be[3] * ac + ce[3] * ab
    bcd = be[3] * cd - ce[3] * bd + de[3] * bc
    cda = ce[3] * da + de[3] * ac + ae[3] * cd
    dab = de[3] * ab + ae[3] * bd + be[3] * da

    alift = sum(ae .* ae)
    blift = sum(be .* be)
    clift = sum(ce .* ce)
    dlift = sum(de .* de)

    det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd)
    permanent = ((absolute(cexdey) + absolute(dexcey)) * absolute(bez) +
                 (absolute(dexbey) + absolute(bexdey)) * absolute(cez) +
                 (absolute(bexcey) + absolute(cexbey)) * absolute(dez)) * alift +
                ((absolute(dexaey) + absolute(aexdey)) * absolute(cez) +
                 (absolute(aexcey) + absolute(cexaey)) * absolute(dez) +
                 (absolute(cexdey) + absolute(dexcey)) * absolute(aez)) * blift +
                ((absolute(aexbey) + absolute(bexaey)) * absolute(dez) +
                 (absolute(bexdey) + absolute(dexbey)) * absolute(aez) +
                 (absolute(dexaey) + absolute(aexdey)) * absolute(bez)) * clift +
                ((absolute(bexcey) + absolute(cexbey)) * absolute(aez) +
                 (absolute(cexaey) + absolute(aexcey)) * absolute(bez) +
                 (absolute(aexbey) + absolute(bexaey)) * absolute(cez)) * dlift
    errbound = isperrboundA * permanent;
    if (det > errbound) || (-det > errbound)
        return det
    else
        println("WARNING: possible error in floating point arithmetic!!!")
        return det
    end
    # TODO: implement adaptive predicates
end


text:
__precompile__(true)

module SortMerge

using Printf, SparseArrays, ProgressMeter, StatsBase

import Base.show
import Base.sortperm
import Base.zip

import Base.length, Base.size, Base.getindex,
       Base.firstindex, Base.lastindex, Base.IndexStyle

export sortmerge, sortperm, nmatch, countmatch, multimatch, simple_join

# --------------------------------------------------------------------
# Source structure
#
struct Source
    size::Int
    sortperm::Vector{Int}
    cmatch::SparseVector{Int,Int}
end

# --------------------------------------------------------------------
# Matched structure
#
struct Matched <: AbstractVector{Vector{Int}}
    nsrc::Int
    nrow::Int
    sources::Vector{Source}
    matched::Matrix{Int}
end

# Methods
length(mm::Matched) = mm.nsrc
size(mm::Matched) = (mm.nsrc,)
getindex(mm::Matched, inds) = mm.matched[:, inds]
firstindex(mm::Matched) = 1
lastindex(mm::Matched) = mm.nsrc
IndexStyle(::Type{Matched}) = IndexLinear()

nmatch(mm::Matched) = mm.nrow
countmatch(mm::Matched, source::Int) = mm.sources[source].cmatch
sortperm(  mm::Matched, source::Int) = mm.sources[source].sortperm

zip(mm::Matched) = zip(map(i -> mm.matched[:,i], 1:mm.nsrc)...)

function subset(mm::Matched, selected::Vector{Int})
    match = fill(0, length(selected), mm.nsrc)
    sources = Vector{Source}()
    for i in 1:mm.nsrc
        match[:, i] = mm[i][selected]
        cm = countmap(match[:, i])
        ii = collect(keys(cm))
        cc = collect(values(cm))
        push!(sources, Source(mm.sources[i].size, mm.sources[i].sortperm,
                              sparsevec(ii, cc, mm.sources[i].size)))
    end
    return Matched(mm.nsrc, length(selected), sources, match)
end

function multimatch(mm::Matched, source::Int, multi::Int; group=false)
    @assert multi >= 1
    index = findall(countmatch(mm, source) .== multi)
    index_out = sortmerge(index, mm.matched[:,source])[2]
    matrix = mm.matched[index_out,:]

    out = Vector{Matched}()
    ngroups = (group  ?  length(index)  :  1)
    for igroup in 1:ngroups
        jj = (group  ?
              findall(matrix[:,source] .== index[igroup])  :
              collect(1:length(index_out)))

        sources = Vector{Source}()
        for i in 1:mm.nsrc
            push!(sources, Source(mm.sources[i].size, Vector{Int}(),
                                  sparsevec(unique(mm.matched[index_out[jj], i]),
                                            multi, length(mm.sources[i].cmatch))))
        end
        push!(out, Matched(mm.nsrc, length(jj), sources, mm.matched[index_out[jj],:]))
    end

    (group)  &&  (return out)
    return out[1]
end

default_lt(v, i, j) = (v[i] < v[j])
default_sd(A, B, i, j) = sign(A[i] - B[j])

# sortmerge(j::NTuple{2, Matched}, A, B, args...; sd=default_sd) =
#   sortmerge(A, B, args..., sort1=j[1].sortperm, sort2=j[2].sortperm, sd=sd)
function sortmerge(A, B, args...;
                   sd=default_sd,
                   sort1=Vector{Int}(),
                   sort2=Vector{Int}(),
                   lt1=default_lt,
                   lt2=default_lt,
                   sorted=false)

    size1 = size(A)[1]
    size2 = size(B)[1]

    if length(sort1) == 0
        sort1 = collect(1:size1)
        if !sorted
            sort1 = sortperm(sort1, lt=(i, j) -> (lt1(A, i, j)))
        end
    end
    if length(sort2) == 0
        sort2 = collect(1:size2)
        if !sorted
            sort2 = sortperm(sort2, lt=(i, j) -> (lt2(B, i, j)))
        end
    end

    i2a = 1
    match1 = Array{Int}(undef, 0)
    match2 = Array{Int}(undef, 0)
    cm1 = fill(0, size1)
    cm2 = fill(0, size2)

    desc = "SortMerge "
    barlen = ProgressMeter.tty_width(desc, stderr, false)
    (barlen > 50)  &&  (barlen = 50)
    prog = Progress(size1, desc=desc, dt=0.5, color=:light_black, barlen=barlen,
                    barglyphs=BarGlyphs('|','█', ['▏','▎','▍','▌','▋','▊','▉'],' ','|',))
    for i1 in 1:size1
        update!(prog, i1)
        for i2 in i2a:size2
            j1 = sort1[i1]
            j2 = sort2[i2]
            if length(args) > 0  # This improves performances
                dd = Int(sd(A, B, j1, j2, args...))
            else
                dd = Int(sd(A, B, j1, j2))
            end

            if     dd == -1; break
            elseif dd ==  1; i2a += 1
            elseif dd ==  0
                push!(match1, j1)
                push!(match2, j2)
                cm1[j1] += 1
                cm2[j2] += 1
            end
        end
    end
    finish!(prog)

    ii = findall(cm1 .> 0); side1 = Source(size1, sort1, sparsevec(ii, cm1[ii], length(cm1)))
    ii = findall(cm2 .> 0); side2 = Source(size2, sort2, sparsevec(ii, cm2[ii], length(cm2)))
    mm = Matched(2, length(match1), [side1, side2], [match1 match2])
    return mm
end

show(io::IO, mime::MIME"text/plain", mm::Matched) = show(io, mm)
show(mm::Matched) = show(stdout, mm)
function show(io::IO, mm::Matched)
    for i in 1:length(mm.sources)
        ss = mm.sources[i]
        uu = length(unique(mm.matched[:,i]))
        maxmult = 0
        (length(ss.cmatch) > 0)  &&  (maxmult = maximum(ss.cmatch))
        @printf(io, "Input %1d: %12d / %12d  (%6.2f%%)  -  max mult. %d\n",
                i, uu, ss.size, 100. * uu/float(ss.size), maxmult)
    end
    @printf(io, "Output : %12d\n", size(mm.matched)[1])
    nothing
end

simple_join(A, B, match::Function) =
    sortmerge(A, B, sorted=true,
              sd=(A, B, i, j) -> (match(A[i], B[j])  ?  0  :  999))

end # module

ext: jl


text:
function jump_to_cbf(problem, name, filename)
    # these are internal functions
    c, A, b, var_cones, con_cones = eval(Main, :(JuMP.conicdata($problem)))
    vartypes = eval(Main, :(JuMP.vartypes_without_fixed($problem)))

    dat = mpbtocbf(name, c, A, b, con_cones, var_cones, vartypes, :Min)

    writecbfdata(filename, dat)
end


text:
module MsgPack

const INT_FP   = (0x00, 0x7f)
const MAP_F    = (0x80, 0x8f)
const ARR_F    = (0x90, 0x9f)
const STR_F    = (0xa0, 0xbf)

const NIL      = 0xc0
const UNUSED   = 0xc1
const FALSE    = 0xc2
const TRUE     = 0xc3
const BIN_8    = 0xc4
const BIN_16   = 0xc5
const BIN_32   = 0xc6
const EXT_8    = 0xc7
const EXT_16   = 0xc8
const EXT_32   = 0xc9
const FLOAT_32 = 0xca
const FLOAT_64 = 0xcb
const UINT_8   = 0xcc
const UINT_16  = 0xcd
const UINT_32  = 0xce
const UINT_64  = 0xcf
const INT_8    = 0xd0
const INT_16   = 0xd1
const INT_32   = 0xd2
const INT_64   = 0xd3

const EXT_F    = (0xd4, 0xd8)

const STR_8    = 0xd9
const STR_16   = 0xda
const STR_32   = 0xdb
const ARR_16   = 0xdc
const ARR_32   = 0xdd
const MAP_16   = 0xde
const MAP_32   = 0xdf

const INT_FN   = (0xe0, 0xff)

const TYPECODE_TO_TYPE = Dict{Int8, Any}()

# For custom encoding.
function encode(x::T)::Vector{UInt8} where T
    tmp = [getfield(x, name) for name in fieldnames(T)]
    MsgPack.pack(tmp)
end

function decode(::Type{T}, x::Vector{UInt8})::T where T
    args = MsgPack.unpack(x)
    T(args...)
end

"""
- Use function overloadings instead of dict.
- Because we want to return the typecode of the parent type if self is not defined.
- E.g., in parameteric types `A{T}`.
"""
get_typecode(T::DataType) = error("Type $T is not registered with a typecode")

function register(::Type{T}, typecode::Integer) where T
    if haskey(TYPECODE_TO_TYPE, typecode)
        error("Type Code $typecode was already registered")
    end
    @eval get_typecode(x::$T)::Int8 = $typecode
    TYPECODE_TO_TYPE[typecode] = T
end

struct Ext
    typecode::Int8
    data::Vector{UInt8}

    function Ext(t::Integer, d::Vector{UInt8}; impltype::Bool=false)
        # -128 to -1 reserved for implementation
        if -128 <= t <= -1
            impltype || error("MsgPack Ext typecode -128 through -1 reserved by implementation")
        elseif !(0 <= t <= 127)
            error("MsgPack Ext typecode must be in the range [-128, 127]")
        end

        new(t, d)
    end
end

Base.:(==)(a::Ext, b::Ext) = a.typecode == b.typecode && a.data == b.data

readn(s::IO, t::DataType) = ntoh(read(s, t))
readi(s::IO, t::DataType) = Int64(readn(s, t))

"""
- Return `Int64` if in range.
- Otherwise fallback to `UInt64`.
"""
function read_u64(s::IO, t::DataType)
    v = UInt64(readn(s, t))
    if v > 2^63 - 1
        v
    else
        Int64(v)
    end
end

const UNPACK_OTHERS = Dict{UInt8, Function}(
    BIN_8  => s -> unpack_bin(s, readn(s, UInt8)),
    BIN_16 => s -> unpack_bin(s, readn(s, UInt16)),
    BIN_32 => s -> unpack_bin(s, readn(s, UInt32)),

    EXT_8  => s -> unpack_ext(s, readn(s, UInt8)),
    EXT_16 => s -> unpack_ext(s, readn(s, UInt16)),
    EXT_32 => s -> unpack_ext(s, readn(s, UInt32)),

    FLOAT_32 => s -> readn(s, Float32),
    FLOAT_64 => s -> readn(s, Float64),

    UINT_8  => s -> readi(s, UInt8),
    UINT_16 => s -> readi(s, UInt16),
    UINT_32 => s -> readi(s, UInt32),
    UINT_64 => s -> read_u64(s, UInt64),

    INT_8  => s -> readi(s, Int8),
    INT_16 => s -> readi(s, Int16),
    INT_32 => s -> readi(s, Int32),
    INT_64 => s -> readi(s, Int64),

    STR_8  => s -> unpack_str(s, readn(s, UInt8)),
    STR_16 => s -> unpack_str(s, readn(s, UInt16)),
    STR_32 => s -> unpack_str(s, readn(s, UInt32)),
    ARR_16 => s -> unpack_arr(s, readn(s, UInt16)),
    ARR_32 => s -> unpack_arr(s, readn(s, UInt32)),
    MAP_16 => s -> unpack_map(s, readn(s, UInt16)),
    MAP_32 => s -> unpack_map(s, readn(s, UInt32))
)

function unpack_others(s::IO, t::UInt8)
    f = try
        UNPACK_OTHERS[t]
    catch
        error("Wrong first byte $t")
    end
    f(s)
end

unpack(s) = unpack(IOBuffer(s))

function unpack(s::IO)
    b = read(s, UInt8)

    if b <= INT_FP[2]
        # positive fixint
        Int64(b)
    elseif b <= MAP_F[2]
        # fixmap
        unpack_map(s, b ⊻ MAP_F[1])
    elseif b <= ARR_F[2]
        # fixarray
        unpack_arr(s, b ⊻ ARR_F[1])
    elseif b <= STR_F[2]
        # fixstr
        unpack_str(s, b ⊻ STR_F[1])
    elseif EXT_F[1] <= b <= EXT_F[2]
        # fixext
        unpack_ext(s, 2^(b - EXT_F[1]))
    elseif b == NIL
        nothing
    elseif b == UNUSED
        error("unused")
    elseif b == FALSE
        false
    elseif b == TRUE
        true
    elseif b <= MAP_32
        unpack_others(s, b)
    else
        # negative fixint
        Int64(reinterpret(Int8, b))
    end
end

# Make sure key is evaluated before value.
unpack_map(s::IO, n::Integer) = Dict(unpack(s) => unpack(s) for _ in 1:n)
unpack_arr(s::IO, n::Integer) = [unpack(s) for _ in 1:n]
unpack_str(s::IO, n::Integer) = String(read(s, n))

function unpack_ext(s::IO, n::Integer)
    typecode = read(s, Int8)
    data = read(s, n)

    T = get(TYPECODE_TO_TYPE, typecode, nothing)
    if T === nothing
        Ext(typecode, data, impltype=true)
    else
        decode(T, data)::T
    end
end

unpack_bin(s::IO, n::Integer) = read(s, n)

function wh(io::IO, head, v)
    write(io, head)
    write(io, hton(v))
end

function pack(v)
    s = IOBuffer()
    pack(s, v)
    take!(s)
end

pack(s::IO, ::Nothing) = write(s, NIL)
pack(s::IO, v::Bool) = v ? write(s, TRUE) : write(s, FALSE)

function pack(s::IO, v::Integer)
    if v < 0
        if v >= -32
            write(s, Int8(v))
        elseif v >= -2^7
            wh(s, INT_8, Int8(v))
        elseif v >= -2^15
            wh(s, INT_16, Int16(v))
        elseif v >= -2^31
            wh(s, INT_32, Int32(v))
        elseif v >= -2^63
            wh(s, INT_64, Int64(v))
        else
            error("MsgPack signed int overflow")
        end
    else
        if v <= 127
            write(s, UInt8(v))
        elseif v <= 2^8 - 1
            wh(s, UINT_8, UInt8(v))
        elseif v <= 2^16 - 1
            wh(s, UINT_16, UInt16(v))
        elseif v <= 2^32 - 1
            wh(s, UINT_32, UInt32(v))
        elseif v <= UInt64(2)^64 - 1
            wh(s, UINT_64, UInt64(v))
        else
            error("MsgPack unsigned int overflow")
        end
    end
end

pack(s::IO, v::Float32) = wh(s, FLOAT_32, v)
pack(s::IO, v::Float64) = wh(s, FLOAT_64, v)

# str format
function pack(s::IO, v::AbstractString)
    n = sizeof(v)
    if n < 2^5
        write(s, STR_F[1] | UInt8(n))
    ## Note: with this section commented out, we do not have
    ##       the most compact format for a string.  However,
    ##       the string is still in spec, and some other
    ##       msgpack libaries (*ahem* Python) can't decode
    ##       strings created with this rule.
    #elseif n < 2^8
    #    wh(s, STR_8, UInt8(n))
    elseif n < 2^16
        wh(s, STR_16, UInt16(n))
    elseif n < 2^32
        wh(s, STR_32, UInt32(n))
    else
        error("MsgPack str overflow: ", n)
    end
    write(s, v)
end

# ext format
function pack(s::IO, v::Ext)
    n = sizeof(v.data)
    if n == 1
        write(s, 0xd4)
    elseif n == 2
        write(s, 0xd5)
    elseif n == 4
        write(s, 0xd6)
    elseif n == 8
        write(s, 0xd7)
    elseif n == 16
        write(s, 0xd8)
    elseif n < 2^8
        wh(s, EXT_8, UInt8(n))
    elseif n < 2^16
        wh(s, EXT_16, UInt16(n))
    elseif n < 2^32
        wh(s, EXT_32, UInt32(n))
    else
        error("MsgPack ext overflow: ", n)
    end
    write(s, v.typecode)
    write(s, v.data)
end

# bin format
function pack(s::IO, v::Vector{UInt8})
    n = length(v)
    if n < 2^8
        wh(s, BIN_8, UInt8(n))
    elseif n < 2^16
        wh(s, BIN_16, UInt16(n))
    elseif n < 2^32
        wh(s, BIN_32, UInt32(n))
    else
        error("MsgPack bin overflow: ", n)
    end
    write(s, v)
end

# Simple arrays
function pack(s::IO, v::Union{Vector, Tuple})
    n = length(v)
    if n < 2^4
        write(s, ARR_F[1] | UInt8(n))
    elseif n < 2^16
        wh(s, ARR_16, UInt16(n))
    elseif n < 2^32
        wh(s, ARR_32, UInt32(n))
    else
        error("MsgPack array overflow: ", n)
    end

    for x in v
        pack(s, x)
    end
end

# Maps
function pack(s::IO, v::Dict)
    n = length(v)
    if n < 2^4
        write(s, MAP_F[1] | UInt8(n))
    elseif n < 2^16
        wh(s, MAP_16, UInt16(n))
    elseif n < 2^32
        wh(s, MAP_32, UInt32(n))
    else
        error("MsgPack map overflow: ", n)
    end

    for (k, x) in v
        pack(s, k)
        pack(s, x)
    end
end

# Custom
function pack(s::IO, v::T) where T
    typecode = get_typecode(v)
    data = encode(v)::Vector{UInt8}
    pack(s, Ext(typecode, data))
end

end


text:
using Plots
using Random
using NearestNeighbors

global my_seed = 1001
Random.seed!(my_seed)

function inv_dist_itp_test(idxs, dists, p, values)
    avg_inv_dist = 0
    f = 0
    for i in 1:length(idxs)
        idx = idxs[i]
        inv_dist = 1/dists[i]^p
        avg_inv_dist += inv_dist
        f += values[idx]*inv_dist
    end
    f = f/avg_inv_dist
end

function inv_dist_itp(idxs::AbstractVector,
                      dists::AbstractVector,
                      p::AbstractFloat,
                      values::AbstractVector)
    avg_inv_dist = 0
    f = 0
    for i in 1:length(idxs)
        idx = idxs[i]
        if idx > 0
            inv_dist = 1/dists[i]^p
            avg_inv_dist += inv_dist
            f += values[idx]*inv_dist
        elseif idx == -5
            # lower boundary
            f+=0
        elseif idx == -6
            # upper boundary
            f+=0
        end
    end
    f = f/avg_inv_dist
end

function sample_function(x, y)
    return exp(-x^2 - y^2)
end

n_samples = 100

samples = rand(2, n_samples).*2 .-1
x_samples = samples[1,:]
y_samples = samples[2,:]

# Create a tree
tree = KDTree(samples)

k = 15

function_values = sample_function.(x_samples, y_samples)

raster_size = 500
x_range = LinRange(-1,1,raster_size)
y_range = LinRange(-1,1,raster_size)

gr()

function plot_the_itp(p)

    raster_values = Matrix{Float64}(undef, (raster_size, raster_size))

    for i in 1:raster_size
        for j in 1:raster_size
            raster_site = [x_range[i], y_range[j]]
            nearest_neighbours = knn(tree, raster_site, k)
            raster_values[i, j] = inv_dist_itp_test(nearest_neighbours[1],
                                                    nearest_neighbours[2],
                                                    p,
                                                    function_values)
        end
    end
    heatmap(x_range,
            y_range,
            raster_values,
            title="p=$p, $k neighbours",
            dpi=350,
            aspect_ratio=:equal,
            xlim=[-1,1])
    savefig("../img/InverseDistance/inv_dist_$(Int(p*100))_$(k)")
end

p = [3]
for i in p
    plot_the_itp(i)
end

true_values = Matrix{Float64}(undef, (raster_size, raster_size))
for i in 1:raster_size
    for j in 1:raster_size
        true_values[i, j] = sample_function(x_range[i], y_range[j])
    end
end

heatmap(x_range,
        y_range,
        true_values,
        title="True Function",
        dpi=350,
        aspect_ratio=:equal,
        xlim=[-1,1])
savefig("../img/InverseDistance/inv_dist_comparison")

scatter(x_samples,
        y_samples,
        title="samples",
        dpi=350,
        aspect_ratio=:equal,
        xlim=[-1,1],
        xlabel=raw"$x$",
        ylabel=raw"$y$")
savefig("../img/InverseDistance/sample_points")

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
#################################
################################# MSE tables
#################################
#
all_best_mse = OrderedCollections.OrderedDict()
#
all_best_mse["ARM_SGP"] = OrderedCollections.OrderedDict()
all_best_mse["ARM_SGP"]["qt_mean"] = 0.2099693151607237
all_best_mse["ARM_SGP"]["updraft_area"] = 335.45875367066867
all_best_mse["ARM_SGP"]["updraft_w"] = 140.92677349078969
all_best_mse["ARM_SGP"]["updraft_qt"] = 28.70836894541848
all_best_mse["ARM_SGP"]["updraft_thetal"] = 172.01631026298847
all_best_mse["ARM_SGP"]["u_mean"] = 1.3235797273549681e-5
all_best_mse["ARM_SGP"]["tke_mean"] = 1313.6493467020196
all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00010624950402975059
all_best_mse["ARM_SGP"]["ql_mean"] = 180.68590465680364
all_best_mse["ARM_SGP"]["qi_mean"] = "NA"
all_best_mse["ARM_SGP"]["thetal_mean"] = 9.770061871604239e-5
all_best_mse["ARM_SGP"]["Hvar_mean"] = 9311.054396714339
all_best_mse["ARM_SGP"]["QTvar_mean"] = 5802.412820315867
#
all_best_mse["Bomex"] = OrderedCollections.OrderedDict()
all_best_mse["Bomex"]["qt_mean"] = 0.08776796077692245
all_best_mse["Bomex"]["updraft_area"] = 130.01639628691083
all_best_mse["Bomex"]["updraft_w"] = 16.057586429511503
all_best_mse["Bomex"]["updraft_qt"] = 6.890128846374632
all_best_mse["Bomex"]["updraft_thetal"] = 69.78417329537506
all_best_mse["Bomex"]["v_mean"] = 67.12281585199479
all_best_mse["Bomex"]["u_mean"] = 0.3201493689027362
all_best_mse["Bomex"]["tke_mean"] = 73.89824775455659
all_best_mse["Bomex"]["temperature_mean"] = 3.252413195141947e-5
all_best_mse["Bomex"]["ql_mean"] = 6.313648778401822
all_best_mse["Bomex"]["qi_mean"] = "NA"
all_best_mse["Bomex"]["thetal_mean"] = 3.2858173111125033e-5
all_best_mse["Bomex"]["Hvar_mean"] = 6752.843931363311
all_best_mse["Bomex"]["QTvar_mean"] = 2345.9671498346493
#
all_best_mse["DryBubble"] = OrderedCollections.OrderedDict()
all_best_mse["DryBubble"]["updraft_area"] = 0.0
all_best_mse["DryBubble"]["updraft_w"] = 0.0
all_best_mse["DryBubble"]["updraft_thetal"] = 0.0
all_best_mse["DryBubble"]["u_mean"] = 0.0
all_best_mse["DryBubble"]["tke_mean"] = 0.0
all_best_mse["DryBubble"]["temperature_mean"] = 0.0
all_best_mse["DryBubble"]["thetal_mean"] = 0.0
all_best_mse["DryBubble"]["Hvar_mean"] = 0.0
#
all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict()
all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.023198692272030544
all_best_mse["DYCOMS_RF01"]["ql_mean"] = 10.347286087600702
all_best_mse["DYCOMS_RF01"]["updraft_area"] = 30.41218602111556
all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.1304092648269375
all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.1537875302070426
all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.18817055108242
all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.002982833065277946
all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.07894988376751147
all_best_mse["DYCOMS_RF01"]["tke_mean"] = 21.327285272822852
all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 7.606763938668142e-5
all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 7.644917953938347e-5
all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1288.1645081619606
all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 518.9286165737209
#
all_best_mse["GABLS"] = OrderedCollections.OrderedDict()
all_best_mse["GABLS"]["updraft_thetal"] = 0.0
all_best_mse["GABLS"]["v_mean"] = 0.0
all_best_mse["GABLS"]["u_mean"] = 0.0
all_best_mse["GABLS"]["tke_mean"] = 0.0
all_best_mse["GABLS"]["temperature_mean"] = 0.0
all_best_mse["GABLS"]["thetal_mean"] = 0.0
all_best_mse["GABLS"]["Hvar_mean"] = 0.0
all_best_mse["GABLS"]["QTvar_mean"] = 0.0
all_best_mse["GABLS"]["qt_mean"] = 0.0
#
all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict()
all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.0
all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.0
all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 0.0
all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 0.0
all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.0
all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 0.0
all_best_mse["life_cycle_Tan2018"]["v_mean"] = 0.0
all_best_mse["life_cycle_Tan2018"]["u_mean"] = 0.0
all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 0.0
all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 0.0
all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 0.0
all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 0.0
all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 0.0
#
all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict()
all_best_mse["Nieuwstadt"]["updraft_area"] = 98.61903956765104
all_best_mse["Nieuwstadt"]["updraft_w"] = 12.620181519729721
all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.60674820900034
all_best_mse["Nieuwstadt"]["u_mean"] = 13.627556919779142
all_best_mse["Nieuwstadt"]["tke_mean"] = 312.92521835084375
all_best_mse["Nieuwstadt"]["temperature_mean"] = 9.606319451874596e-6
all_best_mse["Nieuwstadt"]["thetal_mean"] = 9.918979982485869e-6
all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1280.538121964536
#
all_best_mse["Rico"] = OrderedCollections.OrderedDict()
all_best_mse["Rico"]["qt_mean"] = 1.461216267731136
all_best_mse["Rico"]["updraft_area"] = 476.1668575277046
all_best_mse["Rico"]["updraft_w"] = 103.62523342741915
all_best_mse["Rico"]["updraft_qt"] = 12.084264582170741
all_best_mse["Rico"]["updraft_thetal"] = 133.77973866180497
all_best_mse["Rico"]["v_mean"] = 0.661591213386313
all_best_mse["Rico"]["u_mean"] = 0.5942324640298271
all_best_mse["Rico"]["tke_mean"] = 90.06304177194619
all_best_mse["Rico"]["temperature_mean"] = 0.0006561578714960826
all_best_mse["Rico"]["ql_mean"] = 69.79123138670636
all_best_mse["Rico"]["qi_mean"] = "NA"
all_best_mse["Rico"]["qr_mean"] = 765.9753444576647
all_best_mse["Rico"]["thetal_mean"] = 0.0006459491204637934
all_best_mse["Rico"]["Hvar_mean"] = 20043.952395704495
all_best_mse["Rico"]["QTvar_mean"] = 3977.214921455993
#
all_best_mse["Soares"] = OrderedCollections.OrderedDict()
all_best_mse["Soares"]["qt_mean"] = 0.1257904686397003
all_best_mse["Soares"]["updraft_area"] = 98.21922321246183
all_best_mse["Soares"]["updraft_w"] = 11.423147273156957
all_best_mse["Soares"]["updraft_qt"] = 23.615121270329094
all_best_mse["Soares"]["updraft_thetal"] = 65.7219396386948
all_best_mse["Soares"]["u_mean"] = 94.2923650437216
all_best_mse["Soares"]["tke_mean"] = 241.22739244519957
all_best_mse["Soares"]["temperature_mean"] = 1.0913731129957992e-5
all_best_mse["Soares"]["thetal_mean"] = 1.0343924750695743e-5
all_best_mse["Soares"]["Hvar_mean"] = 1118.8459567581533
#
all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict()
all_best_mse["TRMM_LBA"]["qt_mean"] = 4.370141509071347
all_best_mse["TRMM_LBA"]["updraft_area"] = 7588.792583634186
all_best_mse["TRMM_LBA"]["updraft_w"] = 30830.63975659635
all_best_mse["TRMM_LBA"]["updraft_qt"] = 263.1587569415405
all_best_mse["TRMM_LBA"]["updraft_thetal"] = 1961.8451637527064
all_best_mse["TRMM_LBA"]["v_mean"] = 285.9190262131904
all_best_mse["TRMM_LBA"]["u_mean"] = 113.61185296026021
all_best_mse["TRMM_LBA"]["tke_mean"] = 25299.781769411315
all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.001282306577221344
all_best_mse["TRMM_LBA"]["ql_mean"] = 10591.308689650668
all_best_mse["TRMM_LBA"]["qi_mean"] = "NA"
all_best_mse["TRMM_LBA"]["qr_mean"] = "NA"
all_best_mse["TRMM_LBA"]["qs_mean"] = "NA"
all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.021142841574681895
all_best_mse["TRMM_LBA"]["Hvar_mean"] = 4196.816226656944
all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2772.6720734920077
#
all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict()
all_best_mse["LES_driven_SCM"]["qt_mean"] = 3.5617162347368296
all_best_mse["LES_driven_SCM"]["v_mean"] = 1.3028760898968443
all_best_mse["LES_driven_SCM"]["u_mean"] = 0.44955547910863786
all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.001277413538359007
all_best_mse["LES_driven_SCM"]["ql_mean"] = 55700.59768552013
all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0015156834284607626
#
#################################
#################################
#################################


text:
module VIE
    using ..BEAST
    Mod = BEAST

    """
        singlelayer(;gamma, alpha, beta, tau)
        singlelayer(;wavenumber, alpha, beta, tau)

    Bilinear form given by:

    ```math
        α ∬_{Ω×Ω} j(x)⋅τ(y)⋅k(y) G_{γ}(x,y) + β ∬_{Ω×Ω} div j(x) τ(y)⋅k(y)⋅grad G_{γ}(x,y)
    ```

    with ``G_{γ} = e^{-γ|x-y|} / 4π|x-y|``
    and  ``τ(y)`` contrast dyadic
    """

    function singlelayer(;
            gamma=nothing,
            wavenumber=nothing,
            alpha=nothing,
            beta=nothing,
            tau=nothing)

        if (gamma == nothing) && (wavenumber == nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if (gamma != nothing) && (wavenumber != nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if gamma == nothing
            if iszero(real(wavenumber))
                gamma = -imag(wavenumber)
            else
                gamma = im*wavenumber
            end
        end

        @assert gamma != nothing

        alpha == nothing && (alpha = wavenumber*wavenumber)
        beta  == nothing && (beta  = 1.0)
        tau   == nothing && (tau   = x->1.0)

        Mod.VIESingleLayer(gamma, alpha, beta, tau)
    end

    function singlelayer2(;
            gamma=nothing,
            wavenumber=nothing,
            alpha=nothing,
            beta=nothing,
            tau=nothing)

        if (gamma == nothing) && (wavenumber == nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if (gamma != nothing) && (wavenumber != nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if gamma == nothing
            if iszero(real(wavenumber))
                gamma = -imag(wavenumber)
            else
                gamma = im*wavenumber
            end
        end

        @assert gamma != nothing

        alpha == nothing && (alpha = 1.0)
        beta  == nothing && (beta  = 1.0)
        tau   == nothing && (tau   = x->1.0)

        Mod.VIESingleLayer2(gamma, alpha, beta, tau)
    end

    """
    boundary(;gamma, alpha, tau)
    boundary(;wavenumber, alpha, tau)

    Bilinear form given by:

    ```math
        α ∬_{Ω×Ω} T_n j(x) grad G_{γ}(x,y)⋅τ(y)⋅k(y)
    ```

    with ``G_{γ} = e^{-γ|x-y|} / 4π|x-y|``
    and  ``τ(y)`` contrast dyadic
    """

    function boundary(;
            gamma=nothing,
            wavenumber=nothing,
            alpha=nothing,
            tau=nothing)

        if (gamma == nothing) && (wavenumber == nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if (gamma != nothing) && (wavenumber != nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if gamma == nothing
            if iszero(real(wavenumber))
                gamma = -imag(wavenumber)
            else
                gamma = im*wavenumber
            end
        end

        @assert gamma != nothing

        alpha == nothing && (alpha = 1.0)
        tau   == nothing && (tau   = x->1.0)

        Mod.VIEBoundary(gamma, alpha, tau)
    end

    function boundary2(;
            gamma=nothing,
            wavenumber=nothing,
            alpha=nothing,
            tau=nothing)

        if (gamma == nothing) && (wavenumber == nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if (gamma != nothing) && (wavenumber != nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if gamma == nothing
            if iszero(real(wavenumber))
                gamma = -imag(wavenumber)
            else
                gamma = im*wavenumber
            end
        end

        @assert gamma != nothing

        alpha == nothing && (alpha = 1.0)
        tau   == nothing && (tau   = x->1.0)

        Mod.VIEBoundary2(gamma, alpha, tau)
    end

    """
    doublelayer(;gamma, alpha, beta, tau)
    doublelayer(;wavenumber, alpha, beta, tau)

    Bilinear form given by:

    ```math
        α ∬_{Ω×Ω} j(x) ⋅ grad G_{γ}(x,y) × τ(y)⋅k(y)
    ```

    with ``G_{γ} = e^{-γ|x-y|} / 4π|x-y|``
    and  ``τ(y)`` contrast dyadic
    """

    function doublelayer(;
            gamma=nothing,
            wavenumber=nothing,
            alpha=nothing,
            beta=nothing,
            tau=nothing)

        if (gamma == nothing) && (wavenumber == nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if (gamma != nothing) && (wavenumber != nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if gamma == nothing
            if iszero(real(wavenumber))
                gamma = -imag(wavenumber)
            else
                gamma = im*wavenumber
            end
        end

        @assert gamma != nothing

        alpha == nothing && (alpha = 1.0)
        tau   == nothing && (tau   = x->1.0)

        Mod.VIEDoubleLayer(gamma, alpha, tau)
    end

    planewave(;
        direction    = error("missing arguement `direction`"),
        polarization = error("missing arguement `polarization`"),
        wavenumber   = error("missing arguement `wavenumber`"),
        amplitude    = one(real(typeof(wavenumber)))) =
    Mod.PlaneWaveVIE(direction, polarization, wavenumber, amplitude)

    farfield(;
        wavenumber = error("missing argument: `wavenumber`"),
        tau=error("missing argument: `tau`")) =
            Mod.VIEFarField3D(wavenumber=wavenumber, tau= tau)

end


text:
# generateur d'Instances pour la RODD projet 3

N = 8
Nm = 4
Nf = 4
C = 1
G = 5
A = 2
T = 50
init = 0.001

# Ouvrir le fichier "output.txt" dans lequel on pourra écrire
fout = open("output.txt", "w")

println(fout, "title")
println(fout, "title")
println(fout, "title")
println(fout, "title")

println(fout, "N = " * string(N))
println(fout, "Nm = " * string(Nm))
println(fout, "Nf = " * string(Nf))
println(fout, "C = " * string(C))
println(fout, "G = " * string(G))
println(fout, "A = " * string(A))
println(fout, "T = " * string(T))
println(fout, "init = " * string(init))

for i in 1:N
  for j in 1:G
    myLine = string(i) * " " * string(1) * " " * string(j) * " " * string(1) * " " * string(rand(1:A))
    println(fout, myLine)
    myLine = string(i) * " " * string(1) * " " * string(j) * " " * string(2) * " " * string(rand(1:A))
    println(fout, myLine)
  end
end

close(fout)

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
using GeoStatsBase

"""Interpolate missing values in GeoArray."""
function interpolate!(ga::GeoArray, solver::T, band=1, symbol=:z) where T <: AbstractSolver
    # Irregular grid
    # TODO Use unstructured GeoStats method
    if is_rotated(ga)
        error("Can't interpolate warped grid yet.")

    # Regular grid
    else
        rg = RegularGrid(size(ga.A)[1:2], Tuple(ga.f.translation), (abs(ga.f.linear[1]), abs(ga.f.linear[4])))
        problemdata = georef(NamedTuple{(symbol,)}((ga.A[:,:,band],)), rg)
    end

    problem = EstimationProblem(problemdata, rg, symbol, mapper=CopyMapper())
    solution = solve(problem, solver)

    ga.A[:, :, band] .= reshape(solution[symbol][:mean], size(ga.A)[1:2])
    ga
end

export interpolate!


text:
#randomgrowth2.jl
using PyPlot
using ODE

function randomgrowth2()
    num_trials = 2000
    g = 1
    q = 0.7
    Ns = 80
    # Ns = [50, 100, 200]
    clf()
    # [-1.771086807411  08131947928329]
    # Gap = c d^N
    for jj = 1:length(Ns)
        N = Ns[jj]
        M = round(N*g)
        B = zeros(num_trials)
        @time begin
            for ii = 1:num_trials
                # G = G[M, N, q]
                # B[ii] = G[M,N]
                B[ii] = G(M, N, q)
            end
        end
        C = (B - N*om(g,q)) / (sigma_growth(g,q)*N^(1/3))
        d = 0.2
        subplot(1,length(Ns),jj)
        plt.hist(C - exp(-N*0.0025 - 3), normed = true)
        xlim(-6, 4)

        ## Theory
        t0 = 4
        tn = -6
        dx = 0.005
        deq = (t, y) -> [y[2]; t*y[1]+2*y[1]^3; y[4]; y[1]^2]
        y0 = [airy(t0); airy(1,t0); 0; airy(t0)^2]      # boundary conditions
        t, y = ode23(deq, y0, t0:-dx:tn)                # solve
        F2 = Float64[exp(-y[i][3]) for i = 1:length(y)]                               # the distribution
        f2 = gradient(F2, t)                            # the density

        # add(p, Curve(t, f2, "color", "red", "linewidth", 3))
        # Winston.display(p)
        subplot(1,length(Ns),jj)
        plot(t, f2, "r", linewidth = 3)
        ylim(0, 0.6)
        println(mean(C))
    end
end

function G(N, M, q)
    # computes matrix G[N,M] for a given q
    GG = floor(log(rand(N,M))/log(q))
    #float(rand(N,M) < q)

    # Compute the edges: the boundary
    for ii = 2:N
        GG[ii, 1] = GG[ii, 1] + GG[ii-1, 1]
        GG[1, ii] = GG[1, ii] + GG[1, ii-1]
    end
    # Compute the inside
    for ii = 2:N
        for jj = 2:M
            GG[ii, jj] = GG[ii, jj] + max(GG[ii, jj-1], GG[ii-1, jj])
        end
    end
    ## Plot
    # imagesc((0,N),(M,0),GG)
    # sleep(.01)

    return  GG[N, M]
end

# helper functions sigma.m and om, which describe mean and
# standard deviation of G[N, M]

# Computes G[qN,N]/N for N->inf
om(g,q) = ((1 + sqrt(g * q)) ^ 2) / (1 - q) - 1

sigma_growth(r, q) =
    q^(1/6) * r^(1/6) / (1 - q) *
    (sqrt(r) + sqrt(q))^(2/3) *
    (1 + sqrt(q * r))^(2/3)

randomgrowth2()


text:
@testset "4.6.0 (a csc)^m (b trg)^n" begin
    (a, b, c, d, e, f, m, n, p, x, ) = @variables a b c d e f m n p x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(b*csc(c+d*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(b*csc(c+d*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*csc(c+d*x)^n=#

    @test_int [csc(a + b*x)^1, x, 1, -(atanh(cos(a + b*x))/b)]
    @test_int [csc(a + b*x)^2, x, 2, -(cot(a + b*x)/b)]
    @test_int [csc(a + b*x)^3, x, 2, -(atanh(cos(a + b*x))/(2*b)) - (cot(a + b*x)*csc(a + b*x))/(2*b)]
    @test_int [csc(a + b*x)^4, x, 2, -(cot(a + b*x)/b) - cot(a + b*x)^3/(3*b)]
    @test_int [csc(a + b*x)^5, x, 3, -((3*atanh(cos(a + b*x)))/(8*b)) - (3*cot(a + b*x)*csc(a + b*x))/(8*b) - (cot(a + b*x)*csc(a + b*x)^3)/(4*b)]
    @test_int [csc(a + b*x)^6, x, 2, -(cot(a + b*x)/b) - (2*cot(a + b*x)^3)/(3*b) - cot(a + b*x)^5/(5*b)]
    @test_int [csc(a + b*x)^7, x, 4, -((5*atanh(cos(a + b*x)))/(16*b)) - (5*cot(a + b*x)*csc(a + b*x))/(16*b) - (5*cot(a + b*x)*csc(a + b*x)^3)/(24*b) - (cot(a + b*x)*csc(a + b*x)^5)/(6*b)]
    @test_int [csc(a + b*x)^8, x, 2, -(cot(a + b*x)/b) - cot(a + b*x)^3/b - (3*cot(a + b*x)^5)/(5*b) - cot(a + b*x)^7/(7*b)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*csc(c+d*x))^(n/2)=#

    @test_int [csc(a + b*x)^(7/2), x, 4, -((6*cos(a + b*x)*sqrt(csc(a + b*x)))/(5*b)) - (2*cos(a + b*x)*csc(a + b*x)^(5/2))/(5*b) - (6*sqrt(csc(a + b*x))*Elliptic.E((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/(5*b)]
    @test_int [csc(a + b*x)^(5/2), x, 3, -((2*cos(a + b*x)*csc(a + b*x)^(3/2))/(3*b)) + (2*sqrt(csc(a + b*x))*Elliptic.F((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/(3*b)]
    @test_int [csc(a + b*x)^(3/2), x, 3, -((2*cos(a + b*x)*sqrt(csc(a + b*x)))/b) - (2*sqrt(csc(a + b*x))*Elliptic.E((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/b]
    @test_int [csc(a + b*x)^(1/2), x, 2, (2*sqrt(csc(a + b*x))*Elliptic.F((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/b]
    @test_int [1/csc(a + b*x)^(1/2), x, 2, (2*sqrt(csc(a + b*x))*Elliptic.E((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/b]
    @test_int [1/csc(a + b*x)^(3/2), x, 3, -((2*cos(a + b*x))/(3*b*sqrt(csc(a + b*x)))) + (2*sqrt(csc(a + b*x))*Elliptic.F((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/(3*b)]
    @test_int [1/csc(a + b*x)^(5/2), x, 3, -((2*cos(a + b*x))/(5*b*csc(a + b*x)^(3/2))) + (6*sqrt(csc(a + b*x))*Elliptic.E((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/(5*b)]
    @test_int [1/csc(a + b*x)^(7/2), x, 4, -((2*cos(a + b*x))/(7*b*csc(a + b*x)^(5/2))) - (10*cos(a + b*x))/(21*b*sqrt(csc(a + b*x))) + (10*sqrt(csc(a + b*x))*Elliptic.F((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/(21*b)]

    @test_int [(c*csc(a + b*x))^(7/2), x, 4, -((6*c^3*cos(a + b*x)*sqrt(c*csc(a + b*x)))/(5*b)) - (2*c*cos(a + b*x)*(c*csc(a + b*x))^(5/2))/(5*b) - (6*c^4*Elliptic.E((1/2)*(a - pi/2 + b*x), 2))/(5*b*sqrt(c*csc(a + b*x))*sqrt(sin(a + b*x)))]
    @test_int [(c*csc(a + b*x))^(5/2), x, 3, -((2*c*cos(a + b*x)*(c*csc(a + b*x))^(3/2))/(3*b)) + (2*c^2*sqrt(c*csc(a + b*x))*Elliptic.F((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/(3*b)]
    @test_int [(c*csc(a + b*x))^(3/2), x, 3, -((2*c*cos(a + b*x)*sqrt(c*csc(a + b*x)))/b) - (2*c^2*Elliptic.E((1/2)*(a - pi/2 + b*x), 2))/(b*sqrt(c*csc(a + b*x))*sqrt(sin(a + b*x)))]
    @test_int [(c*csc(a + b*x))^(1/2), x, 2, (2*sqrt(c*csc(a + b*x))*Elliptic.F((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/b]
    @test_int [1/(c*csc(a + b*x))^(1/2), x, 2, (2*Elliptic.E((1/2)*(a - pi/2 + b*x), 2))/(b*sqrt(c*csc(a + b*x))*sqrt(sin(a + b*x)))]
    @test_int [1/(c*csc(a + b*x))^(3/2), x, 3, -((2*cos(a + b*x))/(3*b*c*sqrt(c*csc(a + b*x)))) + (2*sqrt(c*csc(a + b*x))*Elliptic.F((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/(3*b*c^2)]
    @test_int [1/(c*csc(a + b*x))^(5/2), x, 3, -((2*cos(a + b*x))/(5*b*c*(c*csc(a + b*x))^(3/2))) + (6*Elliptic.E((1/2)*(a - pi/2 + b*x), 2))/(5*b*c^2*sqrt(c*csc(a + b*x))*sqrt(sin(a + b*x)))]
    @test_int [1/(c*csc(a + b*x))^(7/2), x, 4, -((2*cos(a + b*x))/(7*b*c*(c*csc(a + b*x))^(5/2))) - (10*cos(a + b*x))/(21*b*c^3*sqrt(c*csc(a + b*x))) + (10*sqrt(c*csc(a + b*x))*Elliptic.F((1/2)*(a - pi/2 + b*x), 2)*sqrt(sin(a + b*x)))/(21*b*c^4)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*csc(c+d*x))^(n/3)=#

    @test_int [csc(a + b*x)^(4/3), x, 2, -((3*cos(a + b*x)*csc(a + b*x)^(1/3)*HypergeometricFunctions._₂F₁(-(1/6), 1/2, 5/6, sin(a + b*x)^2))/(b*sqrt(cos(a + b*x)^2)))]
    @test_int [csc(a + b*x)^(2/3), x, 2, (3*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/6, 1/2, 7/6, sin(a + b*x)^2))/(b*sqrt(cos(a + b*x)^2)*csc(a + b*x)^(1/3))]
    @test_int [csc(a + b*x)^(1/3), x, 2, (3*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/3, 1/2, 4/3, sin(a + b*x)^2))/(2*b*sqrt(cos(a + b*x)^2)*csc(a + b*x)^(2/3))]
    @test_int [1/csc(a + b*x)^(1/3), x, 2, (3*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/2, 2/3, 5/3, sin(a + b*x)^2))/(4*b*sqrt(cos(a + b*x)^2)*csc(a + b*x)^(4/3))]
    @test_int [1/csc(a + b*x)^(2/3), x, 2, (3*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/2, 5/6, 11/6, sin(a + b*x)^2))/(5*b*sqrt(cos(a + b*x)^2)*csc(a + b*x)^(5/3))]
    @test_int [1/csc(a + b*x)^(4/3), x, 2, (3*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/2, 7/6, 13/6, sin(a + b*x)^2))/(7*b*sqrt(cos(a + b*x)^2)*csc(a + b*x)^(7/3))]

    @test_int [(c*csc(a + b*x))^(4/3), x, 2, -((3*c*cos(a + b*x)*(c*csc(a + b*x))^(1/3)*HypergeometricFunctions._₂F₁(-(1/6), 1/2, 5/6, sin(a + b*x)^2))/(b*sqrt(cos(a + b*x)^2)))]
    @test_int [(c*csc(a + b*x))^(2/3), x, 2, (3*c*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/6, 1/2, 7/6, sin(a + b*x)^2))/(b*sqrt(cos(a + b*x)^2)*(c*csc(a + b*x))^(1/3))]
    @test_int [(c*csc(a + b*x))^(1/3), x, 2, (3*c*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/3, 1/2, 4/3, sin(a + b*x)^2))/(2*b*sqrt(cos(a + b*x)^2)*(c*csc(a + b*x))^(2/3))]
    @test_int [1/(c*csc(a + b*x))^(1/3), x, 2, (3*c*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/2, 2/3, 5/3, sin(a + b*x)^2))/(4*b*sqrt(cos(a + b*x)^2)*(c*csc(a + b*x))^(4/3))]
    @test_int [1/(c*csc(a + b*x))^(2/3), x, 2, (3*c*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/2, 5/6, 11/6, sin(a + b*x)^2))/(5*b*sqrt(cos(a + b*x)^2)*(c*csc(a + b*x))^(5/3))]
    @test_int [1/(c*csc(a + b*x))^(4/3), x, 2, (3*c*cos(a + b*x)*HypergeometricFunctions._₂F₁(1/2, 7/6, 13/6, sin(a + b*x)^2))/(7*b*sqrt(cos(a + b*x)^2)*(c*csc(a + b*x))^(7/3))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*csc(c+d*x))^n*with*n*symbolic=#

    @test_int [csc(a + b*x)^n, x, 2, (cos(a + b*x)*csc(a + b*x)^(-1 + n)*HypergeometricFunctions._₂F₁(1/2, (1 - n)/2, (3 - n)/2, sin(a + b*x)^2))/(b*(1 - n)*sqrt(cos(a + b*x)^2))]

    @test_int [(c*csc(a + b*x))^n, x, 2, (c*cos(a + b*x)*(c*csc(a + b*x))^(-1 + n)*HypergeometricFunctions._₂F₁(1/2, (1 - n)/2, (3 - n)/2, sin(a + b*x)^2))/(b*(1 - n)*sqrt(cos(a + b*x)^2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(b*csc(c+d*x)^p)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*csc(c+d*x)^p)^(n/2)*with*p*positive*integer=#

    @test_int [(csc(x)^2)^(7/2), x, 5, (-(5/16))*asinh(cot(x)) - (5/16)*cot(x)*sqrt(csc(x)^2) - (5/24)*cot(x)*(csc(x)^2)^(3/2) - (1/6)*cot(x)*(csc(x)^2)^(5/2)]
    @test_int [(csc(x)^2)^(5/2), x, 4, (-(3/8))*asinh(cot(x)) - (3/8)*cot(x)*sqrt(csc(x)^2) - (1/4)*cot(x)*(csc(x)^2)^(3/2)]
    @test_int [(csc(x)^2)^(3/2), x, 3, (-(1/2))*asinh(cot(x)) - (1/2)*cot(x)*sqrt(csc(x)^2)]
    @test_int [(csc(x)^2)^(1/2), x, 2, -asinh(cot(x))]
    @test_int [1/(csc(x)^2)^(1/2), x, 2, -(cot(x)/sqrt(csc(x)^2))]
    @test_int [1/(csc(x)^2)^(3/2), x, 3, -(cot(x)/(3*(csc(x)^2)^(3/2))) - (2*cot(x))/(3*sqrt(csc(x)^2))]
    @test_int [1/(csc(x)^2)^(5/2), x, 4, -(cot(x)/(5*(csc(x)^2)^(5/2))) - (4*cot(x))/(15*(csc(x)^2)^(3/2)) - (8*cot(x))/(15*sqrt(csc(x)^2))]
    @test_int [1/(csc(x)^2)^(7/2), x, 5, -(cot(x)/(7*(csc(x)^2)^(7/2))) - (6*cot(x))/(35*(csc(x)^2)^(5/2)) - (8*cot(x))/(35*(csc(x)^2)^(3/2)) - (16*cot(x))/(35*sqrt(csc(x)^2))]

    @test_int [(a*csc(x)^2)^(7/2), x, 6, (-(5/16))*a^(7/2)*atanh((sqrt(a)*cot(x))/sqrt(a*csc(x)^2)) - (5/16)*a^3*cot(x)*sqrt(a*csc(x)^2) - (5/24)*a^2*cot(x)*(a*csc(x)^2)^(3/2) - (1/6)*a*cot(x)*(a*csc(x)^2)^(5/2)]
    @test_int [(a*csc(x)^2)^(5/2), x, 5, (-(3/8))*a^(5/2)*atanh((sqrt(a)*cot(x))/sqrt(a*csc(x)^2)) - (3/8)*a^2*cot(x)*sqrt(a*csc(x)^2) - (1/4)*a*cot(x)*(a*csc(x)^2)^(3/2)]
    @test_int [(a*csc(x)^2)^(3/2), x, 4, (-(1/2))*a^(3/2)*atanh((sqrt(a)*cot(x))/sqrt(a*csc(x)^2)) - (1/2)*a*cot(x)*sqrt(a*csc(x)^2)]
    @test_int [(a*csc(x)^2)^(1/2), x, 3, (-sqrt(a))*atanh((sqrt(a)*cot(x))/sqrt(a*csc(x)^2))]
    @test_int [1/(a*csc(x)^2)^(1/2), x, 2, -(cot(x)/sqrt(a*csc(x)^2))]
    @test_int [1/(a*csc(x)^2)^(3/2), x, 3, -(cot(x)/(3*(a*csc(x)^2)^(3/2))) - (2*cot(x))/(3*a*sqrt(a*csc(x)^2))]
    @test_int [1/(a*csc(x)^2)^(5/2), x, 4, -(cot(x)/(5*(a*csc(x)^2)^(5/2))) - (4*cot(x))/(15*a*(a*csc(x)^2)^(3/2)) - (8*cot(x))/(15*a^2*sqrt(a*csc(x)^2))]
    @test_int [1/(a*csc(x)^2)^(7/2), x, 5, -(cot(x)/(7*(a*csc(x)^2)^(7/2))) - (6*cot(x))/(35*a*(a*csc(x)^2)^(5/2)) - (8*cot(x))/(35*a^2*(a*csc(x)^2)^(3/2)) - (16*cot(x))/(35*a^3*sqrt(a*csc(x)^2))]

    @test_int [(a*csc(x)^3)^(5/2), x, 7, (-(154/585))*a^2*cot(x)*sqrt(a*csc(x)^3) - (22/117)*a^2*cot(x)*csc(x)^2*sqrt(a*csc(x)^3) - (2/13)*a^2*cot(x)*csc(x)^4*sqrt(a*csc(x)^3) - (154/195)*a^2*cos(x)*sqrt(a*csc(x)^3)*sin(x) + (154/195)*a^2*sqrt(a*csc(x)^3)*Elliptic.E(pi/4 - x/2, 2)*sin(x)^(3/2)]
    @test_int [(a*csc(x)^3)^(3/2), x, 5, (-(10/21))*a*cos(x)*sqrt(a*csc(x)^3) - (2/7)*a*cot(x)*csc(x)*sqrt(a*csc(x)^3) - (10/21)*a*sqrt(a*csc(x)^3)*Elliptic.F(pi/4 - x/2, 2)*sin(x)^(3/2)]
    @test_int [(a*csc(x)^3)^(1/2), x, 4, -2*cos(x)*sqrt(a*csc(x)^3)*sin(x) + 2*sqrt(a*csc(x)^3)*Elliptic.E(pi/4 - x/2, 2)*sin(x)^(3/2)]
    @test_int [1/(a*csc(x)^3)^(1/2), x, 4, -((2*cot(x))/(3*sqrt(a*csc(x)^3))) - (2*Elliptic.F(pi/4 - x/2, 2))/(3*sqrt(a*csc(x)^3)*sin(x)^(3/2))]
    @test_int [1/(a*csc(x)^3)^(3/2), x, 5, -((14*cos(x))/(45*a*sqrt(a*csc(x)^3))) - (14*Elliptic.E(pi/4 - x/2, 2))/(15*a*sqrt(a*csc(x)^3)*sin(x)^(3/2)) - (2*cos(x)*sin(x)^2)/(9*a*sqrt(a*csc(x)^3))]
    @test_int [1/(a*csc(x)^3)^(5/2), x, 7, -((26*cot(x))/(77*a^2*sqrt(a*csc(x)^3))) - (26*Elliptic.F(pi/4 - x/2, 2))/(77*a^2*sqrt(a*csc(x)^3)*sin(x)^(3/2)) - (78*cos(x)*sin(x))/(385*a^2*sqrt(a*csc(x)^3)) - (26*cos(x)*sin(x)^3)/(165*a^2*sqrt(a*csc(x)^3)) - (2*cos(x)*sin(x)^5)/(15*a^2*sqrt(a*csc(x)^3))]

    @test_int [(a*csc(x)^4)^(7/2), x, 3, -2*a^3*cos(x)^2*cot(x)*sqrt(a*csc(x)^4) - 3*a^3*cos(x)^2*cot(x)^3*sqrt(a*csc(x)^4) - (20/7)*a^3*cos(x)^2*cot(x)^5*sqrt(a*csc(x)^4) - (5/3)*a^3*cos(x)^2*cot(x)^7*sqrt(a*csc(x)^4) - (6/11)*a^3*cos(x)^2*cot(x)^9*sqrt(a*csc(x)^4) - (1/13)*a^3*cos(x)^2*cot(x)^11*sqrt(a*csc(x)^4) - a^3*cos(x)*sqrt(a*csc(x)^4)*sin(x)]
    @test_int [(a*csc(x)^4)^(5/2), x, 3, (-(4/3))*a^2*cos(x)^2*cot(x)*sqrt(a*csc(x)^4) - (6/5)*a^2*cos(x)^2*cot(x)^3*sqrt(a*csc(x)^4) - (4/7)*a^2*cos(x)^2*cot(x)^5*sqrt(a*csc(x)^4) - (1/9)*a^2*cos(x)^2*cot(x)^7*sqrt(a*csc(x)^4) - a^2*cos(x)*sqrt(a*csc(x)^4)*sin(x)]
    @test_int [(a*csc(x)^4)^(3/2), x, 3, (-(2/3))*a*cos(x)^2*cot(x)*sqrt(a*csc(x)^4) - (1/5)*a*cos(x)^2*cot(x)^3*sqrt(a*csc(x)^4) - a*cos(x)*sqrt(a*csc(x)^4)*sin(x)]
    @test_int [(a*csc(x)^4)^(1/2), x, 3, (-cos(x))*sqrt(a*csc(x)^4)*sin(x)]
    @test_int [1/(a*csc(x)^4)^(1/2), x, 3, -(cot(x)/(2*sqrt(a*csc(x)^4))) + (x*csc(x)^2)/(2*sqrt(a*csc(x)^4))]
    @test_int [1/(a*csc(x)^4)^(3/2), x, 5, -((5*cot(x))/(16*a*sqrt(a*csc(x)^4))) + (5*x*csc(x)^2)/(16*a*sqrt(a*csc(x)^4)) - (5*cos(x)*sin(x))/(24*a*sqrt(a*csc(x)^4)) - (cos(x)*sin(x)^3)/(6*a*sqrt(a*csc(x)^4))]
    @test_int [1/(a*csc(x)^4)^(5/2), x, 7, -((63*cot(x))/(256*a^2*sqrt(a*csc(x)^4))) + (63*x*csc(x)^2)/(256*a^2*sqrt(a*csc(x)^4)) - (21*cos(x)*sin(x))/(128*a^2*sqrt(a*csc(x)^4)) - (21*cos(x)*sin(x)^3)/(160*a^2*sqrt(a*csc(x)^4)) - (9*cos(x)*sin(x)^5)/(80*a^2*sqrt(a*csc(x)^4)) - (cos(x)*sin(x)^7)/(10*a^2*sqrt(a*csc(x)^4))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*((b*csc(c+d*x))^p)^n*with*n*symbolic=#

    @test_int [((b*csc(c + d*x))^p)^n, x, 3, (cos(c + d*x)*((b*csc(c + d*x))^p)^n*HypergeometricFunctions._₂F₁(1/2, (1/2)*(1 - n*p), (1/2)*(3 - n*p), sin(c + d*x)^2)*sin(c + d*x))/(d*(1 - n*p)*sqrt(cos(c + d*x)^2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a*(b*csc(c+d*x))^p)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a*(b*csc(c+d*x))^p)^n*with*n*symbolic=#

    @test_int [(a*(b*csc(c + d*x))^p)^n, x, 3, (cos(c + d*x)*(a*(b*csc(c + d*x))^p)^n*HypergeometricFunctions._₂F₁(1/2, (1/2)*(1 - n*p), (1/2)*(3 - n*p), sin(c + d*x)^2)*sin(c + d*x))/(d*(1 - n*p)*sqrt(cos(c + d*x)^2))]

    #= ::Title:: =#
    #=Integrands*of*the*form*(a*csc(e+f*x))^m*(b*Trg(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a*csc(e+f*x))^m*(b*csc(e+f*x))^n=#

    @test_int [(a*csc(e + f*x))^m*(b*csc(e + f*x))^n, x, 3, (a*cos(e + f*x)*(a*csc(e + f*x))^(-1 + m)*(b*csc(e + f*x))^n*HypergeometricFunctions._₂F₁(1/2, (1/2)*(1 - m - n), (1/2)*(3 - m - n), sin(e + f*x)^2))/(f*(1 - m - n)*sqrt(cos(e + f*x)^2))]
    end

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.0 (a csc)^m (b trg)^n.jl


text:
using StaticArrays: MMatrix, MVector
using StructHelpers: @batteries

export Cell, Dataset, SpacegroupType, basis_vectors

"""
    Cell(lattice, positions, types, magmoms=zeros(length(types)))

The basic input data type of `Spglib`.

Lattice parameters `lattice` are given by a ``3×3`` matrix with floating point values,
where ``𝐚``, ``𝐛``, and ``𝐜`` are given as columns.
Fractional atomic positions `positions` are given
by a ``3×N`` matrix with floating point values, where ``N`` is the number of atoms.
Numbers to distinguish atomic species `types` are given by a list of ``N`` integers.
The collinear polarizations `magmoms` only work with `get_symmetry` and are given
as a list of ``N`` floating point values.
"""
struct Cell{N,L,P,T,M}
    lattice::MMatrix{3,3,L,9}
    positions::MMatrix{3,N,P}
    types::MVector{N,T}
    magmoms::MVector{N,M}
end
function Cell(lattice, positions, types, magmoms = zeros(length(types)))
    if lattice isa AbstractVector
        lattice = hcat(lattice...)
    end
    if positions isa AbstractVector
        positions = hcat(positions...)
    end
    N, L, P, T, M =
        length(types), eltype(lattice), eltype(positions), eltype(types), eltype(magmoms)
    return Cell{N,L,P,T,M}(lattice, positions, types, magmoms)
end

@batteries Cell eq = true hash = true

"""
    basis_vectors(cell::Cell)

Return the three basis vectors from `cell`.
"""
function basis_vectors(cell::Cell)
    lattice = cell.lattice
    return lattice[:, 1], lattice[:, 2], lattice[:, 3]
end

# This is an internal function, do not export!
function _expand_cell(cell::Cell)
    @unpack lattice, positions, types, magmoms = cell
    # Reference: https://github.com/mdavezac/spglib.jl/blob/master/src/spglib.jl#L32-L35 and https://github.com/spglib/spglib/blob/444e061/python/spglib/spglib.py#L953-L975
    clattice = Base.cconvert(Matrix{Cdouble}, lattice) |> transpose
    cpositions = Base.cconvert(Matrix{Cdouble}, positions)
    ctypes = Cint[findfirst(isequal(u), unique(types)) for u in types]
    if magmoms !== nothing
        magmoms = Base.cconvert(Vector{Cdouble}, magmoms)
    end
    return Cell(clattice, cpositions, ctypes, magmoms)
end

# This is an internal type, do not export!
struct SpglibSpacegroupType
    number::Cint
    international_short::NTuple{11,Cchar}
    international_full::NTuple{20,Cchar}
    international::NTuple{32,Cchar}
    schoenflies::NTuple{7,Cchar}
    hall_symbol::NTuple{17,Cchar}
    choice::NTuple{6,Cchar}
    pointgroup_international::NTuple{6,Cchar}
    pointgroup_schoenflies::NTuple{4,Cchar}
    arithmetic_crystal_class_number::Cint
    arithmetic_crystal_class_symbol::NTuple{7,Cchar}
end

"""
    SpglibSpacegroupType(number, international_short, international_full, international, schoenflies, hall_symbol, choice, pointgroup_international, pointgroup_schoenflies, arithmetic_crystal_class_number, arithmetic_crystal_class_symbol)

Represent `SpglibSpacegroupType`, see its [official documentation](https://spglib.github.io/spglib/api.html#spg-get-spacegroup-type).
"""
struct SpacegroupType
    number::Int
    international_short::String
    international_full::String
    international::String
    schoenflies::String
    hall_symbol::String
    choice::String
    pointgroup_international::String
    pointgroup_schoenflies::String
    arithmetic_crystal_class_number::Int
    arithmetic_crystal_class_symbol::String
end

struct SpglibError <: Exception
    msg::AbstractString
end

"""
    get_version()

Obtain the version number of `spglib`.

This is the mergence of `spg_get_major_version`, `spg_get_minor_version`, and `spg_get_micro_version` in its C-API.
"""
function get_version()
    major = ccall((:spg_get_major_version, libsymspg), Cint, ())
    minor = ccall((:spg_get_minor_version, libsymspg), Cint, ())
    micro = ccall((:spg_get_micro_version, libsymspg), Cint, ())
    return VersionNumber(major, minor, micro)
end

tuple2matrix(t::NTuple{9}) = hcat(Iterators.partition(t, 3)...)

function rotsFromTuple(rotsTuple::AbstractVector{NTuple{9,Int32}}, nop::Integer)
    r = Array{Int64,3}(undef, 3, 3, nop)
    for i in 1:nop
        r[:, :, i] = reshape(collect(rotsTuple[i]), 3, 3)
    end
    return r
end

function transFromTuple(transTuple::AbstractVector{NTuple{3,Float64}}, nop::Integer)
    t = Matrix{Float64}(undef, 3, nop)
    for i in 1:nop
        t[:, i] = collect(transTuple[i])
    end
    return t
end

const LETTERS = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"

@static if get_version() >= v"1.15.0"
    include("new.jl")
else
    include("old.jl")
end

function Base.show(io::IO, cell::Cell{N}) where {N}
    println(io, "lattice:")
    print(io, " ")
    display(cell.lattice)
    println(io, "$N atomic positions:")
    print(io, " ")
    display(cell.positions)
    println(io, "$N atoms:")
    println(io, " ", cell.types)
end


text:
using Test, YAAD, YAAD.TestUtils
using DiffRules, SpecialFunctions, NaNMath

@testset "broadcast math.jl" begin

# exclude function whose domain is not in [0, 1)
exclusion = Symbol[
    :asec,
    :acosh,
    :acsc,
    :acoth,
    :asecd,
    :acscd,
]

for (mod, name, nargs) in keys(DiffRules.DEFINED_DIFFRULES)
    f_ex_head = Expr(:., mod, QuoteNode(name))

    if nargs == 1
        name in exclusion && continue

        @eval @test gradcheck(x->Broadcast.broadcasted($f_ex_head, x), Variable(rand(10, 10)))
    end
end

@test gradcheck(x->Broadcast.broadcasted(asec, x), Variable(rand(10, 10) .+ 2.0))
@test gradcheck(x->Broadcast.broadcasted(acosh, x), Variable(rand(10, 10) .+ 2.0))
@test gradcheck(x->Broadcast.broadcasted(acsc, x), Variable(rand(10, 10) .+ 2.0))
@test gradcheck(x->Broadcast.broadcasted(acoth, x), Variable(rand(10, 10) .+ 2.0))
@test gradcheck(x->Broadcast.broadcasted(asecd, x), Variable(rand(10, 10) .+ 2.0))
@test gradcheck(x->Broadcast.broadcasted(acscd, x), Variable(rand(10, 10) .+ 2.0))

end # math.jl

gradcheck(x->Broadcast.broadcasted(Base.csch, x), Variable(rand(10, 10)))


text:
###
### **arb_mat.h** -- matrices over the real numbers
###

### Types, macros and constants

### Memory management
arbcall"void arb_mat_init(arb_mat_t mat, slong r, slong c)"
arbcall"void arb_mat_clear(arb_mat_t mat)"
arbcall"slong arb_mat_allocated_bytes(const arb_mat_t x)"
arbcall"void arb_mat_window_init(arb_mat_t window, const arb_mat_t mat, slong r1, slong c1, slong r2, slong c2)"
arbcall"void arb_mat_window_clear(arb_mat_t window)"

### Conversions
arbcall"void arb_mat_set(arb_mat_t dest, const arb_mat_t src)"
#ni arbcall"void arb_mat_set_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src)"
#ni arbcall"void arb_mat_set_round_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src, slong prec)"
#ni arbcall"void arb_mat_set_fmpq_mat(arb_mat_t dest, const fmpq_mat_t src, slong prec)"

### Random generation
#ns arbcall"void arb_mat_randtest(arb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits)"

### Input and output
arbcall"void arb_mat_printd(const arb_mat_t mat, slong digits)"
#ns arbcall"void arb_mat_fprintd(FILE * file, const arb_mat_t mat, slong digits)"

### Comparisons
arbcall"int arb_mat_equal(const arb_mat_t mat1, const arb_mat_t mat2)"
arbcall"int arb_mat_overlaps(const arb_mat_t mat1, const arb_mat_t mat2)"
arbcall"int arb_mat_contains(const arb_mat_t mat1, const arb_mat_t mat2)"
#ni arbcall"int arb_mat_contains_fmpz_mat(const arb_mat_t mat1, const fmpz_mat_t mat2)"
#ni arbcall"int arb_mat_contains_fmpq_mat(const arb_mat_t mat1, const fmpq_mat_t mat2)"
arbcall"int arb_mat_eq(const arb_mat_t mat1, const arb_mat_t mat2)"
arbcall"int arb_mat_ne(const arb_mat_t mat1, const arb_mat_t mat2)"
arbcall"int arb_mat_is_empty(const arb_mat_t mat)"
arbcall"int arb_mat_is_square(const arb_mat_t mat)"
arbcall"int arb_mat_is_exact(const arb_mat_t mat)"
arbcall"int arb_mat_is_zero(const arb_mat_t mat)"
arbcall"int arb_mat_is_finite(const arb_mat_t mat)"
arbcall"int arb_mat_is_triu(const arb_mat_t mat)"
arbcall"int arb_mat_is_tril(const arb_mat_t mat)"
arbcall"int arb_mat_is_diag(const arb_mat_t mat)"

### Special matrices
arbcall"void arb_mat_zero(arb_mat_t mat)"
arbcall"void arb_mat_one(arb_mat_t mat)"
arbcall"void arb_mat_ones(arb_mat_t mat)"
arbcall"void arb_mat_indeterminate(arb_mat_t mat)"
arbcall"void arb_mat_hilbert(arb_mat_t mat, slong prec)"
arbcall"void arb_mat_pascal(arb_mat_t mat, int triangular, slong prec)"
arbcall"void arb_mat_stirling(arb_mat_t mat, int kind, slong prec)"
arbcall"void arb_mat_dct(arb_mat_t mat, int type, slong prec)"

### Transpose
arbcall"void arb_mat_transpose(arb_mat_t dest, const arb_mat_t src)"

### Norms
arbcall"void arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A)"
arbcall"void arb_mat_frobenius_norm(arb_t res, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_bound_frobenius_norm(mag_t res, const arb_mat_t A)"

### Arithmetic
arbcall"void arb_mat_neg(arb_mat_t dest, const arb_mat_t src)"
arbcall"void arb_mat_add(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)"
arbcall"void arb_mat_sub(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)"
arbcall"void arb_mat_mul_classical(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"void arb_mat_mul_threaded(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"void arb_mat_mul_block(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"void arb_mat_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)"
arbcall"void arb_mat_mul_entrywise(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_sqr(arb_mat_t res, const arb_mat_t mat, slong prec)"
arbcall"void arb_mat_pow_ui(arb_mat_t res, const arb_mat_t mat, ulong exp, slong prec)"
#ni arbcall"void _arb_mat_addmul_rad_mag_fast(arb_mat_t C, mag_srcptr A, mag_srcptr B, slong ar, slong ac, slong bc)"
arbcall"void arb_mat_approx_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)"

### Scalar arithmetic
arbcall"void arb_mat_scalar_mul_2exp_si(arb_mat_t B, const arb_mat_t A, slong c)"
arbcall"void arb_mat_scalar_addmul_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec)"
#ni arbcall"void arb_mat_scalar_addmul_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec)"
arbcall"void arb_mat_scalar_addmul_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec)"
arbcall"void arb_mat_scalar_mul_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec)"
#ni arbcall"void arb_mat_scalar_mul_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec)"
arbcall"void arb_mat_scalar_mul_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec)"
arbcall"void arb_mat_scalar_div_si(arb_mat_t B, const arb_mat_t A, slong c, slong prec)"
#ni arbcall"void arb_mat_scalar_div_fmpz(arb_mat_t B, const arb_mat_t A, const fmpz_t c, slong prec)"
arbcall"void arb_mat_scalar_div_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec)"

### Gaussian elimination and solving
arbcall"int arb_mat_lu_classical(slong * perm, arb_mat_t LU, const arb_mat_t A, slong prec)"
arbcall"int arb_mat_lu_recursive(slong * perm, arb_mat_t LU, const arb_mat_t A, slong prec)"
arbcall"int arb_mat_lu(slong * perm, arb_mat_t LU, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_solve_tril_classical(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)"
arbcall"void arb_mat_solve_tril_recursive(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)"
arbcall"void arb_mat_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)"
arbcall"void arb_mat_solve_triu_classical(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)"
arbcall"void arb_mat_solve_triu_recursive(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)"
arbcall"void arb_mat_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)"
arbcall"void arb_mat_solve_lu_precomp(arb_mat_t X, const slong * perm, const arb_mat_t LU, const arb_mat_t B, slong prec)"
arbcall"int arb_mat_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"int arb_mat_solve_lu(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"int arb_mat_solve_precond(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"int arb_mat_solve_preapprox(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, const arb_mat_t R, const arb_mat_t T, slong prec)"
arbcall"int arb_mat_inv(arb_mat_t X, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_det_lu(arb_t det, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_det_precond(arb_t det, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_det(arb_t det, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_approx_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)"
arbcall"void arb_mat_approx_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)"
arbcall"int arb_mat_approx_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_approx_solve_lu_precomp(arb_mat_t X, const slong * perm, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"int arb_mat_approx_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"int arb_mat_approx_inv(arb_mat_t X, const arb_mat_t A, slong prec)"

### Cholesky decomposition and solving
arbcall"int _arb_mat_cholesky_banachiewicz(arb_mat_t A, slong prec)"
arbcall"int arb_mat_cho(arb_mat_t L, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_solve_cho_precomp(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, slong prec)"
arbcall"int arb_mat_spd_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)"
arbcall"void arb_mat_inv_cho_precomp(arb_mat_t X, const arb_mat_t L, slong prec)"
arbcall"int arb_mat_spd_inv(arb_mat_t X, const arb_mat_t A, slong prec)"
arbcall"int _arb_mat_ldl_inplace(arb_mat_t A, slong prec)"
arbcall"int _arb_mat_ldl_golub_and_van_loan(arb_mat_t A, slong prec)"
arbcall"int arb_mat_ldl(arb_mat_t res, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_solve_ldl_precomp(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, slong prec)"
arbcall"void arb_mat_inv_ldl_precomp(arb_mat_t X, const arb_mat_t L, slong prec)"

### Characteristic polynomial and companion matrix
arbcall"void _arb_mat_charpoly(arb_ptr poly, const arb_mat_t mat, slong prec)"
arbcall"void arb_mat_charpoly(arb_poly_t poly, const arb_mat_t mat, slong prec)"
arbcall"void _arb_mat_companion(arb_mat_t mat, arb_srcptr poly, slong prec)"
arbcall"void arb_mat_companion(arb_mat_t mat, const arb_poly_t poly, slong prec)"

### Special functions
arbcall"void arb_mat_exp_taylor_sum(arb_mat_t S, const arb_mat_t A, slong N, slong prec)"
arbcall"void arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)"
arbcall"void arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec)"
arbcall"void _arb_mat_diag_prod(arb_t res, const arb_mat_t mat, slong a, slong b, slong prec)"
arbcall"void arb_mat_diag_prod(arb_t res, const arb_mat_t mat, slong prec)"

### Sparsity structure
#ni arbcall"void arb_mat_entrywise_is_zero(fmpz_mat_t dest, const arb_mat_t src)"
#ni arbcall"void arb_mat_entrywise_not_is_zero(fmpz_mat_t dest, const arb_mat_t src)"
arbcall"slong arb_mat_count_is_zero(const arb_mat_t mat)"
arbcall"slong arb_mat_count_not_is_zero(const arb_mat_t mat)"

### Component and error operations
arbcall"void arb_mat_get_mid(arb_mat_t B, const arb_mat_t A)"
arbcall"void arb_mat_add_error_mag(arb_mat_t mat, const mag_t err)"

### Eigenvalues and eigenvectors


text:
module PkgDownloader

# Flag responsible for printing status info
VERBOSE = true

include("common/errors.jl")
include("utils/utils.jl")

# Write your package code here.
greet() = print("Hello World!")

end

max_stars_repo_name: julbinb/julia-pkgsrc-downloader

max_issues_repo_name: julbinb/julia-pkgsrc-downloader

max_forks_repo_name: julbinb/julia-pkgsrc-downloader


text:
"""
Evaluate an orthogonal polynomial using the three term recurrence relation.
The recurrence relation is assumed to have the form:
'''
    p_{n+1}(x) = (A_n x + B_n) * p_n(x) - C_n * p_{n-1}(x)
    p_{-1} = 0; p_0 = p0
'''
with the coefficients implemented by the `rec_An`, `rec_Bn` and `rec_Cn`
functions and with the initial value implemented with the `p0` function.
This is the convention followed by the DLMF, see `http://dlmf.nist.gov/18.9#i`.
"""
function recurrence_eval(b::OPS, idx::PolynomialDegree, x)
	T = codomaintype(b)
    z0 = T(p0(b))
    z1 = convert(T, rec_An(b, 0) * x + rec_Bn(b, 0))*z0

    d = degree(idx)
    if d == 0
        return z0
    end
    if d == 1
        return z1
    end

    z = z1
    for i = 1:d-1
        z = (rec_An(b, i)*x + rec_Bn(b, i)) * z1 - rec_Cn(b, i) * z0
        z0 = z1
        z1 = z
    end
    z
end

recurrence_eval(b::OPS, idx::LinearIndex, x) = recurrence_eval(b, native_index(b, idx), x)

"Iterate over the values of an orthogonal polynomial sequence using the recurrence relation."
struct OPSValueIterator{T,O<:OPS{T}} <: DictionaryValueIterator{T}
	dict	::	O
	x		::	T
end

OPSValueIterator(O::OPS{T}, x) where {T} = OPSValueIterator(O, convert(T, x))

pointvalues(ops::OPS, x) = OPSValueIterator(ops, x)

# The state vector is (z0,z1,d) as follows:
# - z0 is the value of the polynomial of degree d-1
# - z1 is the value of the polynomial of degree d
function iterate(iter::OPSValueIterator)
	ops = dictionary(iter)
	val = p0(ops)
	if length(ops) > 0
		(1,val), (zero(val),val,0)
	else
		nothing
	end
end

# Here we have to compute the polynomial p_{d+1}(x) using p_d and p_{d-1}
function iterate(iter::OPSValueIterator, state)
	z0, z1, degree = state
	dict = dictionary(iter)
	x = point(iter)
	if degree < length(dict)-1
		if degree == 0
			# We avoid evaluating rec_Cn with degree=0
			z = (rec_An(dict, degree)*x + rec_Bn(dict, degree)) * z1
			(degree+2,z), (z1, z, degree+1)
		else
			z = (rec_An(dict, degree)*x + rec_Bn(dict, degree)) * z1 - rec_Cn(dict, degree) * z0
			(degree+2,z), (z1, z, degree+1)
		end
	else
		nothing
	end
end

"""
Evaluate all the orthonormal polynomials in the sequence in the given
point `x`, using the three-term recurrence relation.

The implementation is based on Gautschi, 2004, Theorem 1.29, p. 12.
"""
function recurrence_eval_orthonormal!(result, b::OPS, x)
    @assert length(result) == length(b)

    α, β = monic_recurrence_coefficients(b)
    # Explicit formula for the constant orthonormal polynomial
    result[1] = 1/sqrt(first_moment(b))
    if length(result) > 1
        # We use an explicit formula for the first degree polynomial too
        result[2] = 1/sqrt(β[2]) * (x-α[1]) * result[1]
        # The rest follows (1.3.13) in the book of Gautschi
        for i = 2:length(b)-1
            result[i+1] = 1/sqrt(β[i+1]) * ( (x-α[i])*result[i] - sqrt(β[i])*result[i-1])
        end
    end
    result
end

"""
Evaluate the derivative of an orthogonal polynomial, based on taking the
derivative of the three-term recurrence relation (see `recurrence_eval`).
"""
function recurrence_eval_derivative(b::OPS, idx::PolynomialDegree, x)
	T = codomaintype(b)
    z0 = one(p0(b))
    z1 = convert(T, rec_An(b, 0) * x + rec_Bn(b, 0))*z0
    z0_d = zero(T)
    z1_d = convert(T, rec_An(b, 0))

    d = degree(idx)
    if d == 0
        return z0_d
    end
    if d == 1
        return z1_d
    end

    z = z1
    z_d = z1_d
    for i = 1:d-1
        z = (rec_An(b, i)*x + rec_Bn(b, i)) * z1 - rec_Cn(b, i) * z0
        z_d = (rec_An(b, i)*x + rec_Bn(b, i)) * z1_d + rec_An(b, i)*z1 - rec_Cn(b, i) * z0_d
        z0 = z1
        z1 = z
        z0_d = z1_d
        z1_d = z_d
    end
    z_d
end

recurrence_eval_derivative(b::OPS, idx::LinearIndex, x) =
    recurrence_eval_derivative(b, native_index(b, idx), x)

"""
Return the recurrence coefficients `α_n` and `β_n` for the monic orthogonal
polynomials (i.e., with leading order coefficient `1`). The recurrence relation
is
```
p_{n+1}(x) = (x-α_n)p_n(x) - β_n p_{n-1}(x).
```
The result is a vector with indices starting at `1` (such that `α_n` above is
given by `α[n+1]`). The last value is `α_{n-1} = α[n]`.
"""
function monic_recurrence_coefficients(b::OPS)
    T = codomaintype(b)
    # n is the maximal degree polynomial
    n = length(b)
    α = zeros(T, n)
    β = zeros(T, n)

    # We keep track of the leading order coefficients of p_{k-1}, p_k and
    # p_{k+1} in γ_km1, γ_k and γ_kp1 respectively. The recurrence relation
    # becomes, with p_k(x) = γ_k q_k(x):
    #
    # γ_{k+1} q_{k+1}(x) = (A_k x + B_k) γ_k q_k(x) - C_k γ_{k-1} q_{k-1}(x).
    #
    # From this we derive the formulas below.

    # We start with k = 1 and recall that p_0(x) = 1 and p_1(x) = A_0 x + B_0,
    # hence γ_0 = 1 and γ_1 = A_0.
    γ_km1 = one(T)
    γ_k = rec_An(b, 0)
    α[1] = -rec_Bn(b, 0) / rec_An(b, 0)
    # see equation (1.3.6) from Gautschi's book, "Orthogonal Polynomials and Computation"
    β[1] = first_moment(b)

    # We can now loop for k from 1 to n.
    for k in 1:n-1
        γ_kp1 = rec_An(b, k) * γ_k
        α[k+1] = -rec_Bn(b, k) * γ_k / γ_kp1
        β[k+1] =  rec_Cn(b, k) * γ_km1 / γ_kp1
        γ_km1 = γ_k
        γ_k = γ_kp1
    end
    α, β
end

"""
Evaluate the three-term recurrence relation for monic orthogonal polynomials
```
p_{n+1}(x) = (x-α_n)p_n(x) - β_n p_{n-1}(x).
```
"""
function monic_recurrence_eval(α, β, idx, x)
    @assert idx > 0
    T = eltype(α)
    # n is the maximal degree polynomial, we want to evaluate p_n(x)
    n = length(α)

    # We store the values p_{k-1}(x), p_k(x) and p_{k+1}(x) in z_km1, z_k and
    # z_kp1 respectively.
    z_km1 = zero(T)
    z_k = one(T)
    z_kp1 = z_k
    if idx == 1
        z_k
    else
        for k in 0:idx-2
            z_kp1 = (x-α[k+1])*z_k - β[k+1]*z_km1
            z_km1 = z_k
            z_k = z_kp1
        end
        z_k
    end
end


text:
function psiNextLebesgueMeasure(p::Vector,psi::Vector)
    N = length(psi)
    psik = zero(psi)
    if length(p)==1
        p = cat(1,zero(p[1]),p)
        c = sqrt(1-abs2(p[end]))
    else
        c = sqrt((1-abs2(p[end]))/(1-abs2(p[end-1])))
    end
    for j=2:(N)
        psik[j] = c*sum([p[end]^(j-i)*(psi[i-1]-conj(p[end-1])*psi[i]) for i=2:j])
    end
    return psik
end

function psiBasis(p::Vector,u::Vector)
    N = length(u)
    B = zeros(Complex128,(N,length(p)+1))
    B[:,1] = u
    for k=1:length(p)
        B[:,k+1] = psiNextLebesgueMeasure([p[1:k]],B[:,k])
    end
    return B[2:end,2:end]
end

function dictObjective(x::Vector,g::Vector,target::Vector,u::Vector)
    p = x[1]*exp(im*2*pi*x[2])
    trial = psiNextLebesgueMeasure([p],u)[2:end]
    return abs2(dot(target,trial))
end

function unitStep(k)
    return k>=0 ? 1 : 0
end

function invBlaschkeElem(p::Complex128,K::Int)
    return [((conj(p)^(-i))*unitStep(i-1)-conj(p)^(-i-1))::Complex128 for i=0:(K-1)]
end

function unitDiskGrid(n::Int,delta::Float64)
    N = nextpow2(n)
    M = N>>1
    v = linspace(delta,one(Float64)-delta,M)
    Zjk = zeros(Complex128,N,N)
    for j=1:(M)
             #bottom, running counter-clockwise
            Zjk[(M+j),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(5*pi/4,7*pi/4,2*j))
            #top, running clockwise
            Zjk[(M-j+1),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(3*pi/4,1*pi/4,2*j))
    end
    for j=1:(M)
            #right, running clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M+j] = v[j]*exp(im*linspace(pi/4,-pi/4,2*j)) #right
            #left, running counter-clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M-j+1] = v[j]*exp(im*linspace(3*pi/4,5*pi/4,2*j)) #left
    end
    return reshape(Zjk,N^2)
end

function afdT(target::Vector,u::Vector,Npoles::Int,Nstarts::Int,lmethod,ltol)
    Zstarts = unitDiskGrid(int(sqrt(Nstarts)),0.3)
    opt = Opt(lmethod,2)
    ftol_rel!(opt, ltol) #local stopping criteria
    lower_bounds!(opt,[0.01; 0])
    upper_bounds!(opt,[0.9;  1-eps(Float64)])
    p = zeros(Complex128,Npoles)
    R = copy(target)
    for j=1:Npoles
        ovec = zeros(Float64,length(Zstarts))
        avec = zeros(Float64,(length(Zstarts),2))
        max_objective!(opt,(x,g)->dictObjective(x,g,R,u))
        for k=1:length(Zstarts)
            beta1 = abs(Zstarts[k]);
            beta2 = (angle(Zstarts[k])+pi)/(2*pi);
            beta0 = [beta1; beta2]
            (ovec[k],avec[k,:],e) = optimize(opt,beta0)
        end
        aopt = avec[indmax(ovec),:]
        p[j] = aopt[1]*exp(im*2*pi*aopt[2])
        B = psiBasis([p[1:j]],u)
        R = (eye(length(target))-B*pinv(B))*target
        Q = invBlaschkeElem(p[j],(length(target))*4)
        R = oneSidedConv(R,Q)[1:length(target)]
    end
    return p
end

function oneSidedConv(u::Vector,v::Vector)
    U = fft(u)
    V = fft(v)
    if length(U)<=length(V)
        return ifft(V.*cat(1,U,zeros(Complex128,length(V)-length(U))))
    else
        return ifft(U.*cat(1,V,zeros(Complex128,length(U)-length(V))))
    end
end


text:
abstract type MathProgCallbackData end
export MathProgCallbackData

# callback has signature:
# function callback(d::MathProgCallbackData)

@define_interface begin
    setlazycallback!
    setcutcallback!
    setheuristiccallback!
    setinfocallback!
    cbgetmipsolution
    cbgetlpsolution
    cbgetobj
    cbgetbestbound
    cbgetexplorednodes
    cbgetstate
    cbaddcut!
    cbaddcutlocal!
    cbaddlazy!
    cbaddlazylocal!
    cbaddsolution!
    cbsetsolutionvalue!
end


text:
##############################################################
#          single quadruplet (item, examinee, response, w)
##############################################################

function _posterior(
    parameters::Parameters2PL,
    latent::Latent1D,
    response::Response,
    w::PolyaGammaSample,
)
    a = parameters.a
    b = parameters.b
    latent_val = latent.val
    r_val = response.val
    w_val = w.val

    sigma2 = [(latent_val - b)^2 * w.val, a^2 * w_val]
    sigma2 = 1 ./ (sigma2 .+ (1 ./ Distributions.var.(parameters.prior.v)))
    mu = [(latent_val - b) * (r_val - 0.5), -a * ((r_val - 0.5) - (a * latent_val * w_val))]
    mu =
        sigma2 .* (
            mu + (
                Distributions.mean.(parameters.prior.v) ./
                Distributions.var.(parameters.prior.v)
            )
        )
    return Distributions.Product([
        Distributions.TruncatedNormal(mu[1], sqrt(sigma2[1]), 0.0, Inf),
        Distributions.Normal(mu[2], sqrt(sigma2[2])),
    ])
end

function posterior(
    item::AbstractItem,
    examinee::AbstractExaminee,
    response::AbstractResponse,
    w::PolyaGammaSample,
)
    return _posterior(item.parameters, examinee.latent, response, w)
end

function __posterior(
    parameters::Parameters2PL,
    latents::Vector{Latent1D},
    responses::Vector{Float64},
    W::Vector{Float64},
)
    prior = parameters.prior.v
    a = parameters.a
    b = parameters.b
    sigma_2_prior = 1 ./ Distributions.var.(prior)
    #W = map( e -> Distributions.rand(PolyaGammaDevRoye1Sampler(1.0, item.parameters.a *(e.latent.val - item.parameters.b))), examinees)
    sigma2 = mapreduce((e, w) -> [(e.val - b)^2, a^2] .* w, +, latents, W)
    sigma2 = 1 ./ (sigma2 + sigma_2_prior)
    mu = mapreduce(
        (e, w, r) -> [
            (e.val - b) * (r - 0.5),
            - a * ((r - 0.5) - (a * e.val * w))
        ],
        +,
        latents,
        W,
        responses,
    )
    mu = sigma2 .* (mu + (Distributions.mean.(prior) .* sigma_2_prior))
    return Distributions.Product([
        Distributions.TruncatedNormal(mu[1], sqrt(sigma2[1]), 0.0, Inf),
        Distributions.Normal(mu[2], sqrt(sigma2[2])),
    ])
end

function _posterior(
    parameters::Parameters2PL,
    latents::Vector{Latent1D}, #must be sorted by e.idx
    responses::Vector{Response}, #only responses of item sorted by e.idx
    W::Vector{PolyaGammaSample}, #sorted by e.idx
)
    prior = parameters.prior.v
    a = parameters.a
    b = parameters.b
    sigma_2_prior = 1 ./ Distributions.var.(prior)
    #W = map( e -> Distributions.rand(PolyaGammaDevRoye1Sampler(1.0, item.parameters.a *(e.latent.val - item.parameters.b))), examinees)
    sigma2 = mapreduce((e, w) -> [(e.val - b)^2, a^2] .* w.val, +, latents, W)
    sigma2 = 1 ./ (sigma2 + sigma_2_prior)
    mu = mapreduce(
        (e, w, r) -> [
            (e.val - b) * (r.val - 0.5),
            -a * ((r.val - 0.5) - (a * e.val * w.val))
        ],
        +,
        latents,
        W,
        responses,
    )
    mu = sigma2 .* (mu + (Distributions.mean.(prior) .* sigma_2_prior))
    return Distributions.Product([
        Distributions.TruncatedNormal(mu[1], sqrt(sigma2[1]), 0.0, Inf),
        Distributions.Normal(mu[2], sqrt(sigma2[2])),
    ])
end

# function __posterior_2pl_1d(
#     a::Float64,
#     b::Float64, #must be sorted by e.idx
#     mu_prior::Vector{Float64},
#     _1_sigma_2_prior::Vector{Float64},
#     latents_val::Vector{Float64},
#     responses_val::Vector{Float64}, #only responses of item sorted by e.idx
#     w_val::Vector{Float64}) #sorted by e.idx
#     sigma2 = 1 ./ (mapreduce((l, w) -> [(l - b)^2, a^2] .* w, +, latents_val, w_val) + _1_sigma_2_prior)
#     mu = mapreduce(
#         (l, w, r) -> [
#             (l - b) * (r - 0.5),
#             -a * ((r - 0.5) - (a * l * w))
#         ],
#         +,
#         latents_val,
#         w_val,
#         responses_val,
#     )
#     mu = sigma2 .* (mu + (mu_prior .* _1_sigma_2_prior))
#     return Distributions.Product([
#         Distributions.TruncatedNormal(mu[1], sqrt(sigma2[1]), 0.0, Inf),
#         Distributions.Normal(mu[2], sqrt(sigma2[2])),
#     ])
# end

# function _posterior(
#     parameters::Parameters2PL,
#     latents::Vector{Latent1D}, #must be sorted by e.idx
#     responses::Vector{Response}, #only responses of item sorted by e.idx
#     W::Vector{PolyaGammaSample}, #sorted by e.idx
# )
#     return __posterior_2pl_1d(
#         parameters.a,
#         parameters.b,
#         Distributions.mean.(parameters.prior.v),
#         1 ./ Distributions.var.(parameters.prior.v),
#         map(l -> l.val, latents),
#         map(r -> r.val, responses),
#         map(w -> w.val, W)
#         )
# end

function posterior(
    item::AbstractItem,
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{<:AbstractResponse},
    W::Vector{PolyaGammaSample};
    already_sorted = false
)
    if !already_sorted
        sort!(examinees, by = e -> e.idx)
        sort!(responses, by = r -> r.examinee_idx)
    end
    return _posterior(item.parameters, map(e -> e.latent, examinees), responses, W)
end

function update_posterior!(
    item::AbstractItem,
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{<:AbstractResponse},
    W::Vector{PolyaGammaSample};
    already_sorted = false,
)
    if !already_sorted
        sort!(examinees, by = e -> e.idx)
        sort!(responses, by = r -> r.examinee_idx)
    end
    item.parameters.posterior = posterior(item, examinees, responses, W; already_sorted = true)
end


text:
include("@aml_create/get_arg_xml_.jl")
include("@aml_create/get_struct_xml_.jl")

"""
@aml creator function
"""
function aml_create(expr::Expr, args_param, args_defaultvalue, args_type, args_var, args_name, args_function, args_literaltype, struct_name, struct_nodetype, struct_function, is_struct_mutable, args_custom_creator, args_custom_extractor, args_custom_updater, T)

    expr.head == :struct || error("Invalid usage of aml_create")

    # defining outter constructors
    # Only define a constructor if the type has fields, otherwise we'll get a stack
    # overflow on construction
    if !isempty(args_var)

        ################################################################
        # Arguments methods

        # Non-aml arguments are ignored
        isaml_args = (!).(ismissing.(args_name)) # missing arg_name means not aml

        amlargs_name = args_name[isaml_args]
        amlargs_function = args_function[isaml_args]
        amlargs_var = args_var[isaml_args]
        amlargs_param = args_param[isaml_args]
        amlargs_defaultvalue = args_defaultvalue[isaml_args]
        amlargs_literaltype = args_literaltype[isaml_args]
        amlargs_type  = args_type[isaml_args]
        push!(isaml_args, false) # separating end
        amlargs_custom_creator = args_custom_creator[isaml_args]
        custom_creator_end = args_custom_creator[end]
        amlargs_custom_extractor = args_custom_extractor[isaml_args]
        custom_extractor_end = args_custom_extractor[end]
        amlargs_custom_updater = args_custom_updater[isaml_args]
        custom_updater_end = args_custom_updater[end]

        amlargs_num = length(amlargs_var)

        args_xmlcreator=Vector{Expr}(undef,amlargs_num)
        args_xmlextractor=Vector{Expr}(undef,amlargs_num)
        args_xmlupdater=Vector{Expr}(undef,amlargs_num)

        args_varcall = Vector{Expr}(undef,amlargs_num)

        ##########################
        # Each argument of the struct
        for iArg=1:amlargs_num
            argtype = amlargs_type[iArg]
            argvar = amlargs_var[iArg]
            argname = amlargs_name[iArg]
            argliteraltype = amlargs_literaltype[iArg]
            argfunction=amlargs_function[iArg]
            argcustomcreator = amlargs_custom_creator[iArg]
            argcustomextractor = amlargs_custom_extractor[iArg]
            argcustomupdater = amlargs_custom_updater[iArg]
            argsym=QuoteNode(argvar)
            ##########################
            # call Expr - For mutability
            args_varcall[iArg] = :(str.$argvar)
            argvarcall = args_varcall[iArg]
            ##########################
            has_arg_xmlchecker = !ismissing(argfunction)

            inps = (has_arg_xmlchecker, esc(argtype), esc(argvar), argname, argliteraltype, esc(argfunction), argsym, argvarcall)

            args_xmlcreator[iArg]=get_arg_xmlcreator(argcustomcreator, inps...)

            args_xmlextractor[iArg]=get_arg_xmlextractor(argcustomextractor, inps...)

            if is_struct_mutable
                args_xmlupdater[iArg] = get_arg_xmlupdater(argcustomupdater, inps...)
            end

        end # endfor
        ################################################################
        # Remove supertype if it is present
        if T isa Expr && T.head == :<:
            T = T.args[1]
        end
        ################################################################
        # Type name is a single name (symbol)
        if T isa Symbol
            S = T
            iscurly = false
        ################################################################
        # Parametric type structs
        elseif T isa Expr && T.head == :curly
            # if T == S{A<:AA,B<:BB}, define two methods
            #   S(...) = ...
            #   S{A,B}(...) where {A<:AA,B<:BB} = ...
            S = T.args[1]
            P = T.args[2:end]
            Q = [U isa Expr && U.head == :<: ? U.args[1] : U for U in P]
            SQ = :($S{$(Q...)})

            iscurly = true
        ################################################################
        else
            error("Invalid usage of @aml")
        end
        ################################################################
        node_initializer = :( aml = createnode($struct_nodetype, $struct_name) )

        struct_definition =:($expr)

        # args_var is used here to support all the inputs (even non-aml)
        esc_args_var = esc.(args_var)

        struct_xmlchecker = get_struct_xmlchecker(struct_function, esc_args_var)

        # Creator
        struct_xmlcreator = get_struct_xmlcreator(S, amlargs_param, struct_xmlchecker, node_initializer, args_xmlcreator, esc_args_var, custom_creator_end)
        # Extractor
        struct_xmlextractor = get_struct_xmlextractor(S, args_xmlextractor, struct_xmlchecker, esc_args_var, custom_extractor_end)

        if iscurly
            # Creator
            struct_xmlcreator_curly = get_struct_xmlcreator(SQ, P, amlargs_param, struct_xmlchecker, node_initializer, args_xmlcreator, esc_args_var, custom_creator_end)
            # Extractor
            struct_xmlextractor_curly = get_struct_xmlextractor(SQ, P, args_xmlextractor, struct_xmlchecker, esc_args_var, custom_extractor_end)
        else
            struct_xmlcreator_curly = nothing
            struct_xmlextractor_curly = nothing
        end

        nothing_method = :( ($(esc(S)))(::Nothing) = nothing )
        # convertNothingMethod = :(Base.convert(::Type{($(esc(S)))}, ::Nothing) = nothing) # for passing nothing to function without using Union{Nothing, S} in the definition
        self_method = :( ($(esc(S)))(in::$(esc(S))) = $(esc(S))(in.aml) )
        pprint_method = :( AcuteML.pprint(in::$(esc(S))) = pprint(in.aml) )

        struct_xmlupdater = get_struct_xmlupdater(is_struct_mutable, S, args_xmlupdater, struct_function, args_varcall, custom_updater_end)

        out = quote
            Base.@__doc__($(esc(struct_definition)))
            $struct_xmlcreator
            $struct_xmlcreator_curly
            $struct_xmlextractor
            $struct_xmlextractor_curly
            $nothing_method
            # $convertNothingMethod
            $self_method
            $pprint_method
            $struct_xmlupdater
        end
    else
        out = nothing
    end
    return out
end

max_stars_repo_name: aminya/AML


text:
"""

    Interval unions sets of defined by unions of disjoint intervals.
    This file includes constructors, arithmetic (including intervals and scalars)
    and complement functions

    Empty sets and intersecting intervals are appropriately handled in the constructor:

    julia> a = interval(0,2) ∪ interval(3,4)
    [0, 2] ∪ [3, 4]

    julia> b = interval(1,2) ∪ interval(4,5) ∪ ∅
    [1, 2] ∪ [4, 5]

    julia> c = a * b
    [0, 10] ∪ [12, 20]

    julia> complement(c)
    [-∞, 0] ∪ [10, 12] ∪ [20, ∞]

"""

###
#   IntervalUnion constructor. Consists of a vector of intervals
###
struct IntervalUnion{T<:Real} <: AbstractInterval{T}
    v :: Vector{Interval{T}}
end

###
#   Outer constructors
###
function intervalUnion(v)
    x = IntervalUnion(v)
    sort!(x.v)
    remove_empties!(x)
    condense!(x)
    closeGaps!(x)
    return x
end

intervalUnion(num :: Real) = IntervalUnion([interval(num)])

intervalUnion(lo :: Real, hi :: Real) = IntervalUnion([interval(lo,hi)])

intervalUnion(x :: Interval) = intervalUnion([x])
∪(x :: Interval) = intervalUnion(x)

∪(x :: Interval, y :: Interval) = intervalUnion([x; y])
∪(x :: Array{Interval{T}}) where T <:Real = intervalUnion(x)

intervalUnion(x :: Interval, y :: IntervalUnion) = intervalUnion([x; y.v])
∪(x :: Interval, y :: IntervalUnion) = intervalUnion(x,y)

intervalUnion(x :: IntervalUnion, y :: Interval) = intervalUnion([x.v; y])
∪(x :: IntervalUnion, y :: Interval) = intervalUnion(x,y)

intervalUnion(x :: IntervalUnion, y :: IntervalUnion) = intervalUnion([x.v; y.v])
∪(x :: IntervalUnion, y :: IntervalUnion) = intervalUnion(x,y)

# MultiInterval can act like a vector
getindex(x :: IntervalUnion, ind :: Integer) = getindex(x.v,ind)
getindex(x :: IntervalUnion, ind :: Array{ <: Integer}) = getindex(x.v,ind)

# Remove ∅ from IntervalUnion
function remove_empties!(x :: IntervalUnion)
    deleteat!(x.v, (i for i in eachindex(x.v) if isempty(x.v[i])))
    x
end

function closeGaps!(x :: IntervalUnion, maxInts = MAXINTS[1])

    while length(x.v) > maxInts     # Global

        # Complement code
        v = sort(x.v)
        vLo = left.(v)
        vHi = right.(v)
        vLo[1] == -∞ ? popfirst!(vLo) : pushfirst!(vHi, -∞)
        vHi[end] == ∞ ? pop!(vHi) : push!(vLo, ∞)
        xc = interval.(vHi,vLo)

        # Close smallest width of complement
        widths = diam.(xc)
        d, i = findmin(widths)
        merge = hull(x[i-1], x[i])
        popat!(x.v, i)
        popat!(x.v, i-1)
        push!(x.v, merge)
        sort!(x.v)
    end
    return x
end

function condense!(f::Function, x :: IntervalUnion)
    v = x.v
    sort!(v)
    i = 1
    while i < length(v)
        these = (i - 1) .+ findall(collect(!f(v[i] ∩ v[j]) for j in i:length(v)))
        if length(these) > 1
            v[i] = hull(v[these])
            deleteat!(v, @view(these[2:end]))
        end
        i += 1
    end
    return x
end

condense(f, x) = condense!(f, copy(x))

# join intervals that intersect.
condense!(x::IntervalUnion) = condense!(isempty, x)
# join intervals that intersect.
condense(x::IntervalUnion) = condense(isempty, x)
# Join intervals which intersect in more than one point
condense_weak(x::IntervalUnion) = condense(i->isempty(i) || isthin(i), x)

iscondensed(f::Function, x :: IntervalUnion) = all(f(intersect(x.v[i], x.v[j])) for i in eachindex(x.v) for j in 1:i-1)

iscondensed(x::IntervalUnion) = iscondensed(isempty, x)
iscondensed_weak(x::IntervalUnion) = iscondensed(i->isempty(i) || isthin(i), x)