text:
module SimpleFloatFunctions

export midpoint, # an artful midpoint
       modulo, nint,
       square, cube, invsquare, invcube, invsqrt, invcbrt, spread, tld, sld,

const SysFloat = Union{Float64, Float32, Float16}

include("hasnan.jl") # does a vector or an array contain any NaNs

function midpoint(x::T,y::T) where {T}
    x, y = ifelse(abs(y) > abs(x), (y, x), (x, y))
    return x - (x/2 - y/2)
end

#=
   This well-behaved bounded modulo implementation is from
   The pitfalls of verifying floating-point computations
   by David Monniaux, 2008
   http://arxiv.org/abs/cs/0701192v5
=#
function modulo(a::T, lowerbound::T, upperbound::T) where {T}
    delta = upperbound - lowerbound
    a - (floor((a - lowerbound)/delta) * delta)
end

nint(a::T) where {T<:AbstractFloat} = (trunc(a + copysign(one(T)/2,a)))
nint(::Type{I}, a::T) where {T<AbstractFloat, I<:Integer} = I(nint(a))

"""
    square(x)

squares x
"""
function square(x::T) where T<:SysFloat
    return x * x
end

function square(x::T) where T<:Number
    return x * x
end

"""
    cube(x)

cubes x
"""
function cube(x::T) where T<:SysFloat
    hi, lo = one_square(x)
    hi, md = two_prod(hi, x)
    md += lo * x
    hi += md
    return hi
end

function cube(x::T) where T<:Number
    result = x*x*x
    return result
end

"""
    invsquare(x)

1 / x^2
"""
@inline function invsquare(x::T) where T<:Number
    return square(inv(x))
end

"""
    invcube(x)

1 / x^3
"""
@inline function invcube(x::T) where T<:Number
    return cube(inv(x))
end

"""
    invsqrt(x)

1 / x^(1/2)
"""
@inline function invsqrt(x::T) where T<:Number
    return sqrt(inv(x))
end

"""
    invcbrt(x)

1 / x^(1/3)
"""
@inline function invcbrt(x::T) where T<:Number
    return cbrt(inv(x))
end

"""
    spread(x)

spread complements trunc()
the nearest integer to x, away from zero
"""
function spread(x::T) where T<:AbstractFloat
    return signbit(x) ? floor(x) : ceil(x)
end

"""
    tld(x, y)

truncates @ref(trunc) the result of x/y

like cld @ref(cld), fld @ref(fld)
"""
function tld(x::T) where T<:AbstractFloat
    return trunc(x/y)
end

"""
    sld(x, y)

spreads @ref(spred) the result of x/y

like cld @ref(cld), fld @ref(fld)
"""
function sld(x::T) where T<:AbstractFloat
    return spread(x/y)
end

#=
    errorfree transformations (internal use)
=#

"""
    two_sum_hilo(a, b)
*unchecked* requirement `|a| ≥ |b|`
Computes `s = fl(a+b)` and `e = err(a+b)`.
"""
@inline function two_sum_hilo(a::T, b::T) where T<:SysFloat
    s = a + b
    e = b - (s - a)
    return s, e
end

"""
    two_sum(a, b)
Computes `s = fl(a+b)` and `e = err(a+b)`.
"""
@inline function two_sum(a::T, b::T) where T<:SysFloat
    s = a + b
    v = s - a
    e = (a - (s - v)) + (b - v)
    return s, e
end

"""
    two_diff_hilo(a, b)
    t = lo + t
247

*unchecked* requirement `|a| ≥ |b|`
Computes `s = fl(a-b)` and `e = err(a-b)`.
"""
@inline function two_diff_hilo(a::T, b::T) where T<:SysFloat
    s = a - b
    e = (a - s) - b
    s, e
end

"""
    two_diff(a, b)
Computes `s = fl(a-b)` and `e = err(a-b)`.
"""
@inline function two_diff(a::T, b::T) where T<:SysFloat
    s = a - b
    v = s - a
    e = (a - (s - v)) - (b + v)

    s, e
end

"""
    two_prod(a, b)
Computes `s = fl(a*b)` and `e = err(a*b)`.
"""
@inline function two_prod(a::T, b::T) where T<:SysFloat
    p = a * b
    e = fma(a, b, -p)
    p, e
end
cube
"""
    one_square(a)

Computes `s = fl(a*a)` and `e = err(a*a)`.
"""
@inline function one_square(a::T) where T<:SysFloat
    p = a * a
    e = fma(a, a, -p)
    return p, e
end

"""
    one_cube(a)

Computes `s = fl(a*a*a)` and `e = err(a*a*a)`.
"""
@inline function one_cube(a::T) where T<:SysFloat
    hi, lo = one_square(a)
    hihi, hilo = two_prod(hi, a)
    lohi, lolo = two_prod(lo, a)
    hilo, lohi = two_sum_hilo(hilo, lohi)
    hi, lo = two_sum_hilo(hihi, hilo)
    lo += lohi + lolo
    return hi, lo
end

"""
    one_cube_hi(a)

Computes `s = fl(a*a*a)`.
"""
@inline function one_cube_hi(a::T) where T<:SysFloat
    hi, lo = one_square(a)
    hi, md = two_prod(hi, a)
    md += lo * a
    hi += md
    return hi
end

# like doubledouble

# Algorithm 6 from Tight and rigourous error bounds for basic building blocks of double-word arithmetic
function add_dd_dd(xhi::T, xlo::T, yhi::T, ylo::T) where T<:SysFloat
    hi, lo = two_sum(xhi, yhi)
    thi, tlo = two_sum(xlo, ylo)
    c = lo + thi
    hi, lo = two_sum_hilo(hi, c)
    c = tlo + lo
    hi, lo = two_sum_hilo(hi, c)
    return hi, lo
end

function prod_dd_fl(ahi, alo, b)
    hi, lo = two_prod(ahi, b)
    lo += alo * b
    hi, lo = two_sum_hilo(hi, lo)
    return hi, lo
end

function prod_dd_fl_hi(ahi, alo, b)
    hi, lo = two_prod(ahi, b)
    lo += alo * b
    hi += lo
    return hi
end

# Algorithm 12 from Tight and rigourous error bounds for basic building blocks of double-word arithmetic
function prod_dd_dd(xhi::T, xlo::T, yhi::T, ylo::T) where T<:SysFloat
    hi, lo = two_prod(xhi, yhi)
    t = xlo * ylo
    t = fma(xhi, ylo, t)
    t = lo + t
    t = fma(xlo, yhi, t)
    t = lo + t
    hi, lo = two_sum_hilo(hi, t)
    return hi, lo
end

function prod_dd_dd_hi(xhi::T, xlo::T, yhi::T, ylo::T) where T<:SysFloat
    hi, lo = two_prod(xhi, yhi)
    t = xlo * ylo
    t = fma(xhi, ylo, t)
    t = fma(xlo, yhi, t)
    hi += lo + t
    return hi
end

end # module SimpleFloatFunctions


text:
export NonEquilSchwingerKeldyshContour,
    NonEquilTiltedSchwingerKeldyshContour,
    SchwingerKeldyshContour,
    SchwingerKeldyshContourRounded1,
    SchwingerKeldyshContourRounded2,
    SchwingerKeldyshContourRounded3,
    get_timepoint,
    get_derivative_timepoint,
    EucledianContour

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

abstract type AbstractContour end
abstract type AbstractSchwingerKeldyshContour <: AbstractContour end

abstract type AbstractNonEquilSchwingerKeldyshContour <: AbstractContour end

####################################################
###         DEFINING the contour objects
####################################################

################ Eucledian ###################

struct EucledianContour <: AbstractContour
    β
end

"""
    Calculate the time point on the contour from the parametetric variable u in [0,1]
    for the Schwinger-Keldysh contour
"""
function get_timepoint(u,C::EucledianContour)
    return -u * C.β * im
end

############### Non-Equilibrium Schwinger-Keldysh contour ###########
struct NonEquilSchwingerKeldyshContour <: AbstractNonEquilSchwingerKeldyshContour
    tmax::Float32
    ϵ::Float64
end

Base.string(c::NonEquilSchwingerKeldyshContour) = string("tmax: ", c.tmax, "\n", "ϵ: ", c.ϵ)

L1(contour::NonEquilSchwingerKeldyshContour) = contour.tmax
L2(contour::NonEquilSchwingerKeldyshContour) = contour.tmax
T(contour::AbstractNonEquilSchwingerKeldyshContour) = L1(contour) + L2(contour)
t1(contour::AbstractNonEquilSchwingerKeldyshContour) = L1(contour) / T(contour)
t2(contour::AbstractNonEquilSchwingerKeldyshContour) = L2(contour) / T(contour)

############### Non-Equilibrium tilted Schwinger-Keldysh contour ###########
struct NonEquilTiltedSchwingerKeldyshContour <: AbstractNonEquilSchwingerKeldyshContour
    tmax::Float32
    θ::Float64      # Tilt angel
end

Base.string(c::NonEquilTiltedSchwingerKeldyshContour) = print("tmax: ", c.tmax, "\n", "θ: ", c.θ)

L1(contour::NonEquilTiltedSchwingerKeldyshContour) = contour.tmax/cos(contour.θ)
L2(contour::NonEquilTiltedSchwingerKeldyshContour) = contour.tmax/cos(contour.θ)

############### Schwinger-Keldysh contour ###############

struct SchwingerKeldyshContour <: AbstractSchwingerKeldyshContour
    β::Float32
    tmax::Float32
    A_β::Float64
    F_β::Float64
end

Base.string(c::SchwingerKeldyshContour) = string("tmax: ", c.tmax, "\n", "β: ", c.β, "\n", "A_β: ", c.A_β, "\n", "F_β: ", c.F_β)

# Defining functions related to the calculation of the SchwingerKeldyshContour
β1(contour::AbstractSchwingerKeldyshContour) = contour.A_β * contour.F_β * contour.β
β2(contour::AbstractSchwingerKeldyshContour) = (1-contour.A_β) * contour.F_β * contour.β
L1(contour::SchwingerKeldyshContour) = sqrt(β1(contour)^2 + contour.tmax^2)
L2(contour::SchwingerKeldyshContour) = sqrt(β2(contour)^2 + contour.tmax^2)
L3(contour::SchwingerKeldyshContour) = (1-contour.F_β)*contour.β
T(contour::SchwingerKeldyshContour)  = L1(contour) + L2(contour) + L3(contour)
t1(contour::SchwingerKeldyshContour) = L1(contour) / T(contour)
t2(contour::SchwingerKeldyshContour) = L2(contour) / T(contour)
t3(contour::SchwingerKeldyshContour) = L3(contour) / T(contour)

# Default parameters
SchwingerKeldyshContour() =
    SchwingerKeldyshContour(1.0,0.8,0.01,0.5)

############### Schwinger-Keldysh contour Rounded 1 ###############
#
#   The Schwinger-Keldysh contour for one rounded corner, i.e.,
#   rounded corner at the $t_{max}$ point.
#
#
#   NR_RT_Points, NR_ET_Points, β, t_max, A_β and F_β is the same parameters
#   as in normal Schwinger-Keldysh contour
#
#   NR_Rounded1_Points: nr of points in the arc
#   r1: radius of the rounded corner
#
#
struct SchwingerKeldyshContourRounded1 <: AbstractSchwingerKeldyshContour
    β::Float32
    tmax::Float32
    A_β::Float64
    F_β::Float64
    r1::Float64
end

Base.string(c::SchwingerKeldyshContourRounded1) = string("tmax: ", c.tmax, "\n", "β: ", c.β, "\n", "A_β: ", x.A_β, "\n", "F_β:", c.F_β, "\n", "r1:", c.r1)

# Defining functions related to the calculation of the SchwingerKeldyshContour
#β1(contour::SchwingerKeldyshContour) = contour.A_β * contour.F_β * contour.β
#β2(contour::SchwingerKeldyshContourRounded1) = (1-contour.A_β) * contour.F_β * contour.β
β3(contour::AbstractSchwingerKeldyshContour) = (1-contour.F_β) * contour.β

H1(contour::SchwingerKeldyshContourRounded1) = sqrt((contour.tmax-contour.r1)^2 + β1(contour)^2)
L1(contour::SchwingerKeldyshContourRounded1) = sqrt(H1(contour)^2 - contour.r1^2)
γ1(contour::SchwingerKeldyshContourRounded1) = pi/2 - asin(contour.r1/H1(contour)) - acos(β1(contour)/H1(contour))

H2(contour::SchwingerKeldyshContourRounded1) = sqrt((contour.tmax-contour.r1)^2 + β2(contour)^2)
L2(contour::SchwingerKeldyshContourRounded1) = sqrt(H2(contour)^2 - contour.r1^2)
γ2(contour::SchwingerKeldyshContourRounded1) = pi/2 - asin(contour.r1/H2(contour)) - acos(β2(contour)/H2(contour))

Gy1(contour::SchwingerKeldyshContourRounded1) = L1(contour)*sin(γ1(contour))
Gx1(contour::SchwingerKeldyshContourRounded1) = L1(contour)*cos(γ1(contour))
Gy2(contour::SchwingerKeldyshContourRounded1) = L2(contour)*sin(γ2(contour))
Gx2(contour::SchwingerKeldyshContourRounded1) = L2(contour)*cos(γ2(contour))

η1(contour::SchwingerKeldyshContourRounded1) = pi/2 - γ1(contour)
η2(contour::SchwingerKeldyshContourRounded1) = pi/2 - γ2(contour)

CC1(contour::SchwingerKeldyshContourRounded1) = η1(contour)*contour.r1
CC2(contour::SchwingerKeldyshContourRounded1) = η2(contour)*contour.r1
L3(contour::SchwingerKeldyshContourRounded1) = (1-contour.F_β)*contour.β

T(contour::SchwingerKeldyshContourRounded1)  = L1(contour) + CC1(contour) + CC2(contour) + L2(contour) + L3(contour)
t1(contour::SchwingerKeldyshContourRounded1) = L1(contour) / T(contour)
t2(contour::SchwingerKeldyshContourRounded1) = CC1(contour) / T(contour)
t3(contour::SchwingerKeldyshContourRounded1) = CC2(contour) / T(contour)
t4(contour::SchwingerKeldyshContourRounded1) = L2(contour) / T(contour)
t5(contour::SchwingerKeldyshContourRounded1) = L3(contour) / T(contour)

############### Schwinger-Keldysh contour Rounded 2 ###############
#
#   The Schwinger-Keldysh contour for one rounded corner, i.e.,
#   rounded corner at the $t_{max}$ point.
#
#
#   NR_RT_Points, NR_ET_Points, β, t_max, A_β and F_β is the same parameters
#   as in normal Schwinger-Keldysh contour
#
#   NR_Rounded1_Points: nr of points in the arc
#   r1: radius of the rounded corner
#
#
struct SchwingerKeldyshContourRounded2 <: AbstractSchwingerKeldyshContour
    β::Float32
    tmax::Float32
    A_β::Float64
    F_β::Float64
    r1::Float64
    r2::Float64
end

Base.string(c::SchwingerKeldyshContourRounded2) = string("tmax: ", c.tmax, "\n", "β: ", c.β, "\n", "A_β: ", x.A_β, "\n", "F_β:", c.F_β, "\n",
                                                              "r1:", c.r1,"\n", "r2:", c.r2)

# Defining functions related to the calculation of the SchwingerKeldyshContour
#β1(contour::SchwingerKeldyshContourRounded2) = contour.A_β * contour.F_β * contour.β
#β2(contour::SchwingerKeldyshContourRounded2) = (1-contour.A_β) * contour.F_β * contour.β
#β3(contour::SchwingerKeldyshContourRounded2) = (1-contour.F_β) * contour.β

H1(contour::SchwingerKeldyshContourRounded2) = sqrt((contour.tmax-contour.r1)^2 + β1(contour)^2)
L1(contour::SchwingerKeldyshContourRounded2) = sqrt(H1(contour)^2 - contour.r1^2)
γ1(contour::SchwingerKeldyshContourRounded2) = pi/2 - asin(contour.r1/H1(contour)) - acos(β1(contour)/H1(contour))

H2(contour::SchwingerKeldyshContourRounded2) = sqrt((contour.tmax-contour.r1-contour.r2)^2 + β2(contour)^2)
L2(contour::SchwingerKeldyshContourRounded2) = sqrt(H2(contour)^2 - (contour.r1+contour.r2)^2)
γ2(contour::SchwingerKeldyshContourRounded2) = pi/2 - asin((contour.r1+contour.r2)/H2(contour)) - acos(β2(contour)/H2(contour))

Gy1(contour::SchwingerKeldyshContourRounded2) = L1(contour)*sin(γ1(contour))
Gx1(contour::SchwingerKeldyshContourRounded2) = L1(contour)*cos(γ1(contour))
Gy2(contour::SchwingerKeldyshContourRounded2) = L2(contour)*sin(γ2(contour))
Gx2(contour::SchwingerKeldyshContourRounded2) = L2(contour)*cos(γ2(contour))

η1(contour::SchwingerKeldyshContourRounded2) = pi/2 - γ1(contour)
η2(contour::SchwingerKeldyshContourRounded2) = pi/2 - γ2(contour)

CC1_1(contour::SchwingerKeldyshContourRounded2) = η1(contour)*contour.r1
CC1_2(contour::SchwingerKeldyshContourRounded2) = η2(contour)*contour.r1
CC2(contour::SchwingerKeldyshContourRounded2) = η2(contour)*contour.r2
L3(contour::SchwingerKeldyshContourRounded2) = β3(contour)

T(contour::SchwingerKeldyshContourRounded2)  = L1(contour) + CC1_1(contour) + CC1_2(contour) + L2(contour) + CC2(contour) + L3(contour)
t1(contour::SchwingerKeldyshContourRounded2) = L1(contour) / T(contour)
t2(contour::SchwingerKeldyshContourRounded2) = CC1_1(contour) / T(contour)
t3(contour::SchwingerKeldyshContourRounded2) = CC1_2(contour) / T(contour)
t4(contour::SchwingerKeldyshContourRounded2) = L2(contour) / T(contour)
t5(contour::SchwingerKeldyshContourRounded2) = CC2(contour) / T(contour)
t6(contour::SchwingerKeldyshContourRounded2) = L3(contour) / T(contour)

############### Schwinger-Keldysh contour Rounded 3 ###############
#
#   The Schwinger-Keldysh contour for one rounded corner, i.e.,
#   rounded corner at the $t_{max}$ point.
#
#
#   NR_RT_Points, NR_ET_Points, β, t_max, A_β and F_β is the same parameters
#   as in normal Schwinger-Keldysh contour
#
#   NR_Rounded1_Points: nr of points in the arc
#   r1: radius of the rounded corner
#
#
struct SchwingerKeldyshContourRounded3 <: AbstractSchwingerKeldyshContour
    β::Float32
    tmax::Float32
    A_β::Float64
    F_β::Float64
    r1::Float64
    r2::Float64
    r3::Float64
end

Base.string(io::IO, c::SchwingerKeldyshContourRounded3) = string("tmax: ", c.tmax, "\n", "β: ", c.β, "\n", "A_β: ", x.A_β, "\n", "F_β:", c.F_β, "\n",
                                                              "r1:", c.r1, "\n", "r2:", c.r2, "\n", "r3:", c.r3)

# Defining functions related to the calculation of the SchwingerKeldyshContour
#β1(contour::SchwingerKeldyshContourRounded2) = contour.A_β * contour.F_β * contour.β
#β2(contour::SchwingerKeldyshContourRounded2) = (1-contour.A_β) * contour.F_β * contour.β
#β3(contour::SchwingerKeldyshContourRounded2) = (1-contour.F_β) * contour.β

H1(contour::SchwingerKeldyshContourRounded3) = sqrt((contour.tmax-contour.r1-contour.r3)^2 + (β1(contour)+contour.r3)^2)
L1(contour::SchwingerKeldyshContourRounded3) = sqrt(H1(contour)^2 - (contour.r1+contour.r3)^2)
γ1(contour::SchwingerKeldyshContourRounded3) = pi/2 - asin((contour.r1+contour.r3)/H1(contour)) - acos((β1(contour)+contour.r3)/H1(contour))

H2(contour::SchwingerKeldyshContourRounded3) = sqrt((contour.tmax-contour.r1-contour.r2)^2 + β2(contour)^2)
L2(contour::SchwingerKeldyshContourRounded3) = sqrt(H2(contour)^2 - (contour.r1+contour.r2)^2)
γ2(contour::SchwingerKeldyshContourRounded3) = pi/2 - asin((contour.r1+contour.r2)/H2(contour)) - acos(β2(contour)/H2(contour))

Gy1(contour::SchwingerKeldyshContourRounded3) = L1(contour)*sin(γ1(contour))
Gx1(contour::SchwingerKeldyshContourRounded3) = L1(contour)*cos(γ1(contour))
Gy2(contour::SchwingerKeldyshContourRounded3) = L2(contour)*sin(γ2(contour))
Gx2(contour::SchwingerKeldyshContourRounded3) = L2(contour)*cos(γ2(contour))

η1(contour::SchwingerKeldyshContourRounded3) = pi/2 - γ1(contour)
η2(contour::SchwingerKeldyshContourRounded3) = pi/2 - γ2(contour)

CC1_1(contour::SchwingerKeldyshContourRounded3) = η1(contour)*contour.r1
CC1_2(contour::SchwingerKeldyshContourRounded3) = η2(contour)*contour.r1
CC2(contour::SchwingerKeldyshContourRounded3) = η2(contour)*contour.r2

L3(contour::SchwingerKeldyshContourRounded3) = β3(contour) - contour.r3*sin(η1(contour))
CC3(contour::SchwingerKeldyshContourRounded3) = η1(contour)*contour.r3

T(contour::SchwingerKeldyshContourRounded3)  = L1(contour) + CC1_1(contour) + CC1_2(contour) + L2(contour) + CC2(contour) + L3(contour) + CC3(contour)
t1(contour::SchwingerKeldyshContourRounded3) = L1(contour) / T(contour)
t2(contour::SchwingerKeldyshContourRounded3) = CC1_1(contour) / T(contour)
t3(contour::SchwingerKeldyshContourRounded3) = CC1_2(contour) / T(contour)
t4(contour::SchwingerKeldyshContourRounded3) = L2(contour) / T(contour)
t5(contour::SchwingerKeldyshContourRounded3) = CC2(contour) / T(contour)
t6(contour::SchwingerKeldyshContourRounded3) = L3(contour) / T(contour)
t7(contour::SchwingerKeldyshContourRounded3) = CC3(contour) / T(contour)

################################################
### Defining the related functions
################################################

"""
    Calculate the time point on the contour from the parametetric variable u in [0,1]
    for the Non-Equilibrium Schwinger-Keldysh contour
"""
function get_timepoint(u,C::NonEquilSchwingerKeldyshContour)
    if u <= t1(C)
        return complex(C.tmax*(u/t1(C)), 0)
    elseif t1(C) < u <= t1(C) + t2(C)
        return complex(C.tmax - C.tmax*((u-t1(C))/t2(C)))
    end
    return nothing
end

"""
    Calculate the time point on the contour from the parametetric variable u in [0,1]
    for the Non-Equilibrium tilted Schwinger-Keldysh contour
"""
function get_timepoint(u,C::NonEquilTiltedSchwingerKeldyshContour)
    if u <= t1(C)
        return complex(C.tmax*(u/t1(C)), -(C.tmax * C.θ)*(u/t1(C)))
    elseif t1(C) < u <= t1(C) + t2(C)
        return complex(C.tmax - C.tmax*((u-t1(C))/t2(C)),-(C.tmax * C.θ) -  (C.tmax * C.θ)* ((u-t1(C))/t2(C)))
    end
    return nothing
end

"""
    Calculate the time point on the contour from the parametetric variable u in [0,1]
    for the Schwinger-Keldysh contour
"""
function get_timepoint(u,C::SchwingerKeldyshContour)
    if u <= t1(C)
        return complex(C.tmax*(u/t1(C)), -β1(C)*(u/t1(C)))
    elseif t1(C) < u <= t1(C) + t2(C)
        return complex(C.tmax - C.tmax*((u-t1(C))/t2(C)), -β1(C) - β2(C)*((u-t1(C))/t2(C)))
    elseif u > t1(C) + t2(C)
        return complex(0, -β1(C) - β2(C) - L3(C)*((u-t1(C)-t2(C))/t3(C)))
    end
    return nothing
end

"""
    Calculate the time point on the contour from the parametetric variable u in [0,1]
    for the Schwinger-Keldysh contour
"""
function get_timepoint(u,C::SchwingerKeldyshContourRounded1)
    if (u <= t1(C))
        ui = u / t1(C)
        return complex(Gx1(C)*ui, - Gy1(C)*ui);
    elseif t1(C) < u <= t1(C) + t2(C)
        ui = (u - t1(C))/t2(C)
        return complex(C.tmax - C.r1 + C.r1*cos(η1(C) - η1(C)*ui), -β1(C) + C.r1*sin( η1(C) - η1(C)*ui));
    elseif t1(C) + t2(C) < u <= t1(C) + t2(C) + t3(C)
        ui = (u - t1(C) - t2(C))/t3(C)
        return complex(C.tmax - C.r1 + C.r1*cos(-η2(C)*ui), -β1(C)+C.r1*sin(-η2(C)*ui));
    elseif t1(C) + t2(C) + t3(C) < u <= t1(C) + t2(C) + t3(C) + t4(C)
        ui = (u - t1(C) - t2(C) - t3(C))/t4(C)
        return complex(Gx2(C) - Gx2(C)*ui, - (β1(C)+β2(C) - Gy2(C)) - Gy2(C)*ui);
    else
        ui = (u - t1(C) - t2(C) - t3(C) - t4(C))/t5(C)
        return complex(0, - β1(C)-β2(C) - ui*β3(C))
    end

    return nothing
end

"""
    Calculate the time point on the contour from the parametetric variable u in [0,1]
    for the Schwinger-Keldysh contour
"""
function get_timepoint(u,C::SchwingerKeldyshContourRounded2)
    if (u <= t1(C))
        ui = u / t1(C)
        return complex(Gx1(C)*ui, - Gy1(C)*ui);
    elseif t1(C) < u <= t1(C) + t2(C)
        ui = (u - t1(C))/t2(C)
        return complex(C.tmax - C.r1 + C.r1*cos(η1(C) - η1(C)*ui), -β1(C) + C.r1*sin( η1(C) - η1(C)*ui));
    elseif t1(C) + t2(C) < u <= t1(C) + t2(C) + t3(C)
        ui = (u - t1(C) - t2(C))/t3(C)
        return complex(C.tmax - C.r1 + C.r1*cos(-η2(C)*ui), -β1(C)+C.r1*sin(-η2(C)*ui));
    elseif t1(C) + t2(C) + t3(C) < u <= t1(C) + t2(C) + t3(C) + t4(C)
        ui = (u - t1(C) - t2(C) - t3(C))/t4(C)
        return complex(Gx2(C) + C.r2*(1-cos(pi/2-γ2(C))) - Gx2(C)*ui, - (β1(C)+β2(C) - Gy2(C) - C.r2*sin(pi/2 -γ2(C))) - Gy2(C)*ui);
    elseif t1(C) + t2(C) + t3(C) + t4(C) < u <= t1(C) + t2(C) + t3(C) + t4(C) + t5(C)
        ui = (u - t1(C) - t2(C) - t3(C) - t4(C))/t5(C)
        return complex(C.r2*(1 + cos(pi-η2(C)*(1-ui))), -β1(C)-β2(C) + C.r2*sin(pi-η2(C)*(1-ui)) );
    else
        ui = (u - t1(C) - t2(C) - t3(C) - t4(C) - t5(C))/t6(C)
        return complex(0, - β1(C)-β2(C) - ui*β3(C))
    end

    return nothing
end

"""
    Calculate the time point on the contour from the parametetric variable u in [0,1]
    for the Schwinger-Keldysh contour with three rounded corners
"""
function get_timepoint(u,C::SchwingerKeldyshContourRounded3)
    if (u <= t1(C))
        ui = u / t1(C)
        return complex(C.r3*(1-cos(η1(C))) + Gx1(C)*ui, C.r3*(1-sin(η1(C))) - Gy1(C)*ui);
    elseif t1(C) < u <= t1(C) + t2(C)
        ui = (u - t1(C))/t2(C)
        return complex(C.tmax - C.r1 + C.r1*cos(η1(C) - η1(C)*ui), -β1(C) + C.r1*sin( η1(C) - η1(C)*ui));
    elseif t1(C) + t2(C) < u <= t1(C) + t2(C) + t3(C)
        ui = (u - t1(C) - t2(C))/t3(C)
        return complex(C.tmax - C.r1 + C.r1*cos(-η2(C)*ui), -β1(C)+C.r1*sin(-η2(C)*ui));
    elseif t1(C) + t2(C) + t3(C) < u <= t1(C) + t2(C) + t3(C) + t4(C)
        ui = (u - t1(C) - t2(C) - t3(C))/t4(C)
        return complex(Gx2(C) + C.r2*(1-cos(pi/2-γ2(C))) - Gx2(C)*ui, - (β1(C)+β2(C) - Gy2(C) - C.r2*sin(pi/2 -γ2(C))) - Gy2(C)*ui);
    elseif t1(C) + t2(C) + t3(C) + t4(C) < u <= t1(C) + t2(C) + t3(C) + t4(C) + t5(C)
        ui = (u - t1(C) - t2(C) - t3(C) - t4(C))/t5(C)
        return complex(C.r2*(1 + cos(pi-η2(C)*(1-ui))), -β1(C)-β2(C) + C.r2*sin(pi-η2(C)*(1-ui)) );
    elseif t1(C) + t2(C) + t3(C) + t4(C) + t5(C) < u <= t1(C) + t2(C) + t3(C) + t4(C) + t5(C) + t6(C)
        ui = (u - t1(C) - t2(C) - t3(C) - t4(C) - t5(C))/t6(C)
        return complex(0, - β1(C)-β2(C) - ui*L3(C))
    else
        ui = (u - t1(C) - t2(C) - t3(C) - t4(C) - t5(C) - t6(C))/t7(C)
        return complex(C.r3*(1+cos(pi + η1(C)*ui)), -C.β + C.r3*( 1 + sin(pi + η1(C)*ui) ))
    end

    return nothing
end

##################### Calculating the derivative ############################

"""
    Caluclate the derivative of the timepoint in the contour
    given the parameter variable u in [0,1]

    The derivative is taken using the ForwardDiff.derivative

    This will work for all contours
"""
function get_derivative_timepoint(u,C::AbstractContour)
    return complex(
            ForwardDiff.derivative(x->real(get_timepoint(x,C)),u),
            ForwardDiff.derivative(x->imag(get_timepoint(x,C)),u)
        )
end

ext: jl


text:
"""
    PartialAggregation

Functor object for using partial aggregation in an L-shaped algorithm. Create by supplying a [`PartialAggregate`](@ref) object through `aggregate ` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Parameters
- `size::Int`: Number of cuts in each aggregate
...
"""
struct PartialAggregation{T <: AbstractFloat} <: AbstractAggregation
    size::Int
    num_subproblems::Int
    start_id::Int
    collection::CutCollection{T}

    function PartialAggregation(size::Integer, start_id::Integer, num_subproblems::Integer, ::Type{T}) where T <: AbstractFloat
        return new{T}(min(size, num_subproblems), num_subproblems, start_id, CutCollection(T, start_id))
    end
end
"""
    FullAggregation

Functor object for using complete aggregation in an L-shaped algorithm. Create by supplying an [`Aggregate`](@ref) object through `aggregate ` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
FullAggregation(start_id::Integer, num_subproblems::Integer, ::Type{T}) where T <: AbstractFloat = PartialAggregation(num_subproblems, start_id, num_subproblems, T)

function aggregate_cut!(lshaped::AbstractLShaped, aggregation::PartialAggregation, cut::HyperPlane)
    added = passthrough!(lshaped, aggregation, cut)
    add_to_collection!(aggregation.collection, cut, lshaped.x)
    if considered(aggregation.collection) == aggregation.size
        if collection_size(aggregation.collection) == aggregation.size
            aggregated_cut = aggregate(aggregation.collection)
            added |= add_cut!(lshaped, aggregated_cut)
        end
        renew!(aggregation.collection, aggregation.collection.id % num_thetas(lshaped) + 1)
    end
    return added
end

function passthrough!(lshaped::AbstractLShaped, aggregation::PartialAggregation, cut::HyperPlane{H}) where H <: HyperPlaneType
    return add_cut!(lshaped, HyperPlane(cut.δQ, cut.q, aggregation.collection.id, H))
end

function passthrough!(lshaped::AbstractLShaped, aggregation::PartialAggregation, cut::HyperPlane{OptimalityCut})
    return false
end

function aggregate_cut!(cutqueue::CutQueue, aggregation::PartialAggregation, ::MetaData, t::Integer, cut::HyperPlane, x::AbstractArray)
    passthrough!(cutqueue, aggregation, cut, t, x)
    add_to_collection!(aggregation.collection, cut, x)
    if considered(aggregation.collection) == aggregation.size
        if collection_size(aggregation.collection) == aggregation.size
            aggregated_cut = aggregate(aggregation.collection)
            put!(cutqueue, (t, aggregated_cut))
        end
        new_id = aggregation.start_id + (aggregation.collection.id - aggregation.start_id + 1) % ceil(Int, aggregation.num_subproblems/aggregation.size)
        renew!(aggregation.collection, new_id)
    end
    return nothing
end

function passthrough!(cutqueue::CutQueue, aggregation::PartialAggregation, cut::HyperPlane{H}, t::Integer, x::AbstractVector) where H <: HyperPlaneType
    put!(cutqueue, (t, HyperPlane(cut.δQ, cut.q, aggregation.collection.id, H)))
    return nothing
end

function passthrough!(::CutQueue, ::PartialAggregation, ::HyperPlane{OptimalityCut}, ::Integer, ::AbstractVector)
    return nothing
end

function num_thetas(num_subproblems::Integer, aggregation::PartialAggregation)
    return ceil(Int, num_subproblems / aggregation.size)
end

function num_thetas(num_subproblems::Integer, aggregation::PartialAggregation, sp::StochasticPrograms.ScenarioProblems)
    jobsize = ceil(Int, num_subproblems / nworkers())
    n = ceil(Int, jobsize/aggregation.size)
    remainder = num_subproblems - (nworkers() - 1) * jobsize
    aggregationem = ceil(Int, remainder/aggregation.size)
    return n * (nworkers() - 1) + aggregationem
end

function num_thetas(::Integer, aggregation::PartialAggregation, sp::DistributedScenarioProblems)
    return sum([ceil(Int, nscen/aggregation.size) for nscen in sp.scenario_distribution])
end

function flush!(lshaped::AbstractLShaped, aggregation::PartialAggregation)
    added = false
    if collection_size(aggregation.collection) > 0 && collection_size(aggregation.collection) == considered(aggregation.collection)
        aggregated_cut = aggregate(aggregation.collection)
        added |= add_cut!(lshaped, aggregated_cut)
    end
    renew!(aggregation.collection, aggregation.start_id)
    return added
end

function flush!(cutqueue::CutQueue, aggregation::PartialAggregation, ::MetaData, t::Integer, x::AbstractArray)
    if collection_size(aggregation.collection) > 0 && collection_size(aggregation.collection) == considered(aggregation.collection)
        aggregated_cut = aggregate(aggregation.collection)
        put!(cutqueue, (t, aggregated_cut))
    end
    renew!(aggregation.collection, aggregation.start_id)
    return nothing
end

# API
# ------------------------------------------------------------
"""
    PartialAggregate

Factory object for [`PartialAggregation`](@ref). Pass to `aggregate` in the `LShapedSolver` factory function.  See ?PartialAggregation for parameter descriptions.

"""
mutable struct PartialAggregate <: AbstractAggregator
    size::Int
    start_id::Int

    function PartialAggregate(size::Integer)
        return new(size, 1)
    end

    function PartialAggregate(size::Integer, start_id::Int)
        return new(size, start_id)
    end
end

function (aggregator::PartialAggregate)(num_subproblems::Integer, ::Type{T}) where T <: AbstractFloat
    aggregator.size == 1 && return NoAggregation()
    aggregator.size == num_subproblems && return FullAggregation(1, num_subproblems, T)
    return PartialAggregation(aggregator.size, aggregator.start_id, num_subproblems, T)
end

function remote_aggregator(aggregation::PartialAggregation, sp::ScenarioProblems, w::Integer)
    (nscen, extra) = divrem(num_subproblems(sp), nworkers())
    prev = map(2:(w - 1)) do p
        jobsize = nscen + (extra + 2 - p > 0)
        ceil(Int, jobsize / aggregation.size)
    end
    start_id = isempty(prev) ? 1 : sum(prev) + 1
    return PartialAggregate(aggregation.size, start_id)
end

function remote_aggregator(aggregation::PartialAggregation, sp::DistributedScenarioProblems, w::Integer)
    prev = map(2:(w - 1)) do p
        ceil(Int, sp.scenario_distribution[p-1] / aggregation.size)
    end
    start_id = isempty(prev) ? 1 : sum(prev) + 1
    nscen = sp.scenario_distribution[w-1]
    return PartialAggregate(aggregation.size, start_id)
end

function str(aggregator::PartialAggregate)
    return "partial cut aggregation of size $(aggregator.size)"
end

"""
    Aggregate

Factory object for [`FullAggregation`](@ref). Pass to `aggregate` in the `LShapedSolver` factory function.

"""
struct Aggregate <: AbstractAggregator end

function (::Aggregate)(num_subproblems::Integer, ::Type{T}) where T
    return FullAggregation(1, num_subproblems, T)
end

function str(::Aggregate)
    return "full cut aggregation"
end


text:
using Test
using Surrogates

@testset "Radials.jl" begin include("Radials.jl") end
@testset "Kriging.jl" begin include("Kriging.jl") end
@testset "Sampling" begin include("sampling.jl") end
@testset "Optimization" begin include("optimization.jl") end
@testset "LinearSurrogate" begin include("linearSurrogate.jl") end
@testset "Lobachesky" begin include("lobachesky.jl") end
@testset "RandomForestSurrogate" begin include("random_forest.jl") end
#@testset "SVMSurrogate" begin include("SVMSurrogate.jl") end
@testset "NeuralSurrogate" begin include("neuralSurrogate.jl") end
@testset "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end
@testset "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end
@testset "AD_Compatibility" begin include("AD_compatibility.jl") end


text:
using Propagators: Cell, propagator, propagate!, block, constant, switch

# function sqrt_heron(x)
#   g = x/x
#   while !(g^2 ≈ x)
#     g = (x/g+g)/2
#   end
#   return g
# end

function heron_step(x, g, h)
  block(x, g) do
    xdg = Cell()
    gpxdg = Cell()
    two = Cell()
    constant(2)(two)
    propagator(/)(x, g, xdg)
    propagator(+)(g, xdg, gpxdg)
    propagator(/)(gpxdg, two, h)
  end
end

function sqrt_network(x, answer)
  block(x) do
    one = Cell()
    constant(1)(one)
    sqrt_iter(x, one, answer)
  end
end

function sqrt_iter(x, g, answer)
  block(x, g) do
    done = Cell()
    notdone = Cell()
    xifnotdone = Cell()
    gifnotdone = Cell()
    newg = Cell()
    propagator((g, x) -> g^2 ≈ x)(g, x, done)
    propagator(!)(done, notdone)
    switch(done, g, answer)
    switch(notdone, x, xifnotdone)
    switch(notdone, g, gifnotdone)
    heron_step(xifnotdone, gifnotdone, newg)
    sqrt_iter(xifnotdone, newg, answer)
  end
end

# x = Cell()
# y = Cell()
# sqrt_network(x, y)
# x[] = 2
# propagate!()
# y[]


text:
# Searching
# ---------

function findidx(t::LeafNode{K, V, B}, key::AbstractSegment{K}) where {K, V, B}
    if isempty(t)
        return 0
    end
    # TODO: search keys
    return searchsortedfirst(t.entries, key, 1, Int(t.count),
                             Base.ord(basic_isless, identity, false))
end

function findidx(t::InternalNode{K, V, B}, key::AbstractSegment{K}) where {K, V, B}
    if isempty(t.keys)
        return 0
    end
    i = searchsortedfirst(t.keys, key)
    return min(i, length(t.keys))
end

function Base.haskey(t::LeafNode{K, V, B}, key::AbstractSegment{K}) where {K, V, B}
    i = findidx(t, key)
    return 1 <= i <= length(t) && first(t.keys[i]) == first(key) &&
           last(t.keys[i]) == last(key)
end

function Base.haskey(t::InternalNode{K, V, B}, key::AbstractSegment{K}) where {K, V, B}
    i = findidx(t, key)
    if i <= length(t) - 1 && !(key < t.keys[i])
        return haskey(t.children[i+1], key)
    else
        return haskey(t.children[i], key)
    end
end

function Base.haskey(t::SegmentBTree{K, V, B}, key0::Tuple{Any, Any}) where {K, V, B}
    key = Segment{K}(key0[1], key0[2])
    return haskey(t.root, key)
end

function Base.findfirst(t::LeafNode{K, V, B}, key::AbstractSegment{K}, f) where {K, V, B}
    i = findidx(t, key)
    while 1 <= i <= length(t) &&
          first(t.keys[i]) == first(key) &&
          last(t.keys[i]) == last(key)
        if f(t.entries[i], key)
            return t.entries[i]
        end
        i += 1
        if i > length(t) && t.right !== nothing
            t = notnothing(t.right)
            i = 1
        end
    end

    return nothing
end

function Base.findfirst(t::InternalNode{K, V, B}, key::AbstractSegment{K}, f) where {K, V, B}
    i = findidx(t, key)
    if i <= length(t) - 1 && !(key < t.keys[i])
        return findfirst(t.children[i+1], key, f)
    else
        return findfirst(t.children[i], key, f)
    end
end

true_cmp(a, b) = true

function Base.findfirst(t::SegmentBTree{K, V, B}, key::AbstractSegment{K},
                        f=true_cmp) where {K, V, B}
    return findfirst(t.root, key, f)
end

function Base.findfirst(t::SegmentBTree{K, V, B}, key0::Tuple{Any, Any},
                        f=true_cmp) where {K, V, B}
    return findfirst(t, Segment{K}(key0[1], key0[2]), f)
end


text:
"""
    Sauvola <: AbstractImageBinarizationAlgorithm
    Sauvola(; bias = 0.2, window_size=7)

    binarize([T,] img, f::Sauvola)
    binarize!([out,] img, f::Sauvola)

Applies Sauvola--Pietikäinen adaptive image binarization [1] under the
assumption that the input image is textual.

# Output

Return the binarized image as an `Array{Gray{T}}` of size `size(img)`. If
`T` is not specified, it is inferred from `out` and `img`.

# Details

The input image is binarized by varying the threshold across the image, using
a modified version of Niblack's algorithm [2]. Niblack's approach was to define
a threshold ``T`` for each pixel based on the mean ``m`` and standard deviation
``s`` of the intensities of neighboring pixels in a window around it, given by

```math
T(x,y) = m(x,y) + k \\cdot s(x,y),
```

where ``k`` is a user-defined parameter weighting the influence of the standard
deviation on the value of ``T``.

Niblack's algorithm is highly sensitive to variations in the gray values of
background pixels, which often exceed local thresholds and appear as artifacts
in the binarized image. Sauvola and Pietikäinen [1] introduce the dynamic range
``R`` of the standard deviation (i.e. its maximum possible value in the color
space), such that the threshold is given by

```math
T(x,y) = m(x,y) \\cdot \\left[ 1 + k \\cdot \\left( \\frac{s(x,y)}{R} - 1 \\right) \\right]
```

This adaptively amplifies the contribution made by the standard deviation to the
value of ``T``.

The Sauvola--Pietikäinen algorithm is implemented here using an optimization
proposed by Shafait, Keysers and Breuel [3], in which integral images are used
to calculate the values of ``m`` and ``s`` for each pixel in constant time.
Since each of these data structures can be computed in a single pass over the
source image, runtime is significantly improved.

# Arguments

## `img`

An image which is binarized according to a per-pixel adaptive
threshold into background (0) and foreground (1) pixel values.

## `window_size::Integer` (denoted by ``w`` in the publication)

The threshold for each pixel is a function of the distribution of the intensities
of all neighboring pixels in a square window around it. The side length of this
window is ``2w + 1``, with the target pixel in the center position.

## `bias::Real` (denoted by ``k`` in the publication)

A user-defined biasing parameter. This can take negative values, though values
in the range [0.2, 0.5] are typical. According to [1], this algorithm is not too
sensitive to the value of ``k```.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(img, Sauvola(window_size = 9, bias = 0.2))
```

# References

1. J. Sauvola and M. Pietikäinen (2000). "Adaptive document image binarization". *Pattern Recognition* 33 (2): 225-236. [doi:10.1016/S0031-3203(99)00055-2](https://doi.org/10.1016/S0031-3203(99)00055-2)
2. Wayne Niblack (1986). *An Introduction to Image Processing*. Prentice-Hall, Englewood Cliffs, NJ: 115-16.
3. Faisal Shafait, Daniel Keysers and Thomas M. Breuel (2008). "Efficient implementation of local adaptive thresholding techniques using integral images". Proc. SPIE 6815, Document Recognition and Retrieval XV, 681510 (28 January 2008). [doi:10.1117/12.767755](https://doi.org/10.1117/12.767755)
"""
struct Sauvola <: AbstractImageBinarizationAlgorithm
    window_size::Int
    bias::Float32

    function Sauvola(window_size::Integer, bias::Real)
        window_size < 0 && throw(ArgumentError("window_size should be non-negative."))
        new(window_size, bias)
    end
end

Sauvola(; window_size::Int = 7, bias::Real = 0.2) = Sauvola(window_size, bias)

function binarize(algorithm::Sauvola, img::AbstractArray{T,2}) where T <: Colorant
    binarize(algorithm, Gray.(img))
end

function (f::Sauvola)(out::GenericGrayImage, img::GenericGrayImage)
    window_size = f.window_size
    k = f.bias

    size(out) == size(img) || throw(ArgumentError("out and img should have the same shape, instead they are $(size(out)) and $(size(img))"))

    img_raw = channelview(img)
    I = integral_image(img_raw)
    I² = integral_image(img_raw.^2)
    R = 0.5 # dynamic range of standard deviation, in [1] it's set to 128 for 8-bit image

    function threshold(pixel::CartesianIndex{2})
        row₀, col₀, row₁, col₁ = get_window_bounds(img, pixel, window_size)
        m = μ_in_window(I, row₀, col₀, row₁, col₁)
        s = σ_in_window(I², m, row₀, col₀, row₁, col₁)
        return m * (1 + (k * ((s / R) - 1)))
    end

    @simd for pixel in CartesianIndices(img)
        out[pixel] = img[pixel] <= threshold(pixel) ? 0 : 1
    end

    return out
end

(f::Sauvola)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
using DataFrames, FixedEffectModels, Random, CategoricalArrays

# Very simple setup
N = 10000000
K = 100
id1 = rand(1:div(N, K), N)
id2 = rand(1:K, N)
x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)
df = DataFrame(id1 = id1, id2 = id2, x1 = x1, x2 = x2, y = y)
# first time
@time reg(df, @formula(y ~ x1 + x2))
# 14s
@time reg(df, @formula(y ~ x1 + x2))
# 0.582029 seconds (852 allocations: 535.311 MiB, 18.28% gc time)
@time reg(df, @formula(y ~ x1 + x2),  Vcov.cluster(:id2))
#   0.621690 seconds (693 allocations: 768.945 MiB, 7.69% gc time)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1)))
# 1.143941 seconds (245.39 k allocations: 942.937 MiB, 12.93% gc time, 14.99% compilation time)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1)), Vcov.cluster(:id1))
# 1.242207 seconds (245.73 k allocations: 1022.348 MiB, 9.48% gc time, 14.10% compilation time)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))
# 2.255812 seconds (351.74 k allocations: 1.076 GiB, 3.98% gc time, 12.93% compilation time)

# More complicated setup
N = 800000 # number of observations
M = 40000 # number of workers
O = 5000 # number of firms
id1 = rand(1:M, N)
id2 = [rand(max(1, div(x, 8)-10):min(O, div(x, 8)+10)) for x in id1]
x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)
df = DataFrame(id1 = id1, id2 = id2, x1 = x1, x2 = x2, y = y)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))
#   3.048292 seconds (422.51 k allocations: 114.317 MiB, 6.86% compilation time)

+# fixest
n = 10_000_000
nb_dum = [div(n,20), floor(Int, sqrt(n)), floor(Int, n^.33)]
N = nb_dum.^3
id1 = categorical(rand(1:nb_dum[1], n))
id2 = categorical(rand(1:nb_dum[2], n))
id3 = categorical(rand(1:nb_dum[3], n))
X1 = rand(n)
ln_y = 3 .* X1 .+ rand(n)
df = DataFrame(X1 = X1, ln_y = ln_y, id1 = id1, id2 = id2, id3 = id3)
@time reg(df, @formula(ln_y ~ X1 + fe(id1)), Vcov.cluster(:id1))
# 0.869512 seconds (234.23 k allocations: 828.818 MiB, 18.95% compilation time)
@time reg(df, @formula(ln_y ~ X1 + fe(id1) + fe(id2)), Vcov.cluster(:id1))
# 2.192262 seconds (300.08 k allocations: 985.534 MiB, 4.61% gc time, 9.42% compilation time)
@time reg(df, @formula(ln_y ~ X1 + fe(id1) + fe(id2) + fe(id3)), Vcov.cluster(:id1))
#  2.700051 seconds (406.80 k allocations: 1.117 GiB, 3.56% gc time, 10.41% compilation time)


text:
# Before running this, please make sure to activate and instantiate the
# environment with [this `Project.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-housekingcounty/Project.toml) and
# [this `Manifest.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/EX-housekingcounty/Manifest.toml).
# For instance, copy these files to a folder 'EX-housekingcounty', `cd` to it and
#
# ```julia
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# ## Getting started
#
# This tutorial is adapted from [the corresponding MLR3 tutorial](https://mlr3gallery.mlr-org.com/posts/2020-01-30-house-prices-in-king-county/).
#
# ### Loading and  preparing the data

using MLJ
using PrettyPrinting
import DataFrames: DataFrame, select!, Not, describe
import Statistics
using Dates
using PyPlot
using UrlDownload

df = DataFrame(urldownload("https://raw.githubusercontent.com/tlienart/DataScienceTutorialsData.jl/master/data/kc_housing.csv", true))
describe(df)

# We drop unrelated columns

select!(df, Not([:id, :date]))
schema(df)

# Afterwards, we convert the zip code to an unordered factor (`Multiclass`), we also create two binary features `isrenovated` and `has_basement` derived from `yr_renovated` and `sqft_basement`:

coerce!(df, :zipcode => Multiclass)
df.isrenovated  = @. !iszero(df.yr_renovated)
df.has_basement = @. !iszero(df.sqft_basement)
schema(df)

# These created variables should be treated as OrderedFactor,

coerce!(df, :isrenovated => OrderedFactor, :has_basement => OrderedFactor);

# The feature `waterfront` is currently encoded as a string, but it's really just a boolean:

unique(df.waterfront)

# So let's recode it

df.waterfront = (df.waterfront .!= "FALSE")
coerce!(df, :waterfront => OrderedFactor);

# For a number of the remaining features which are treated as `Count` there are few unique values in which case it might make more sense to recode them as OrderedFactor, this can be done with `autotype`:

coerce!(df, autotype(df, :few_to_finite))
schema(df)

# Let's also rescale the column `price` to be in 1000s of dollars:

df.price = df.price ./ 1000;

# For simplicity let's just drop a few additional columns that don't seem to matter much:

select!(df, Not([:yr_renovated, :sqft_basement, :zipcode]));

# ### Basic data visualisation
#
# Let's plot a basic histogram of the prices to get an idea for the distribution:

plt.figure(figsize=(8,6))
plt.hist(df.price, color = "blue", edgecolor = "white", bins=50,
         density=true, alpha=0.5)
plt.xlabel("Price", fontsize=14)
plt.ylabel("Frequency", fontsize=14)

# \figalt{Histogram of the prices}{hist_price.svg}

# Let's see if there's a difference between renovated and unrenovated flats:

plt.figure(figsize=(8,6))
plt.hist(df.price[df.isrenovated .== true], color="blue", density=true,
        edgecolor="white", bins=50, label="renovated", alpha=0.5)
plt.hist(df.price[df.isrenovated .== false], color="red", density=true,
        edgecolor="white", bins=50, label="unrenovated", alpha=0.5)
plt.xlabel("Price", fontsize=14)
plt.ylabel("Frequency", fontsize=14)
plt.legend(fontsize=12)

# \figalt{Histogram of the prices depending on renovation}{hist_price2.svg}
# We can observe that renovated flats seem to achieve higher sales values, and this might thus be a relevant feature.
#
#
# Likewise, this could be done to verify that `condition`, `waterfront` etc are important features.

# ## Fitting a first model

DTR = @load DecisionTreeRegressor pkg=DecisionTree

y, X = unpack(df, ==(:price))
train, test = partition(collect(eachindex(y)), 0.7, shuffle=true, rng=5)
tree = machine(DTR(), X, y)

fit!(tree, rows=train);

# Let's see how it does

rms(y[test], MLJ.predict(tree, rows=test))

# Let's try to do better.

# ### Random forest model

# We might be able to improve upon the RMSE using more powerful learners.

RFR = @load RandomForestRegressor pkg=ScikitLearn

# That model only accepts input in the form of `Count` and so we have to coerce all `Finite` types into `Count`:

coerce!(X, Finite => Count);

# Now we can fit

rf_mdl = RFR()
rf = machine(rf_mdl, X, y)
fit!(rf, rows=train)

rms(y[test], MLJ.predict(rf, rows=test))

# A bit better but it would be best to check this a bit more carefully:

cv3 = CV(; nfolds=3)
res = evaluate(rf_mdl, X, y, resampling=CV(shuffle=true),
               measure=rms, verbosity=0)

# ### GBM
#
# Let's try a different kind of model: Gradient Boosted Decision Trees from the package xgboost and we'll try to tune it too.

XGBR = @load XGBoostRegressor

# It expects a `Table(Continuous)` input so we need to coerce `X` again:

coerce!(X, Count => Continuous)

xgb  = XGBR()
xgbm = machine(xgb, X, y)
fit!(xgbm, rows=train)

rms(y[test], MLJ.predict(xgbm, rows=test))

# Let's try to tune it, first we define ranges for a number of useful parameters:

r1 = range(xgb, :max_depth, lower=3, upper=10)
r2 = range(xgb, :num_round, lower=1, upper=25);

# And now we tune, we use a very coarse resolution because we use so many ranges, `2^7` is already some 128 models...

tm = TunedModel(model=xgb, tuning=Grid(resolution=7),
                resampling=CV(rng=11), ranges=[r1,r2],
                measure=rms)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

rms(y[test], MLJ.predict(mtm, rows=test))

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
"""
    from_float_array(A)

Reinterpret a real array as an array of quaternions

The input array must have an initial dimension whose size is 4, because
successive indices in that dimension will be considered successive components
of the output quaternion.

Note that this returns a view of the original data [via
`reinterpret(reshape,...)`] only if the base type of the input array
`isbitstype`; otherwise, a new array of `Quaternion`s must be created, and the
memory copied.

See also [`to_float_array`](@ref).

"""
function from_float_array(A::AbstractArray{T}) where {T<:Real}
    isbitstype(T) ? from_float_array(Val(true), A) : from_float_array(Val(false), A)
end
function from_float_array(::Val{true}, A::AbstractArray{T}) where {T<:Real}
    reinterpret(reshape, Quaternion{T}, A)
end
function from_float_array(::Val{false}, A::AbstractArray{T}) where {T<:Real}
    @assert size(A, 1)==4 "First dimension of `A` must be 4, not $(size(A, 1))"
    Q = Array{Quaternion{T}}(undef, size(A)[2:end])
    @inbounds for (i, j) in zip(eachindex(Q), Base.Iterators.partition(eachindex(A), 4))
        @views Q[i] = Quaternion{T}(A[j])
    end
    Q
end

"""
    to_float_array(A)

View a quaternion array as an array of real numbers

The output array will have an extra initial dimension whose size is 4, because
successive indices in that dimension correspond to successive components of the
quaternion.

Note that this returns a view of the original data only if the base type of the
input array `isbitstype`; otherwise, a new array of that type must be created,
and the memory copied.

See also [`from_float_array`](@ref).
"""
function to_float_array(A::AbstractArray{<:AbstractQuaternion{T}}) where {T<:Real}
    isbitstype(T) ? to_float_array(Val(true), A) : to_float_array(Val(false), A)
end
function to_float_array(::Val{true}, A::AbstractArray{<:AbstractQuaternion{T}}) where {T<:Real}
    reinterpret(reshape, T, A)
end
function to_float_array(::Val{false}, A::AbstractArray{<:AbstractQuaternion{T}}) where {T<:Real}
    F = Array{T}(undef, (4, size(A)...))
    @inbounds for (i, j) in zip(eachindex(A), Base.Iterators.partition(eachindex(F), 4))
        @views F[j] .= A[i].components[:]
    end
    F
end
to_float_array(q::AbstractQuaternion) = collect(float(q).components)

"""
    to_euler_angles(R)

Open Pandora's Box.

If somebody is trying to make you use Euler angles, tell them no, and walk
away, and go and tell your mum.

You don't want to use Euler angles.  They are awful.  Stay away.  It's one
thing to convert from Euler angles to quaternions; at least you're moving in
the right direction.  But to go the other way?!  It's just not right.

Assumes the Euler angles correspond to the quaternion `R` via

    R = exp(α𝐤/2) * exp(β𝐣/2) * exp(γ𝐤/2)

where 𝐣 and 𝐤 rotate about the fixed ``y`` and ``z`` axes, respectively, so
this reprents an initial rotation (in the positive sense) through an angle
``γ`` about the axis ``z``, followed by a rotation through ``β`` about the axis
``y``, and a final rotation through ``α`` about the axis ``z``.  This is
equivalent to performing an initial rotation through ``α`` about the axis
``z``, followed by a rotation through ``β`` about the *rotated* axis ``y'``,
followed by a rotation through ``γ`` about the *twice-rotated* axis ``z''``.
The angles are naturally assumed to be in radians.

NOTE: Before opening an issue reporting something "wrong" with this function,
be sure to read all of [this
page](https://github.com/moble/quaternion/wiki/Euler-angles-are-horrible),
*especially* the very last section about opening issues or pull requests.

# Returns
- `αβγ::Vector{T}`

# Raises
- `AllHell` if you try to actually use Euler angles, when you could have been
  using quaternions like a sensible person.

# See Also
- [`from_euler_angles`](@ref): Create quaternion from Euler angles
- [`to_euler_phases`](@ref): Convert quaternion to Euler phases
- [`from_euler_phases`](@ref): Create quaternion from Euler phases
"""
function to_euler_angles(q::AbstractQuaternion)
    q = float(q)
    a0 = 2acos(√((q.w^2+q.z^2)/abs2(q)))
    a1 = atan(q.z, q.w)
    a2 = atan(-q.x, q.y)
    [a1+a2, a0, a1-a2]
end

"""
    from_euler_angles(α, β, γ)

Come over from the dark side.

Assumes the Euler angles correspond to the quaternion `R` via

    R = exp(α𝐤/2) * exp(β𝐣/2) * exp(γ𝐤/2)

where 𝐣 and 𝐤 rotate about the fixed ``y`` and ``z`` axes, respectively, so
this reprents an initial rotation (in the positive sense) through an angle
``γ`` about the axis ``z``, followed by a rotation through ``β`` about the axis
``y``, and a final rotation through ``α`` about the axis ``z``.  This is
equivalent to performing an initial rotation through ``α`` about the axis
``z``, followed by a rotation through ``β`` about the *rotated* axis ``y'``,
followed by a rotation through ``γ`` about the *twice-rotated* axis ``z''``.
The angles are naturally assumed to be in radians.

NOTE: Before opening an issue reporting something "wrong" with this function,
be sure to read all of [this
page](https://github.com/moble/quaternion/wiki/Euler-angles-are-horrible),
*especially* the very last section about opening issues or pull requests.

# See Also
- [`to_euler_angles`](@ref): Convert quaternion to Euler angles
- [`to_euler_phases`](@ref): Convert quaternion to Euler phases
- [`from_euler_phases`](@ref): Create quaternion from Euler phases
"""
function from_euler_angles(α, β, γ)
    Rotor(
        cos(β/2)*cos((α+γ)/2),
        -sin(β/2)*sin((α-γ)/2),
        sin(β/2)*cos((α-γ)/2),
        cos(β/2)*sin((α+γ)/2)
    )
end
from_euler_angles(αβγ) = from_euler_angles(αβγ...)

function to_euler_phases!(z::Array{Complex{T}}, R::AbstractQuaternion{T}) where {T}
    a = R[1]^2 + R[4]^2
    b = R[2]^2 + R[3]^2
    sqrta = √a
    sqrtb = √b
    if iszero(sqrta)
        zp = one(Complex{T})
    else
        zp = Complex{T}(R[1], R[4]) / sqrta  # exp[i(α+γ)/2]
    end
    if iszero(sqrtb)
        zm = one(Complex{T})
    else
        zm = Complex{T}(R[3], -R[2]) / sqrtb  # exp[i(α-γ)/2]
    end
    z[1] = zp * zm  # exp[iα]
    z[2] = Complex{T}((a - b), 2 * sqrta * sqrtb) / (a + b)  # exp[iβ]
    z[3] = zp * conj(zm)  # exp[iγ]
    z
end

"""
    to_euler_phases(q)
    to_euler_phases!(z, q)

Convert input quaternion to complex phases of Euler angles

Interpreting the input quaternion as a rotation (though its normalization
scales out), we can define the complex Euler phases from the Euler angles (α,
β, γ) as

    zₐ ≔ exp(i*α)
    zᵦ ≔ exp(i*β)
    zᵧ ≔ exp(i*γ)

These are more useful geometric quantites than the angles themselves — being
involved in computing spherical harmonics and Wigner's 𝔇 matrices — and can be
computed from the components of the corresponding quaternion algebraically
(without the use of transcendental functions).

# Returns
- `z::Vector{Complex{T}}`: complex phases (zₐ, zᵦ, zᵧ) in that order.

# See Also
- [`from_euler_phases`](@ref): Create quaternion from Euler phases
- [`to_euler_angles`](@ref): Convert quaternion to Euler angles
- [`from_euler_angles`](@ref): Create quaternion from Euler angles

"""
function to_euler_phases(R::AbstractQuaternion{T}) where {T}
    z = Array{Complex{T}}(undef, 3)
    to_euler_phases!(z, R)
    z
end

"""
    from_euler_phases(zₐ, zᵦ, zᵧ)
    from_euler_phases(z)

Return the quaternion corresponding to these Euler phases.

Interpreting the input quaternion as a rotation (though its normalization
scales out), we can define the complex Euler phases from the Euler angles (α,
β, γ) as

    zₐ ≔ exp(i*α)
    zᵦ ≔ exp(i*β)
    zᵧ ≔ exp(i*γ)

These are more useful geometric quantites than the angles themselves — being
involved in computing spherical harmonics and Wigner's 𝔇 matrices — and can be
computed from the components of the corresponding quaternion algebraically
(without the use of transcendental functions).

# Parameters

- `z::Vector{Complex{T}}`: complex vector of length 3, representing the complex
  phases (zₐ, zᵦ, zᵧ) in that order.

# Returns
- `R::Quaternion{T}`

# See Also
- [`to_euler_phases`](@ref): Convert quaternion to Euler phases
- [`to_euler_angles`](@ref): Convert quaternion to Euler angles
- [`from_euler_angles`](@ref): Create quaternion from Euler angles

"""
function from_euler_phases(zₐ, zᵦ, zᵧ)
    zb = √(zᵦ)  # exp[iβ/2]
    zp = √(zₐ * zᵧ)  # exp[i(α+γ)/2]
    zm = √(zₐ * conj(zᵧ))  # exp[i(α-γ)/2]
    if abs(zₐ - zp * zm) > abs(zₐ + zp * zm)
        zp *= -1
    end
    Rotor(zb.re * zp.re, -zb.im * zm.im, zb.im * zm.re, zb.re * zp.im)
end
from_euler_phases(z) = from_euler_phases(z...)

"""
    to_spherical_coordinates(q)

Return the spherical coordinates corresponding to this quaternion.

We can treat the quaternion as a transformation taking the ``z`` axis to some
direction ``n̂``.  This direction can be described in terms of spherical
coordinates (θ, ϕ).  Here, we use the standard commonly used in physics: θ
represents the "polar angle" between the ``z`` axis and the direction ``n̂``,
while ϕ represents the "azimuthal angle" between the ``x`` axis and the
projection of ``n̂`` into the ``x``-``y`` plane.  Both angles are given in
radians.

"""
function to_spherical_coordinates(q::Q) where {Q<:AbstractQuaternion}
    q = float(q)
    a0 = 2acos(√((q.w^2+q.z^2)/abs2(q)))
    a1 = atan(q.z, q.w)
    a2 = atan(-q.x, q.y)
    [a0, a1+a2]
end

"""
    from_spherical_coordinates(θ, ϕ)

Return a rotor corresponding to these spherical coordinates.

Considering (θ, ϕ) as a point ``n̂`` on the sphere, we can also construct a
quaternion that rotates the ``z`` axis onto that point.  Here, we use the
standard commonly used in physics: θ represents the "polar angle" between the
``z`` axis and the direction ``n̂``, while ϕ represents the "azimuthal angle"
between the ``x`` axis and the projection of ``n̂`` into the ``x``-``y`` plane.
Both angles must be given in radians.

"""
function from_spherical_coordinates(θ, ϕ)
    sϕ, cϕ = sincos(ϕ/2)
    sθ, cθ = sincos(θ/2)
    Rotor(cθ*cϕ, -sθ*sϕ, sθ*cϕ, cθ*sϕ)
end
from_spherical_coordinates(θϕ) = from_spherical_coordinates(θϕ...)

"""
    from_rotation_matrix(ℛ)

Convert 3x3 rotation matrix to quaternion.

Assuming the 3x3 matrix `ℛ` rotates a vector `v` according to

    v' = ℛ * v,

we can also express this rotation in terms of a quaternion `R` such that

    v' = R * v * R⁻¹.

This function returns that quaternion, using Bar-Itzhack's algorithm to allow
for non-orthogonal matrices.  [J. Guidance, Vol. 23, No. 6,
p. 1085](http://dx.doi.org/10.2514/2.4654)

"""
function from_rotation_matrix(ℛ)
    K = Array{eltype(ℛ), 2}(undef, (4, 4))
    K[1, 1] = (ℛ[1, 1] - ℛ[2, 2] - ℛ[3, 3])/3
    K[1, 2] = (ℛ[2, 1] + ℛ[1, 2])/3
    K[1, 3] = (ℛ[3, 1] + ℛ[1, 3])/3
    K[1, 4] = (ℛ[2, 3] - ℛ[3, 2])/3
    K[2, 2] = (ℛ[2, 2] - ℛ[1, 1] - ℛ[3, 3])/3
    K[2, 3] = (ℛ[3, 2] + ℛ[2, 3])/3
    K[2, 4] = (ℛ[3, 1] - ℛ[1, 3])/3
    K[3, 3] = (ℛ[3, 3] - ℛ[1, 1] - ℛ[2, 2])/3
    K[3, 4] = (ℛ[1, 2] - ℛ[2, 1])/3
    K[4, 4] = (ℛ[1, 1] + ℛ[2, 2] + ℛ[3, 3])/3
    H = Symmetric(K)

    # compute the *dominant* (largest eigenvalue) eigenvector
    eigenvec = eigen(transpose(H), 4:4).vectors[:, 1]

    # convert it into a quaternion
    Rotor(eigenvec[4], -eigenvec[1], -eigenvec[2], -eigenvec[3])
end

"""
    to_rotation_matrix(q)

Convert quaternion to 3x3 rotation matrix.

Assuming the quaternion `R` rotates a vector `v` according to

    v' = R * v * R⁻¹,

we can also express this rotation in terms of a 3x3 matrix `ℛ` such that

    v' = ℛ * v.

This function returns that matrix.

"""
function to_rotation_matrix(q::Q) where {Q<:AbstractQuaternion}
    m = Array{eltype(q)}(undef, 3, 3)
    n = inv(abs2(q))
    m[1, 1] = 1 - 2*(q.y^2 + q.z^2) * n
    m[1, 2] = 2*(q.x*q.y - q.z*q.w) * n
    m[1, 3] = 2*(q.x*q.z + q.y*q.w) * n
    m[2, 1] = 2*(q.x*q.y + q.z*q.w) * n
    m[2, 2] = 1 - 2*(q.x^2 + q.z^2) * n
    m[2, 3] = 2*(q.y*q.z - q.x*q.w) * n
    m[3, 1] = 2*(q.x*q.z - q.y*q.w) * n
    m[3, 2] = 2*(q.y*q.z + q.x*q.w) * n
    m[3, 3] = 1 - 2*(q.x^2 + q.y^2) * n
    m
end


text:
# ---
# title: 1137. N-th Tribonacci Number
# id: problem1137
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Recursion
# link: <https://leetcode.com/problems/n-th-tribonacci-number/description/>
# hidden: true
# ---
#
# The Tribonacci sequence Tn is defined as follows:
#
# T0 = 0, T1 = 1, T2 = 1, and Tn+3 = Tn \+ Tn+1 \+ Tn+2 for n >= 0.
#
# Given `n`, return the value of Tn.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 4
#     Output: 4
#     Explanation:
#     T_3 = 0 + 1 + 1 = 2
#     T_4 = 1 + 1 + 2 = 4
#
#
# **Example 2:**
#
#
#
#     Input: n = 25
#     Output: 1389537
#
#
#
#
# **Constraints:**
#
#   * `0 <= n <= 37`
#   * The answer is guaranteed to fit within a 32-bit integer, ie. `answer <= 2^31 - 1`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
_is_plotable_history(::UnivalueHistory) = false
_is_plotable_history{I,V<:Real}(::QHistory{I,V}) = true
_is_plotable_history{I,V<:Real}(::History{I,V}) = true

_filter_plotable_histories(h::MVHistory) =
    filter((k,v) -> _is_plotable_history(v), h.storage)

@recipe function plot(h::Union{History,QHistory})
    markershape --> :ellipse
    title       --> "Value History"
    get(h)
end

@recipe function plot(h::MVHistory)
    filtered = _filter_plotable_histories(h)
    k_vec = [k for (k, v) in filtered]
    v_vec = [v for (k, v) in filtered]
    if length(v_vec) > 0
        markershape --> :ellipse
        label       --> map(string, k_vec)'
        if get(d, :layout, nothing) != nothing
            title  --> d[:label]
            legend --> false
        else
            title  --> "Multivalue History"
        end
        get_vec = map(get, v_vec)
        [x for (x, y) in get_vec], [y for (x, y) in get_vec]
    else
        throw(ArgumentError("Can't plot an empty history, nor a history with strange types"))
    end
end

@recipe function plot{T<:ValueHistories.UnivalueHistory}(hs::AbstractVector{T})
    for h in hs
        @series begin
            h
        end
    end
end


text:
struct Comparator{T₁ <: KernelDensity.UnivariateKDE,
                  T₂ <: KernelDensity.UnivariateKDE}
    pdf₁::KernelDensity.InterpKDE{T₁}
    pdf₂::KernelDensity.InterpKDE{T₂}

    function Comparator(samples::T₁,
                        targets::T₂,
                        weights::T₃) where {T₁ <: AbstractMatrix,
                                            T₂ <: AbstractVector,
                                            T₃ <: AbstractVector}
        X = samples
        t = targets
        w = weights
        h = vec(w' * X)

        # Separate the random projections by label.
        # TODO Perhaps make more robust by deriving (binary) values from t,
        # rather than assuming they'll always be 1 and not-1.
        h₀ = h[t .!= 1]
        h₁ = h[t .== 1]

        # Set lower and upper bounds on the bandwidth of the KDEs, with some
        # allowance for interpolating values that will fall outside the range of
        # the training data.
        # TODO The constant term here should be a parameter to the model.
        bounds = (min(minimum(h₀), minimum(h₁)) - 1,
                  max(maximum(h₀), maximum(h₁)) + 1)

        # Calculate kernel density estimates for each per-label subset.
        # TODO The npoints value here should be a parameter to the model.
        kde₀ = KernelDensity.kde_lscv(h₀, npoints = 100, boundary = bounds)
        kde₁ = KernelDensity.kde_lscv(h₁, npoints = 100, boundary = bounds)

        # Interpret each KDE as a probability density function.
        pdf₀ = KernelDensity.InterpKDE(kde₀)
        pdf₁ = KernelDensity.InterpKDE(kde₁)

        param(::KernelDensity.InterpKDE{T}) where {T} = T
        new{param(pdf₀),param(pdf₁)}(pdf₀, pdf₁)
    end
end

(c::Comparator)(x::T) where {T <: Number} = (Distributions.pdf(c.pdf₁, x) -
                                             Distributions.pdf(c.pdf₂, x))


text:
module BenchAvgFilter2D

import Random
using BenchmarkTools
using TapirBenchmarks

Random.seed!(1234)

const SUITE = BenchmarkGroup()

for n in [2^9, 2^11]
    SUITE[:n=>n] = s0 = BenchmarkGroup()
    ys0, xs = avgfilter2d_setup(n)
    s0[:impl=>:seq] = @benchmarkable avgfilter2d_seq!($(copy(ys0)), $xs)
    s0[:impl=>:threads] = @benchmarkable avgfilter2d_threads!($(copy(ys0)), $xs)
    s0[:impl=>:tapir_dac] = @benchmarkable avgfilter2d_tapir_dac!($(copy(ys0)), $xs)
    # s0[:impl=>:tapir_seq] = @benchmarkable avgfilter2d_tapir_seq!($(copy(ys0)), $xs)
end

end
BenchAvgFilter2D.SUITE


text:
using Xpress
Solver = Xpress

using BenchmarkTools

BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

const nvars = 1000
const ncons = 1000
const ncoefs = 100

const inds = collect(1:div(nvars,ncoefs):nvars)
const coefs = collect(1.0:Float64(ncoefs))
@assert length(inds) == length(coefs)

const At = sparse(repmat(inds, ncons), repmat(collect(1:ncoefs),1,ncons)'[:], repmat(coefs, ncons), nvars, ncons)
const rhs = fill(24.1,ncons)

# env = Solver.Env()
# m1 = Solver.Model(env,"model")
m1 = Solver.Model()
Solver.add_cvars!(m1, zeros(nvars))
function add100cons1(m)
    for i in 1:ncons
        Solver.add_constr!(m, inds, coefs, '<', 24.1)
    end
end

m2 = Solver.Model()
Solver.add_cvars!(m2, zeros(nvars))
function add100cons2(m)

    add_constrs_t!(m, At, '<', rhs)
    nothing
end

# @benchmark add100cons1(m1)
# @benchmark add100cons2(m2)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GSWC"
version = v"3.0.5"

# Collection of sources required to build GSWCBuilder
sources = [
    "https://github.com/TEOS-10/GSW-C/archive/master.zip" =>
    "4d77c322b9a8ac107b0fee3fe36267e43384ae6cb2532ab8cd883e80d086aa5b",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd GSW-C-master/
$CC $CFLAGS -fPIC -c -O3 -Wall  gsw_oceanographic_toolbox.c gsw_saar.c
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ] || [ $target = "x86_64-apple-darwin14" ]; then   LD=$CC; elif [ $target = "x86_64-unknown-freebsd11.1" ]; then   LD=ld; fi
$LD $LDFLAGS -fPIC -shared -o libgswteos-10.$dlext gsw_oceanographic_toolbox.o gsw_saar.o -lm
mkdir $prefix/lib/
cp libgswteos-10.$dlext $prefix/lib/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

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

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: kouketsu/GSWCBuilder

max_issues_repo_name: kouketsu/GSWCBuilder

max_forks_repo_name: kouketsu/GSWCBuilder


text:
module test_testset_verbose

using Test

if VERSION >= v"1.6.0-DEV.1437"
    @testset "Verbose 1" verbose = true begin
        @testset "Verbose 2" begin
            @test true
        end
    end
end

end # module test_testset_verbose


text:
module RecursiveFactorization

include("./lu.jl")

end # module


text:
function timejulia(i)
    function steps(env, j)
        for _ in 1:j
            s, r, d, info = step!(env, rand(1:3))
            d && (s = reset!(env))
        end
    end

    env = Environment("../envs/Basic", logFile=pwd()*"/../envs/logs/BasicTestJl.log")
    steps(env, 10)
    t = @elapsed steps(env, i)
    close!(env)
    return t
end

function timepython(i)
    py"""
    from mlagents_envs.environment import UnityEnvironment
    from gym_unity.envs import UnityToGymWrapper
    import time
    import random
    import os

    def pysteps(i):
        def steps(env, j):
            start = time.time()
            for _ in range(0, j):
                s, r, d, info = env.step(random.randint(0, 2))
                if d == True:
                    s = env.reset()
            return time.time() - start

        unity_env = UnityEnvironment(
            "../envs/Basic", no_graphics=True, args=["-logFile", os.getcwd()+"/../envs/logs/BasicTestPy.log"]
        )
        env = UnityToGymWrapper(unity_env)
        t = steps(env, i)
        env.close()
        return t
    """
    return py"pysteps"(i)
end

function speedtest(steps)
    @testset "Speed test" begin
        t1 = timejulia(steps)
        t2 = timepython(steps)
        @test t1 < t2 * 1.05
    end;
end


text:
module B_structs

using JuMP
using O_structs

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- structures of data : master ------------------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

export D__type                               # structure of problem data
export R__type                               # structure of relax master problem
export T1_type                               # structure of temporary data (updated each iteration)    -> valid only for algorithm 1 (Stand_Bend)
export T2_type                               # structure of temporary data (updated each iteration)    -> valid only for algorithm 2 (Adapt_Bend)
export T3_type                               # structure of temporary data (updated each iteration)    -> valid only for algorithm 3 (Zaker_Bend)
export H__type                               # structure of data stored during the solution
export M__type                               # structure of exact solutions stored during the solution
export B1_type                               # structure of Benders master problem                     -> valid only for algorithm 1 (Stand_Bend)
export B2_type                               # structure of Benders master problem                     -> valid only for algorithm 2 (Adapt_Bend)
export B3_type                               # structure of Benders master problem                     -> valid only for algorithm 3 (Zaker_Bend)

mutable struct D__type
    ms::ms_type                              # data of master problem sets
    mp::mp_type                              # data of master problem parameters
    ps::ps_type                              # data of subproblem sets
    pp::pp_type                              # data of subproblem parameters
   unc::u_type                               # data of uncertain parameters
  case::Int64                                # case to solve (0, 1, 2, or 3)
  algm::Int64                                # algorithm used (0, 1, 2, or 3)
     q::Float64                              # tolerance damping parameter                            -> valid only for algorithm 3 (Zaker_Bend)
     w::Int64                                # number of subproblems solved each iteration            -> valid only for algorithm 2 (Adapt_Bend)
  kidx::Array{BitArray{1},1}                 # array of binary association with subset ℰₖ             -> valid only for algorithm 2 (Adapt_Bend)
  iidx::Array{Array{Int64,1},1}              # array of subsets {ℰₖ, 𝑘=1,..,w}                        -> valid only for algorithm 2 (Adapt_Bend)
     J::Int64                                # maximum number of iteration
     δ::Float64                              # convergence tolerance
end

mutable struct R__type
     m::JuMP.Model                           # relaxed master problem JuMP model
  vars::Array{JuMP.VariableRef,2}            # reference to variablems {(βᵢ,xᵢ), ∀i}
  cons::Array{JuMP.ConstraintRef,2}          # reference to Benders cuts
end

mutable struct T1_type
     x::Array{Float64,2}                     # decisions {xᵢ, ∀i} set by the master problem @ iter k
     c::Array{Float64,2}                     # parameters {cᵢ, ∀i}
     β::Array{Float64,1}                     # variables value {βᵢ, ∀i} of the master problem @ iter k
     θ::Array{Float64,1}                     # optimal objectives {θᵢ, ∀i} of subproblems @ iter k
     λ::Array{Float64,2}                     # subgradient {λᵢ, ∀i} wrt xᵢ of subproblems @ iter k
     Δ::Float64                              # optimality gap @ iter k
end

mutable struct T2_type
     x::Array{Float64,2}                     # decisions {xᵢ, ∀i} set by the master problem @ iter k
     c::Array{Float64,2}                     # parameters {cᵢ, ∀i}
     β::Array{Float64,1}                     # variables value {βᵢ, ∀i} of the master problem @ iter k
    θl::Array{Float64,1}                     # valid lower bounds on objectives {θᵢ, ∀i} of subproblems @ iter k
    λl::Array{Float64,2}                     # valid subgradients {λᵢ, ∀i} wrt xᵢ of subproblems @ iter k
    θu::Array{Float64,1}                     # valid upper bounds on objectives {θᵢ, ∀i} of subproblems @ iter k
    ϕu::Array{Float64,2}                     # valid subgradients {ϕᵢ, ∀i} wrt cᵢ of subproblems @ iter k
    Ie::Array{Int64,1}                       # index of w subproblems solved exactly @ iter k
     Δ::Float64                              # optimality gap @ iter k
end

mutable struct T3_type
     x::Array{Float64,2}                     # decisions {xᵢ, ∀i} set by the master problem @ iter k
     c::Array{Float64,2}                     # parameters {cᵢ, ∀i}
     β::Array{Float64,1}                     # variables value {βᵢ, ∀i} of the master problem @ iter k
    θp::Array{Float64,1}                     # valid upper bounds on objectives {θᵢ, ∀i} of subproblems @ iter k
    θd::Array{Float64,1}                     # valid lower bounds on objectives {θᵢ, ∀i} of subproblems @ iter k
     λ::Array{Float64,2}                     # valid subgradients {λᵢ, ∀i} wrt xᵢ of subproblems @ iter k
     Δ::Float64                              # optimality gap @ iter k
end

mutable struct H__type
     L::Array{Float64,1}                     # lower bound values (up to iter k)
     U::Array{Float64,1}                     # upper bound values (up to iter k)
     t::Array{Float64,1}                     # time spend to perform each iteration (up to iter k)
     T::Array{Float64,2}                     # time spend to perform each iteration, split by steps (up to iter k)
     k::Int64                                # current iteration number
end

mutable struct M__type
     n::Int64                                # number of exact solution stored
     θ::Array{Float64,1}                     # optimal objective of exact solutions stored
     λ::Array{Float64,2}                     # subgradient wrt x of exact solutions stored
     ϕ::Array{Float64,2}                     # subgradient wrt c of exact solutions stored
     x::Array{Float64,2}                     # values of x of exact solutions stored
     c::Array{Float64,2}                     # values of c of exact solutions stored
end

mutable struct B1_type
   rmp::R__type                              # relaxed master problem
  temp::T1_type                              # temporary data (updated each iteration)
  hist::H__type                              # data stored during the solution
  data::D__type                              # problem data
end

mutable struct B2_type
   rmp::R__type                              # relaxed master problem
  temp::T2_type                              # temporary data (updated each iteration)
  hist::H__type                              # data stored during the solution
    m0::M__type                              # exact solutions stored
  data::D__type                              # problem data
 end

mutable struct B3_type
   rmp::R__type                              # relaxed master problem
  temp::T3_type                              # temporary data (updated each iteration)
  hist::H__type                              # data stored during the solution
  data::D__type                              # problem data
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- structures of data : subproblem --------------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

export E__type                               # structure of exact subproblem
export O__type                               # structure of oracles
export Q__type                               # structure of extended exact solutions stored during the solution
export d__type                               # structure of subproblem data
export t1_type                               # structure of temporary data (updated each iteration)    -> valid only for algorithm 1 (Stand_Bend)
export t2_type                               # structure of temporary data (updated each iteration)    -> valid only for algorithm 2 (Adapt_Bend)
export t3_type                               # structure of temporary data (updated each iteration)    -> valid only for algorithm 3 (Zaker_Bend)
export S1_type                               # structure of Benders subproblem                         -> valid only for algorithm 1 (Stand_Bend)
export S2_type                               # structure of Benders subproblem                         -> valid only for algorithm 2 (Adapt_Bend)
export S3_type                               # structure of Benders subproblem                         -> valid only for algorithm 3 (Zaker_Bend)

mutable struct E__type
      m::JuMP.Model                          # exact subproblem JuMP model
   vars::Array{JuMP.VariableRef,1}           # reference to variables (c₀,ϕ)
   objf::GenericAffExpr{Float64,VariableRef} # reference to objective function
end

mutable struct O__type
      m::JuMP.Model                          # oracle JuMP model
   cons::Array{JuMP.ConstraintRef,1}         # reference to oracle new constraints (new exact solutions)
   vars::Array{JuMP.VariableRef,1}           # reference to variables (ϕ,γ)
   objf::GenericAffExpr{Float64,VariableRef} # reference to objective function
   cexp::GenericAffExpr{Float64,VariableRef} # reference to objective function
   help::Array{Float64,1}                    # rhs values of old constraints (old exact solutions)
      n::Int64                               # number of exact solutions added to the oracle
end

mutable struct Q__type
     ϕ::Array{Float64,1}                     # operational cost dependent on uncertain cost cᵢ
    c0::Float64                              # operational cost independent of uncertain costs c
    yG::Array{Float64,3}                     # generation level of conventional generation
    yI::Array{Float64,3}                     # charging power of storage generation
    yO::Array{Float64,3}                     # discharging power of storage generation
    yL::Array{Float64,3}                     # energy level of storage generation
    yS::Array{Float64,2}                     # load shedding
    x0::Array{Float64,1}                     # values of rhs parameters in the subproblem (eg, capacity)
 end

mutable struct d__type
     ps::ps_type                             # data of subproblem sets
     pp::pp_type                             # data of subproblem parameters
    unc::u_type                              # data of uncertain parameters
   algm::Int64                               # algorithm used (0, 1, 2, or 3)
      J::Int64                               # maximum number of exact solution to store
end

mutable struct t1_type
      x::Array{Float64,2}                    # decisions {xᵢ, ∀i} set by the master problem @ iter k
      c::Array{Float64,2}                    # parameters {cᵢ, ∀i}
      θ::Array{Float64,1}                    # optimal objectives {θᵢ, ∀i} of subproblems @ iter k
      λ::Array{Float64,2}                    # subgradient {λᵢ, ∀i} wrt xᵢ of subproblems @ iter k
      i::Int64                               # index of subproblem to solve
end

mutable struct t2_type
     x::Array{Float64,2}                     # decisions {xᵢ, ∀i} set by the master problem @ iter k
     c::Array{Float64,2}                     # parameters {cᵢ, ∀i}
    θl::Array{Float64,1}                     # valid lower bounds on objectives {θᵢ, ∀i} of subproblems @ iter k
    θu::Array{Float64,1}                     # valid upper bounds on objectives {θᵢ, ∀i} of subproblems @ iter k
    λl::Array{Float64,2}                     # valid subgradients {λᵢ, ∀i} wrt xᵢ of subproblems @ iter k
    ϕu::Array{Float64,2}                     # valid subgradients {ϕᵢ, ∀i} wrt cᵢ of subproblems @ iter k
     i::Int64                                # index of subproblem to solve
end

mutable struct t3_type
     x::Array{Float64,2}                     # decisions {xᵢ, ∀i} set by the master problem @ iter k
     c::Array{Float64,2}                     # parameters {cᵢ, ∀i}
    θp::Array{Float64,1}                     # valid upper bounds on objectives {θᵢ, ∀i} of subproblems @ iter k
    θd::Array{Float64,1}                     # valid lower bounds on objectives {θᵢ, ∀i} of subproblems @ iter k
     λ::Array{Float64,2}                     # valid subgradients {λᵢ, ∀i} wrt xᵢ of subproblems @ iter k
     i::Int64                                # optimality gap @ iter k
end

mutable struct S1_type
     ex::E__type                             # exact subproblem
   data::d__type                             # subproblem data
   temp::t1_type                             # temporary data (updated each iteration)
end

mutable struct S2_type
     ex::E__type                             # exact subproblem
    olb::O__type                             # oracle for lower bound
    oub::O__type                             # oracle for upper bound
      m::M__type                             # exact solutions stored
      q::Array{Q__type,1}                    # extended exact solutions stored
   data::d__type                             # subproblem data
   temp::t2_type                             # temporary data (updated each iteration)
end

mutable struct S3_type
     ex::E__type                             # exact subproblem
   data::d__type                             # subproblem data
   temp::t3_type                             # temporary data (updated each iteration)
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

end


text:
"""
Module grf3D_Threads

Provides CPU (Threads) functions for 3D Gaussian random fields generation with anisotropic exponential and isotropic Gaussian covariance functions.

# Usage
    using ParallelRandomFields.grf3D_Threads

# Functions
- [`grf3D_expon!()`](@ref)
- [`grf3D_gauss!()`](@ref)
- [`generate_grf3D()`](@ref)

To see a description of a function type `?<functionname>`.
"""
module grf3D_Threads

    export grf3D_expon!, grf3D_gauss!, generate_grf3D

    @doc raw"""

        grf3D_expon!(Yf, sf, cl, nh, nx, ny, nz, dx, dy, dz; me=0, co1=0, co2=0, co3=0, do_reset=true)

    Compute the Gaussian random field `Yf` with anisotropic exponential covariance function and returns the wall-time `wtime_it`.

    # Arguments
    - `Yf::Data.Array`: the Gaussian random field array.
    - `sf::Data.Number`: the standard deviation.
    - `cl::Any`: the correlation length tuple `(clx, cly, clz)` in `x`, `y` and `z` dimension.
    - `nh::Int`: the number of harmonics.
    - `nx::Int`: the number of grid points in `x` dimension.
    - `ny::Int`: the number of grid points in `y` dimension.
    - `nz::Int`: the number of grid points in `z` dimension.
    - `dx::Data.Number`: the numerical grid size in `x` dimension.
    - `dy::Data.Number`: the numerical grid size in `y` dimension.
    - `dz::Data.Number`: the numerical grid size in `z` dimension.
    - `me=0::Int`: the optional kwarg mpi process id.
    - `co1=0::Int`: the optional kwarg mpi coord x.
    - `co2=0::Int`: the optional kwarg mpi coord y.
    - `co3=0::Int`: the optional kwarg mpi coord z.
    - `do_reset=true`: the optional kwarg to reset the random seed.
    """
    grf3D_expon!()

    @doc raw"""

        grf3D_gauss!(Yf, sf, cl, nh, k_m, nx, ny, nz, dx, dy, dz; me=0, co1=0, co2=0, co3=0, do_reset=true)

    Compute the Gaussian random field `Yf` with isotropic Gaussian covariance function and returns the wall-time `wtime_it`.

    # Arguments
    - `Yf::Data.Array`: the Gaussian random field array.
    - `sf::Data.Number`: the standard deviation.
    - `cl::Any`: the correlation length tuple `(clx, cly, clz)` in `x`, `y` and `z` dimension.
    - `nh::Int`: the number of harmonics.
    - `k_m::Data.Number`: the maximum value of the wave number.
    - `nx::Int`: the number of grid points in `x` dimension.
    - `ny::Int`: the number of grid points in `y` dimension.
    - `nz::Int`: the number of grid points in `z` dimension.
    - `dx::Data.Number`: the numerical grid size in `x` dimension.
    - `dy::Data.Number`: the numerical grid size in `y` dimension.
    - `dz::Data.Number`: the numerical grid size in `z` dimension.
    - `me=0::Int`: the optional kwarg mpi process id.
    - `co1=0::Int`: the optional kwarg mpi coord x.
    - `co2=0::Int`: the optional kwarg mpi coord y.
    - `co3=0::Int`: the optional kwarg mpi coord z.
    - `do_reset=true`: the optional kwarg to reset the random seed.
    """
    grf3D_gauss!()

    @doc raw"""

        generate_grf3D(lx=100.0, ly=100.0, lz=100.0, sf=1.0, cl=(10.0, 8.0, 5.0), k_m=100.0, nh=10000, nx=64, ny=64, nz=64;
                       cov_typ=\"expon\", do_reset=true, do_viz=true, do_save=false)

    Returns the Gaussian random field with exponential `cov_typ=\"expon\"` or Gaussian `cov_typ=\"gauss\"` covariance function for default values.

    # Arguments
    - `lx=100.0`: the `x` dimension domain extend.
    - `ly=100.0`: the `y` dimension domain extend.
    - `lz=100.0`: the `z` dimension domain extend.
    - `sf=1.0`: the standard deviation.
    - `cl=(10.0, 8.0, 5.0)`: the correlation length; `(clx, cly, clz)` tuple (`cov_typ=\"expon\"`), `cl[1]/2.0` (`cov_typ=\"gauss\"`).
    - `k_m=100.0`: the maximum value of the wave number.
    - `nh=10000`: the number of harmonics.
    - `nx=64`: the numerical grid size in `x` dimension.
    - `ny=64`: the numerical grid size in `y` dimension.
    - `nz=64`: the numerical grid size in `y` dimension.
    - `cov_typ=\"expon\"`: the covariance function `\"expon\"` or `\"gauss\"`.
    - `do_viz=true`: the optional kwarg to enable visualisation.
    - `do_save=false`: the optional kwarg to save the random field to disk.
    - `do_reset=true`: the optional kwarg to reset the random seed.
    """
    generate_grf3D()

    using ParallelStencil
    using ParallelStencil.FiniteDifferences3D
    @init_parallel_stencil(Threads, Float64, 3)

    include(joinpath("shared", "grf3D.jl"))
    include(joinpath("shared", "generate_grf3D.jl"))
end


text:
# This file is a part of DataDeps.jl. License is MIT.
module DataDeps

using HTTP
using Reexport
@reexport using SHA

export DataDep, ManualDataDep
export register, resolve, @datadep_str, preupload_check
export unpack

include("errors.jl")
include("types.jl")

include("util.jl")
include("registration.jl")

include("locations.jl")
include("verification.jl")

include("resolution.jl")
include("resolution_automatic.jl")
include("resolution_manual.jl")

include("preupload.jl")

include("fetch_helpers.jl")
include("post_fetch_helpers.jl")

end # module


text:
#Write a program that prints ‘Hello World’ to the screen.

function hello_world()
    println("Hello World")
end

hello_world()

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
module Test

  using ThreadPools

  function sumd(x)
    sumd = 0.
    for i in 1:size(x,1)
      for j in i+1:size(x,1)
        sumd += sqrt((x[i,1]-x[j,1])^2 + (x[i,2]-x[j,2])^2 + (x[i,3]-x[j,3])^2)
      end
    end
    return sumd
  end

  function serial(n)
    s = 0.
    for i in 1:Threads.nthreads()-1
      x = rand(n,3)
      s += sumd(x)
    end
    return s
  end

  function spawn(n)
    nspawn = Threads.nthreads()-1
    t = Vector{Task}(undef,nspawn)
    for i in 1:nspawn
      x = rand(n,3)
      t[i] = ThreadPools.@tspawnat i+1 sumd(x)
    end
    s = 0.
    for i in 1:nspawn
      s += fetch(t[i])
    end
    return s
  end

end

using BenchmarkTools

n = 1000

println(" Serial: ")
@btime Test.serial(n)

println(" Parallel: ")
@btime Test.spawn(n)


text:
using jInvSeismic.BasicFWI
using jInvSeismic.Utils
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare
using Distributed

using DelimitedFiles

@everywhere begin
	using jInvSeismic.BasicFWI
	using jInv.Mesh
	using jInv.Utils
	using jInv.ForwardShare
end

# setup model and attenuation function
nx = 164; nz = 80
domain = [0.0,10.0,0.0,4.0]
Mr = getRegularMesh(domain,[nx,nz])

# the slowness model
m = 3.0*ones(nx,nz);
m[30:80,20:40] .= 5.0; # velocity...
m = 1.0./(m.^2); # convert to slowness squared

# m = readdlm("SEGmodel2Dsalt.dat");
# m = m*1e-3;
# m = m';
# m = 1.0./m^2

# attenuation for BC
padx = 20; padz = 20
a    = 2.0;
xc = getCellCenteredGrid(Mr)
gamma = getHelmholtzABL(Mr,true,[padx;padz],a);

# parameters for the Helmholtz (units in km)
h = Mr.h;
n = [nx; nz;]
omega = 2*pi*[0.5;1.0;2.0;3.0]
nfreq = length(omega)

# generate sources
q = zeros(tuple(n.+1...)); q[padx+1:4:end-padx-1,1] .= 1e4
Q = sdiag(vec(q))
Q = Q[:,(LinearIndices(sum(Q,dims=2) .!= 0))[findall(sum(Q,dims=2) .!= 0)]]
nsrc = size(Q,2)

# receivers
p = zeros(tuple(n.+1...)); p[padx+1:end-padx-1,1] .= 1
P = sdiag(vec(p))
P = P[:,(LinearIndices(sum(P,dims=2) .!= 0))[findall(sum(P,dims=2) .!= 0)]]
nrec = size(P,2)

pFor  = getBasicFWIparam(omega,gamma,Q,P,Mr)
pForp = getBasicFWIparam(omega,gamma,Q,P,Mr,true)

# inversion mesh and forward mesh are the same here
M2Mp = ones(length(pForp))


text:
# This file includes the bit stream generator

export BitGenerator, SymbolGenerator

"""
    $TYPEDEF

Bit stream generator

# Fields

    $TYPEDFIELDS
"""
struct BitGenerator
    "Generated bits"
    bits::Vector{Bool}
    BitGenerator(nbits::Int) = new(rand(Bool, nbits))
end

"""
    $TYPEDEF

Symbol generator

# Fields
    $TYPEDFIELDS
"""
struct SymbolGenerator
    "Generated symbols"
    symbols::Vector{Int}
    SymbolGenerator(ns::Int, M::Int) = new(rand(1 : M, ns))
end


text:
"""
    rand_unitary(N, dims=(N,))

Generate a Haar distributed random unitary operator for a Hilbert space of size `N`. It is possible to specify the subspace dimensions with the `dims` argument. Returns a dense matrix.

# Example
```jldoctest
julia> U = rand_unitary(4,(2,2));

julia> U'*U ≈ qeye(4,(2,2))
true
```
"""
function rand_unitary(N::Integer, dims::Dims=(N,))
    return Operator(rand_unitary(ComplexF64,N),dims)
end

function rand_unitary(::Type{T},n::Integer) where {T<:Compat.LinearAlgebra.BlasComplex}
    # Generate a Haar distributed random unitary matrix
    # ref https://arxiv.org/pdf/math-ph/0609050.pdf
    A = randn(T,n,n)
    f = qr!(A)
    U = Matrix(Compat.LinearAlgebra.QRCompactWYQ(f.factors,f.T))
    @inbounds for i = 1:n
        U[:,i] .*= sign(f.factors[i,i]) # U = Q*Diagonal(diag(R)./abs.(diag(R)))
    end
    return U
end


text:
# test Pose2DPoint2D constraint evaluation function

using LinearAlgebra, Statistics
using RoME
# using KernelDensityEstimate
using RoMEPlotting, Gadfly
using Test

# to export PDF
using Gadfly
using Cairo

##

@testset "Prepare a 2D factor graph with poses and points..." begin

##

N = 100
fg = initfg()

initCov = Matrix(Diagonal([0.03;0.03;0.001]))
odoCov = Matrix(Diagonal([3.0;3.0;0.01]))

# Some starting position
addVariable!(fg, :x0, Pose2, N=N)
initPosePrior = PriorPose2(MvNormal(zeros(3), initCov))
addFactor!(fg,[:x0], initPosePrior)

# and a second pose
addVariable!(fg, :x1, Pose2, N=N)
ppc = Pose2Pose2(MvNormal([50.0;0.0;pi/2], odoCov))
addFactor!(fg, [:x0; :x1], ppc)

# test evaluation of pose pose constraint
pts = approxConv(fg, :x0x1f1, :x1)

# perform inference
tree = solveTree!(fg)

# check that yaw is working
addVariable!(fg, :x2, Pose2, N=N)
ppc = Pose2Pose2(MvNormal([50.0;0.0;0.0], odoCov))
addFactor!(fg, [:x1;:x2], ppc)

# new landmark
l1 = addVariable!(fg, :l1, Point2, N=N)
# and pose to landmark constraint
rhoZ1 = norm([10.0;0.0])
ppr = Pose2Point2BearingRange{Uniform, Normal}(Uniform(-pi,pi),Normal(rhoZ1,1.0))
addFactor!(fg, [:x0;:l1], ppr)

# add a prior to landmark
pp2 = PriorPoint2(MvNormal([10.0;0.0], Matrix(Diagonal([1.0;1.0]))))

f5 = addFactor!(fg,[:l1], pp2)

# do inference
tree = solveTree!(fg)
# ensureAllInitialized!(fg)
# tree = wipeBuildNewTree!(fg)
# [inferOverTree!(fg, tree, N=N) for i in 1:2]

println("test Pose2D plotting")

plotSLAM2DPoses(fg);
plotSLAM2DLandmarks(fg);
plotSLAM2D(fg);

p1 = getBelief(fg, :l1)
# pts = getVal(fg, :l1)
# p1= kde!(pts)
p1c = getBelief(fg, :x0)

plotKDE( p1 , dimLbls=["x";"y";"z"]) # |> PDF("/tmp/test.pdf")

plotKDE( [p1c;p1] , dimLbls=["x";"y";"z"],c=["red";"black"],levels=3, dims=[1;2])
# plotKDE( [marginal(p1c,[1;2]);marginal(p1,[1;2])] , dimLbls=["x";"y";"z"],c=["red";"black"],levels=3, dims=[1;2])
p1c = deepcopy(p1)

plotKDE( marginal(getBelief(getVariable(fg, :x2)),[1;2]) , dimLbls=["x";"y";"z"])

axis = [[1.5;3.5]';[-1.25;1.25]';[-1.0;1.0]']

# @warn "Reinsert draw test.pdf"
plotKDE( p1, dimLbls=["x";"y";"z"], axis=axis) |> PDF("/tmp/test.pdf",30cm,20cm)

#
Base.rm("/tmp/test.pdf")

##

end

#


text:
# Tom Anderson
# vdiv_filter.jl
# MIT License
# 10/30/2016
# Function library for an NGSpice post processor using the Julia repl

using Winston
import Winston.*

global const m = 1e-3
global const u = 1e-6
global const n = 1e-9
global const p = 1e-12
global const f = 1e-15
global const k = 1e3
global const K = 1e3
global const meg = 1e6
global const Meg = 1e6
global const M = 1e6
global const MEG = 1e6
global const G = 1e9
global const g = 1e9

vdiv(z1, z2) = (z2 + 0im) ./ (z1 + z2)
para(z1, z2) = (z1 + 0im) .* vdiv(z1, z2)

vdiv(z1::Array{Float64,1}, z2::Array{Float64,1}) = (z2 + 0im) ./ (z1 + z2)
para(z1::Array{Float64,1}, z2::Array{Float64,1}) = (z1 + 0im) .* vdiv(z1, z2)
vdiv(z1::Array{Complex,1}, z2::Array{Complex,1}) = z2 ./ (z1 + z2)
para(z1::Array{Complex,1}, z2::Array{Complex,1}) = z1 .* vdiv(z1, z2)
vdiv(z1::Array{Complex,1}, z2::Array{Float64,1}) = z2 ./ (z1 + z2)
para(z1::Array{Complex,1}, z2::Array{Float64,1}) = z1 .* vdiv(z1, z2)
vdiv(z1::Array{Float64,1}, z2::Array{Complex,1}) = z2 ./ (z1 + z2)
para(z1::Array{Float64,1}, z2::Array{Complex,1}) = z1 .* vdiv(z1, z2)

vdiv(z1::Float64, z2::Float64) = (z2 + 0im) / (z1 + z2)
para(z1::Float64, z2::Float64) = (z1 + 0im) * vdiv(z1, z2)
vdiv(z1::Complex, z2::Complex) = z2 / (z1 + z2)
para(z1::Complex, z2::Complex) = z1 * vdiv(z1, z2)
vdiv(z1::Complex, z2::Float64) = z2 / (z1 + z2)
para(z1::Complex, z2::Float64) = z1 * vdiv(z1, z2)
vdiv(z1::Float64, z2::Complex) = z2 / (z1 + z2)
para(z1::Float64, z2::Complex) = z1 * vdiv(z1, z2)

magn(x) =hypot(real(x), imag(x))

db(x) = 20.0 * log10(magn(x))
logsweep(x1, x2, npts) = logspace(log10(x1), log10(x2), npts)
zc(c,freq) = -(1im / (2π * c)) ./ freq
zl(l,freq) = 1im * 2π * l .* freq
zr(r,freq) = ones(Integer(length(freq))) * r
Base.max(x::Array{Float64,1})=maximum(x)
Base.min(x::Array{Float64,1})=minimum(x)

max_stars_repo_name: tomacorp/PyPost

max_issues_repo_name: tomacorp/PyPost

max_forks_repo_name: tomacorp/PyPost


text:
push!(LOAD_PATH, "../src/")
using Documenter, GeneralLedger

makedocs(sitename="GeneralLedger.jl",
	  pages=[
		   "Objective" => "index.md"
		   "Introduction" => "introduction.md"
		   "Tutorials" => [
				  "About GL" => "tutorials/aboutgl.md"
				  "GL Processes" => "tutorials/glprocesses.md"
				  "ERD" => "tutorials/erd.md"
				  "Installing Julia" => "tutorials/installation.md"
				  "ELT vs ETL" => "tutorials/elt.md"
				  "Self-Service Data Analytic" => "tutorials/selfservice.md"
				  "Visualizations, Buttons, sliders, filters, n-D plots, plots vs graphs" => "tutorials/plots.md"
				  "p-value, null hypothesis and real time analytic" => "tutorials/analytic.md"
				  "Time Series, Impact analysis" => "tutorials/timeseries.md"
				  "ML for GL" => "tutorials/mlforgl.md"
				  "NLP for GL" => "tutorials/nlp.md"
				  "Graph Theory / Network Science" => "tutorials/graph.md"
				  ]
			"API" => "api.md"
		   ],
	  )


text:
using RBNF
using PrettyPrint
export lex, parse_qasm, Token, print_ast, string_ast

""" Check https://arxiv.org/pdf/1707.03429.pdf for grammar specification
"""
struct QASMLang end

second((a, b)) = b
second(vec::V) where V <: AbstractArray = vec[2]

struct Struct_bin
    l
    op :: RBNF.Token
    r
end

RBNF.@parser QASMLang begin
    # define ignorances
    ignore{space}

    @grammar
    # define grammars
    mainprogram := ["OPENQASM", ver=real, ';', prog=program]
    program     = statement{*}
    statement   = (decl | gate | opaque | qop | ifstmt | barrier)
    # stmts
    ifstmt      := ["if", '(', l=id, "==", r=nninteger, ')', body=qop]
    opaque      := ["opaque", id=id, ['(', [arglist1=idlist].?, ')'].? , arglist2=idlist, ';']
    barrier     := ["barrier", value=mixedlist]
    decl        := [regtype=("qreg" | "creg"), id=id, '[', int=nninteger, ']', ';']

    # gate
    gate        := [decl=gatedecl, [goplist=goplist].?, '}']
    gatedecl    := ["gate", id=id, ['(', [args=idlist].?, ')'].?, (outs=idlist), '{']

    goplist     = (uop |barrier_ids){*}
    barrier_ids := ["barrier", ids=idlist, ';'] # not impl
    # qop
    qop         = (uop | measure | reset)
    reset       := ["reset", arg=argument, ';'] # not impl
    measure     := ["measure", arg1=argument, "->", arg2=argument, ';'] # not impl

    uop         = (iduop | u | cx)
    iduop      := [gate_name=id, ['(', [args=explist].?, ')'].?, outs=mixedlist, ';']
    u          := ['U', '(', in1=exp, ',', in2=exp, ',', in3 = exp, ')', out=argument, ';']
    cx         := ["CX", out1=argument, ',', out2=argument, ';']

    idlist     := [hd=id, [',', tl=idlist].?]

    mixedlist  := [hd=argument, [',', tl=mixedlist].?]

    argument   := [id=id, ['[', (arg=nninteger), ']'].?]

    explist    := [hd=exp, [',', tl=explist].?]
    pi         := "pi"
    atom       =  (real | nninteger | pi | id | fnexp) | (['(', exp, ')'] % second) | neg
    fnexp      := [fn=fn, '(', arg=exp, ')']
    neg        := ['-', value=exp]

    exp        = [l=mul,  [op=('+' |'-'), r=exp].?] => _.op === nothing ? _.l : Struct_bin(_.l, _.op, _.r)
    mul        = [l=atom, [op=('*' | '/'), r=mul].?] => _.op === nothing ? _.l : Struct_bin(_.l, _.op, _.r)
    fn         = ("sin" | "cos" | "tan" | "exp" | "ln" | "sqrt")

    # define tokens
    @token
    id        := r"\G[a-z]{1}[A-Za-z0-9_]*"
    real      := r"\G([0-9]+\.[0-9]*|[0-9]*\.[0.9]+)([eE][-+]?[0-9]+)?"
    nninteger := r"\G([1-9]+[0-9]*|0)"
    space     := r"\G\s+"
end

Token = RBNF.Token

function lex(src :: String)
    RBNF.runlexer(QASMLang, src)
end

function parse_qasm(tokens :: Vector{Token{A} where A})
    ast, ctx = RBNF.runparser(mainprogram, tokens)
    ast
end

print_ast = PrettyPrint.pprint
string_ast = PrettyPrint.pformat