text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

using Base.Math: nan_dom_err

struct AmosException <: Exception
    id::Int32
end

function Base.showerror(io::IO, ex::AmosException)
    print(io, "AmosException with id $(ex.id): ")
    if ex.id == 0
        print(io, "normal return, computation complete.")
    elseif ex.id == 1
        print(io, "input error.")
    elseif ex.id == 2
        print(io, "overflow.")
    elseif ex.id == 3
        print(io, "input argument magnitude large, less than half machine accuracy loss by argument reduction.")
    elseif ex.id == 4
        print(io, "input argument magnitude too large, complete loss of accuracy by argument reduction.")
    elseif ex.id == 5
        print(io, "algorithm termination condition not met.")
    else
        print(io, "invalid error flag.")
    end
end

## Airy functions
function _airy(z::Complex{Float64}, id::Int32, kode::Int32)
    ai1, ai2 = Ref{Float64}(), Ref{Float64}()
    ae1, ae2 = Ref{Int32}(), Ref{Int32}()

    ccall((:zairy_,libopenspecfun), Cvoid,
          (Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
           Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
           real(z), imag(z), id, kode,
           ai1, ai2, ae1, ae2)

    if ae2[] == 0 || ae2[] == 3 # ignore underflow and less than half machine accuracy loss
        return complex(ai1[], ai2[])
    else
        throw(AmosException(ae2[]))
    end
end

function _biry(z::Complex{Float64}, id::Int32, kode::Int32)
    ai1, ai2 = Ref{Float64}(), Ref{Float64}()
    ae1 = Ref{Int32}()

    ccall((:zbiry_,libopenspecfun), Cvoid,
          (Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
           Ref{Float64}, Ref{Float64}, Ref{Int32}),
           real(z), imag(z), id, kode,
           ai1, ai2, ae1)

    if ae1[] == 0 || ae1[] == 3 # ignore less than half machine accuracy loss
        return complex(ai1[], ai2[])
    else
        throw(AmosException(ae1[]))
    end
end

"""
    airyai(x)

Airy function of the first kind ``\\operatorname{Ai}(x)``.

External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function)

See also: [`airyaix`](@ref), [`airyaiprime`](@ref), [`airybi`](@ref)
"""
function airyai end
airyai(z::Complex{Float64}) = _airy(z, Int32(0), Int32(1))

"""
    airyaiprime(x)

Derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x)``.

External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function)

See also: [`airyaiprimex`](@ref), [`airyai`](@ref), [`airybi`](@ref)
"""
function airyaiprime end
airyaiprime(z::Complex{Float64}) =  _airy(z, Int32(1), Int32(1))

"""
    airybi(x)

Airy function of the second kind ``\\operatorname{Bi}(x)``.

External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function)

See also: [`airybix`](@ref), [`airybiprime`](@ref),  [`airyai`](@ref)
"""
function airybi end
airybi(z::Complex{Float64}) = _biry(z, Int32(0), Int32(1))

"""
    airybiprime(x)

Derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x)``.

External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function)

See also: [`airybiprimex`](@ref), [`airybi`](@ref), [`airyai`](@ref)
"""
function airybiprime end
airybiprime(z::Complex{Float64}) = _biry(z, Int32(1), Int32(1))

"""
    airyaix(x)

Scaled Airy function of the first kind ``\\operatorname{Ai}(x) e^{\\frac{2}{3} x
\\sqrt{x}}``.  Throws `DomainError` for negative `Real` arguments.

External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function)

See also: [`airyai`](@ref), [`airyaiprime`](@ref), [`airybi`](@ref)
"""
function airyaix end
airyaix(z::Complex{Float64}) = _airy(z, Int32(0), Int32(2))

"""
    airyaiprimex(x)

Scaled derivative of the Airy function of the first kind ``\\operatorname{Ai}'(x)
e^{\\frac{2}{3} x \\sqrt{x}}``.  Throws `DomainError` for negative `Real` arguments.

External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function)

See also: [`airyaiprime`](@ref), [`airyai`](@ref), [`airybi`](@ref)
"""
function airyaiprimex end
airyaiprimex(z::Complex{Float64}) =  _airy(z, Int32(1), Int32(2))

"""
    airybix(x)

Scaled Airy function of the second kind ``\\operatorname{Bi}(x) e^{- \\left| \\operatorname{Re} \\left( \\frac{2}{3} x \\sqrt{x} \\right) \\right|}``.

External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function)

See also: [`airybi`](@ref), [`airybiprime`](@ref), [`airyai`](@ref)
"""
function airybix end
airybix(z::Complex{Float64}) = _biry(z, Int32(0), Int32(2))

"""
    airybiprimex(x)

Scaled derivative of the Airy function of the second kind ``\\operatorname{Bi}'(x) e^{- \\left| \\operatorname{Re} \\left( \\frac{2}{3} x \\sqrt{x} \\right) \\right|}``.

External links: [DLMF](https://dlmf.nist.gov/9.2), [Wikipedia](https://en.wikipedia.org/wiki/Airy_function)

See also: [`airybiprime`](@ref), [`airybi`](@ref), [`airyai`](@ref)
"""
function airybiprimex end
airybiprimex(z::Complex{Float64}) = _biry(z, Int32(1), Int32(2))

for afn in (:airyai, :airyaiprime, :airybi, :airybiprime,
            :airyaix, :airyaiprimex, :airybix, :airybiprimex)
    @eval begin
        $afn(z::Complex) = $afn(float(z))
        $afn(z::Complex{Float16}) = Complex{Float16}($afn(Complex{Float32}(z)))
        $afn(z::Complex{Float32}) = Complex{Float32}($afn(Complex{Float64}(z)))
        $afn(z::Complex{<:AbstractFloat}) = throw(MethodError($afn,(z,)))
    end
    if afn in (:airyaix, :airyaiprimex)
        @eval $afn(x::Real) = x < 0 ? throw(DomainError(x, "`x` must be nonnegative.")) : real($afn(complex(float(x))))
    else
        @eval $afn(x::Real) = real($afn(complex(float(x))))
    end
end

function airyai(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_ai, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
    return z
end

## Bessel functions

# besselj0, besselj1, bessely0, bessely1
for jy in ("j","y"), nu in (0,1)
    jynu = Expr(:quote, Symbol(jy,nu))
    jynuf = Expr(:quote, Symbol(jy,nu,"f"))
    bjynu = Symbol("bessel",jy,nu)
    if jy == "y"
        @eval begin
            $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm),  Float64, (Float64,), x), x)
            $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x)
            $bjynu(x::Float16) = Float16($bjynu(Float32(x)))
        end
    else
        @eval begin
            $bjynu(x::Float64) = ccall(($jynu,libm),  Float64, (Float64,), x)
            $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
            $bjynu(x::Float16) = Float16($bjynu(Float32(x)))
        end
    end
    @eval begin
        $bjynu(x::Real) = $bjynu(float(x))
        $bjynu(x::Complex) = $(Symbol("bessel",jy))($nu,x)
    end
end

function _besselh(nu::Float64, k::Int32, z::Complex{Float64}, kode::Int32)
    ai1, ai2 = Ref{Float64}(), Ref{Float64}()
    ae1, ae2 = Ref{Int32}(), Ref{Int32}()

    ccall((:zbesh_,libopenspecfun), Cvoid,
           (Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}, Ref{Int},
            Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
            real(z), imag(z), nu, kode, k, 1,
            ai1, ai2, ae1, ae2)

    if ae2[] == 0 || ae2[] == 3
        return complex(ai1[], ai2[])
    else
        throw(AmosException(ae2[]))
    end
end

function _besseli(nu::Float64, z::Complex{Float64}, kode::Int32)
    ai1, ai2 = Ref{Float64}(), Ref{Float64}()
    ae1, ae2 = Ref{Int32}(), Ref{Int32}()

    ccall((:zbesi_,libopenspecfun), Cvoid,
          (Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
           Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
           real(z), imag(z), nu, kode, 1,
           ai1, ai2, ae1, ae2)

    if ae2[] == 0 || ae2[] == 3
        return complex(ai1[], ai2[])
    else
        throw(AmosException(ae2[]))
    end
end

function _besselj(nu::Float64, z::Complex{Float64}, kode::Int32)
    ai1, ai2 = Ref{Float64}(), Ref{Float64}()
    ae1, ae2 = Ref{Int32}(), Ref{Int32}()

    ccall((:zbesj_,libopenspecfun), Cvoid,
          (Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
           Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
           real(z), imag(z), nu, kode, 1,
           ai1, ai2, ae1, ae2)

    if ae2[] == 0 || ae2[] == 3
        return complex(ai1[], ai2[])
    else
        throw(AmosException(ae2[]))
    end
end

function _besselk(nu::Float64, z::Complex{Float64}, kode::Int32)
    ai1, ai2 = Ref{Float64}(), Ref{Float64}()
    ae1, ae2 = Ref{Int32}(), Ref{Int32}()

    ccall((:zbesk_,libopenspecfun), Cvoid,
          (Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
           Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
           real(z), imag(z), nu, kode, 1,
           ai1, ai2, ae1, ae2)

    if ae2[] == 0 || ae2[] == 3
        return complex(ai1[], ai2[])
    else
        throw(AmosException(ae2[]))
    end
end

function _bessely(nu::Float64, z::Complex{Float64}, kode::Int32)
    ai1, ai2 = Ref{Float64}(), Ref{Float64}()
    ae1, ae2 = Ref{Int32}(), Ref{Int32}()
    wrk1, wrk2 = Ref{Float64}(), Ref{Float64}()

    ccall((:zbesy_,libopenspecfun), Cvoid,
          (Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
           Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Int32}),
           real(z), imag(z), nu, kode, 1,
           ai1, ai2, ae1, wrk1, wrk2, ae2)

    if ae2[] == 0 || ae2[] == 3
        return complex(ai1[], ai2[])
    else
        throw(AmosException(ae2[]))
    end
end

"""
    besselh(nu, [k=1,] x)

Bessel function of the third kind of order `nu` (the Hankel function). `k` is either 1 or 2,
selecting [`hankelh1`](@ref) or [`hankelh2`](@ref), respectively.
`k` defaults to 1 if it is omitted.

External links: [DLMF](https://dlmf.nist.gov/10.2.5) and [DLMF](https://dlmf.nist.gov/10.2.6), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1)

See also: [`besselhx`](@ref) for an exponentially scaled variant.
"""
function besselh end

function besselh(nu::Float64, k::Integer, z::Complex{Float64})
    if nu < 0
        s = (k == 1) ? 1 : -1
        return _besselh(-nu,Int32(k),z,Int32(1)) * complex(cospi(nu),-s*sinpi(nu))
    end
    return _besselh(nu,Int32(k),z,Int32(1))
end

function besselh(nu::Float64, k::Integer, x::Float64)
    # Given that x is real, Jnu(x) and Ynu(x) are also real.
    if k == 1
        return complex(besselj(nu, x), bessely(nu, x))
    elseif k == 2
        return complex(besselj(nu, x), -bessely(nu, x))
    else
        # We emulate ZBESH's behaviour
        throw(AmosException(1))
    end
end

"""
    besselhx(nu, [k=1,] z)

Compute the scaled Hankel function ``\\exp(∓iz) H_ν^{(k)}(z)``, where
``k`` is 1 or 2, ``H_ν^{(k)}(z)`` is `besselh(nu, k, z)`, and ``∓`` is
``-`` for ``k=1`` and ``+`` for ``k=2``.  `k` defaults to 1 if it is omitted.

The reason for this function is that ``H_ν^{(k)}(z)`` is asymptotically
proportional to ``\\exp(∓iz)/\\sqrt{z}`` for large ``|z|``, and so the
[`besselh`](@ref) function is susceptible to overflow or underflow
when `z` has a large imaginary part.  The `besselhx` function cancels this
exponential factor (analytically), so it avoids these problems.

External links: [DLMF](https://dlmf.nist.gov/10.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1)

See also: [`besselh`](@ref)
"""
function besselhx end

function besselhx(nu::Float64, k::Integer, z::Complex{Float64})
    if nu < 0
        s = (k == 1) ? 1 : -1
        return _besselh(-nu,Int32(k),z,Int32(2)) * complex(cospi(nu),-s*sinpi(nu))
    end
    return _besselh(nu,Int32(k),z,Int32(2))
end

function besseli(nu::Float64, z::Complex{Float64})
    if nu < 0
        if isinteger(nu)
            return _besseli(-nu,z,Int32(1))
        else
            return _besseli(-nu,z,Int32(1)) - 2_besselk(-nu,z,Int32(1))*sinpi(nu)/pi
        end
    else
        return _besseli(nu,z,Int32(1))
    end
end

function besselix(nu::Float64, z::Complex{Float64})
    if nu < 0
        if isinteger(nu)
            return _besseli(-nu,z,Int32(2))
        else
            return _besseli(-nu,z,Int32(2)) - 2_besselk(-nu,z,Int32(2))*exp(-abs(real(z))-z)*sinpi(nu)/pi
        end
    else
        return _besseli(nu,z,Int32(2))
    end
end

function besselj(nu::Float64, z::Complex{Float64})
    if nu < 0
        if isinteger(nu)
            return _besselj(-nu,z,Int32(1))*cospi(nu)
        else
            return _besselj(-nu,z,Int32(1))*cospi(nu) + _bessely(-nu,z,Int32(1))*sinpi(nu)
        end
    else
        return _besselj(nu,z,Int32(1))
    end
end

besselj(nu::Cint, x::Float64) = ccall((:jn, libm), Float64, (Cint, Float64), nu, x)
besselj(nu::Cint, x::Float32) = ccall((:jnf, libm), Float32, (Cint, Float32), nu, x)

function besseljx(nu::Float64, z::Complex{Float64})
    if nu < 0
        if isinteger(nu)
            return _besselj(-nu,z,Int32(2))*cospi(nu)
        else
            return _besselj(-nu,z,Int32(2))*cospi(nu) + _bessely(-nu,z,Int32(2))*sinpi(nu)
        end
    else
        return _besselj(nu,z,Int32(2))
    end
end

besselk(nu::Float64, z::Complex{Float64}) = _besselk(abs(nu), z, Int32(1))

besselkx(nu::Float64, z::Complex{Float64}) = _besselk(abs(nu), z, Int32(2))

function bessely(nu::Cint, x::Float64)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    end
    ccall((:yn, libm), Float64, (Cint, Float64), nu, x)
end
function bessely(nu::Cint, x::Float32)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    end
    ccall((:ynf, libm), Float32, (Cint, Float32), nu, x)
end

function bessely(nu::Float64, z::Complex{Float64})
    if nu < 0
        return _bessely(-nu,z,Int32(1))*cospi(nu) - _besselj(-nu,z,Int32(1))*sinpi(nu)
    else
        return _bessely(nu,z,Int32(1))
    end
end

function besselyx(nu::Float64, z::Complex{Float64})
    if nu < 0
        return _bessely(-nu,z,Int32(2))*cospi(nu) - _besselj(-nu,z,Int32(2))*sinpi(nu)
    else
        return _bessely(nu,z,Int32(2))
    end
end

"""
    besseli(nu, x)

Modified Bessel function of the first kind of order `nu`, ``I_\\nu(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.25.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1)

See also: [`besselix(nu,x)`](@ref SpecialFunctions.besselix), [`besselj(nu,x)`](@ref SpecialFunctions.besselj), [`besselk(nu,x)`](@ref SpecialFunctions.besselk)
"""
function besseli(nu::Real, x::AbstractFloat)
    if x < 0 && !isinteger(nu)
        throw(DomainError(x, "`x` must be nonnegative and `nu` must be an integer."))
    end
    real(besseli(float(nu), complex(x)))
end

"""
    besselix(nu, x)

Scaled modified Bessel function of the first kind of order `nu`, ``I_\\nu(x) e^{- | \\operatorname{Re}(x) |}``.

External links: [DLMF](https://dlmf.nist.gov/10.25.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1)

See also: [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselj(nu,x)`](@ref SpecialFunctions.besselj), [`besselk(nu,x)`](@ref SpecialFunctions.besselk)
"""
function besselix(nu::Real, x::AbstractFloat)
    if x < 0 && !isinteger(nu)
        throw(DomainError(x, "`x` must be nonnegative and `nu` must be an integer."))
    end
    real(besselix(float(nu), complex(x)))
end

"""
    besselj(nu, x)

Bessel function of the first kind of order `nu`, ``J_\\nu(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1)

See also: [`besseljx(nu,x)`](@ref SpecialFunctions.besseljx), [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselk(nu,x)`](@ref SpecialFunctions.besselk)
"""
function besselj(nu::Real, x::AbstractFloat)
    if isinteger(nu)
        if typemin(Cint) <= nu <= typemax(Cint)
            return besselj(Cint(nu), x)
        end
    elseif x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    end
    real(besselj(float(nu), complex(x)))
end

"""
    besseljx(nu, x)

Scaled Bessel function of the first kind of order `nu`, ``J_\\nu(x) e^{- | \\operatorname{Im}(x) |}``.

External links: [DLMF](https://dlmf.nist.gov/10.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1)

See also: [`besselj(nu,x)`](@ref SpecialFunctions.besselj), [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselk(nu,x)`](@ref SpecialFunctions.besselk)
"""
function besseljx(nu::Real, x::AbstractFloat)
    if x < 0 && !isinteger(nu)
        throw(DomainError(x, "`x` must be nonnegative and `nu` must be an integer."))
    end
    real(besseljx(float(nu), complex(x)))
end

"""
    besselk(nu, x)

Modified Bessel function of the second kind of order `nu`, ``K_\\nu(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.25.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1)

See also: See also: [`besselkx(nu,x)`](@ref SpecialFunctions.besselkx), [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselj(nu,x)`](@ref SpecialFunctions.besselj)
"""
function besselk(nu::Real, x::AbstractFloat)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    elseif x == 0
        return oftype(x, Inf)
    end
    real(besselk(float(nu), complex(x)))
end

"""
    besselkx(nu, x)

Scaled modified Bessel function of the second kind of order `nu`, ``K_\\nu(x) e^x``.

External links: [DLMF](https://dlmf.nist.gov/10.25.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1)

See also: [`besselk(nu,x)`](@ref SpecialFunctions.besselk), [`besseli(nu,x)`](@ref SpecialFunctions.besseli), [`besselj(nu,x)`](@ref SpecialFunctions.besselj)
"""
function besselkx(nu::Real, x::AbstractFloat)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    elseif x == 0
        return oftype(x, Inf)
    end
    real(besselkx(float(nu), complex(x)))
end

"""
    bessely(nu, x)

Bessel function of the second kind of order `nu`, ``Y_\\nu(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1)

See also [`besselyx(nu,x)`](@ref SpecialFunctions.besselyx) for a scaled variant.
"""
function bessely(nu::Real, x::AbstractFloat)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    elseif isinteger(nu) && typemin(Cint) <= nu <= typemax(Cint)
        return bessely(Cint(nu), x)
    end
    real(bessely(float(nu), complex(x)))
end

"""
    besselyx(nu, x)

Scaled Bessel function of the second kind of order `nu`,
``Y_\\nu(x) e^{- | \\operatorname{Im}(x) |}``.

External links: [DLMF](https://dlmf.nist.gov/10.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1)

See also [`bessely(nu,x)`](@ref SpecialFunctions.bessely)
"""
function besselyx(nu::Real, x::AbstractFloat)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    end
    real(besselyx(float(nu), complex(x)))
end

for f in ("i", "ix", "j", "jx", "k", "kx", "y", "yx")
    bfn = Symbol("bessel", f)
    @eval begin
        $bfn(nu::Real, x::Real) = $bfn(nu, float(x))
        function $bfn(nu::Real, z::Complex)
            Tf = promote_type(float(typeof(nu)),float(typeof(real(z))))
            $bfn(Tf(nu), Complex{Tf}(z))
        end
        $bfn(nu::Float16, x::Complex{Float16}) = Complex{Float16}($bfn(Float32(nu), Complex{Float32}(x)))
        $bfn(nu::Float32, x::Complex{Float32}) = Complex{Float32}($bfn(Float64(nu), Complex{Float64}(x)))
        $bfn(k::T, z::Complex{T}) where {T<:AbstractFloat} = throw(MethodError($bfn,(k,z)))
    end
end

for bfn in (:besselh, :besselhx)
    @eval begin
        $bfn(nu, z) = $bfn(nu, 1, z)
        $bfn(nu::Real, k::Integer, x::Real) = $bfn(float(nu), k, float(x))
        $bfn(nu::AbstractFloat, k::Integer, x::AbstractFloat) = $bfn(float(nu), k, complex(x))
        function $bfn(nu::Real, k::Integer, z::Complex)
            Tf = promote_type(float(typeof(nu)),float(typeof(real(z))))
            $bfn(Tf(nu), k, Complex{Tf}(z))
        end
        $bfn(nu::Float16, k::Integer, x::Complex{Float16}) = Complex{Float16}($bfn(Float32(nu), k, Complex{Float32}(x)))
        $bfn(nu::Float32, k::Integer, x::Complex{Float32}) = Complex{Float32}($bfn(Float64(nu), k, Complex{Float64}(x)))
        $bfn(nu::T, k::Integer, z::Complex{T}) where {T<:AbstractFloat} = throw(MethodError($bfn,(nu,k,z)))
    end
end

besselh(nu::Float16, k::Integer, x::Float16) = Complex{Float16}(besselh(Float32(nu), k, Float32(x)))
besselh(nu::Float32, k::Integer, x::Float32) = Complex{Float32}(besselh(Float64(nu), k, Float64(x)))

"""
    besselj0(x)

Bessel function of the first kind of order 0, ``J_0(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1)

See also: [`besselj(nu,x)`](@ref SpecialFunctions.besselj)
"""
function besselj0(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_j0, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
    return z
end

"""
    besselj1(x)

Bessel function of the first kind of order 1, ``J_1(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J%CE%B1)

See also: [`besselj(nu,x)`](@ref SpecialFunctions.besselj)
"""
function besselj1(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_j1, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
    return z
end

function besselj(n::Integer, x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_jn, :libmpfr), Int32, (Ref{BigFloat}, Clong, Ref{BigFloat}, Int32), z, n, x, ROUNDING_MODE[])
    return z
end

"""
    bessely0(x)

Bessel function of the second kind of order 0, ``Y_0(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1)

See also: [`bessely(nu,x)`](@ref SpecialFunctions.bessely)
"""
function bessely0(x::BigFloat)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    end
    z = BigFloat()
    ccall((:mpfr_y0, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
    return z
end

"""
    bessely1(x)

Bessel function of the second kind of order 1, ``Y_1(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.2.3), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_second_kind:_Y%CE%B1)

See also: [`bessely(nu,x)`](@ref SpecialFunctions.bessely)
"""
function bessely1(x::BigFloat)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    end
    z = BigFloat()
    ccall((:mpfr_y1, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
    return z
end

function bessely(n::Integer, x::BigFloat)
    if x < 0
        throw(DomainError(x, "`x` must be nonnegative."))
    end
    z = BigFloat()
    ccall((:mpfr_yn, :libmpfr), Int32, (Ref{BigFloat}, Clong, Ref{BigFloat}, Int32), z, n, x, ROUNDING_MODE[])
    return z
end

"""
    hankelh1(nu, x)

Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.2.5), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1)

See also: [`hankelh1x`](@ref SpecialFunctions.hankelh1x)
"""
hankelh1(nu, z) = besselh(nu, 1, z)

"""
    hankelh2(nu, x)

Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x)``.

External links: [DLMF](https://dlmf.nist.gov/10.2.6), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1)

See also: [`hankelh2x(nu,x)`](@ref SpecialFunctions.hankelh2x)
"""
hankelh2(nu, z) = besselh(nu, 2, z)

"""
    hankelh1x(nu, x)

Scaled Bessel function of the third kind of order `nu`, ``H^{(1)}_\\nu(x) e^{-x i}``.

External links: [DLMF](https://dlmf.nist.gov/10.2.5), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1)

See also: [`hankelh1`](@ref SpecialFunctions.hankelh1)
"""
hankelh1x(nu, z) = besselhx(nu, 1, z)

"""
    hankelh2x(nu, x)

Scaled Bessel function of the third kind of order `nu`, ``H^{(2)}_\\nu(x) e^{x i}``.

External links: [DLMF](https://dlmf.nist.gov/10.2.6), [Wikipedia](https://en.wikipedia.org/wiki/Bessel_function#Hankel_functions:_H(1)%CE%B1,_H(2)%CE%B1)

See also: [`hankelh2(nu,x)`](@ref SpecialFunctions.hankelh2)
"""
hankelh2x(nu, z) = besselhx(nu, 2, z)


text:
"""
	$(SIGNATURES)

Grid for possible no of courses attempted at the start of each period.
Also contains info about no of courses that can be taken in each period.

This can be done in arbitrary units (e.g. one model course = one data course)

nCompleted = no of courses taken at start of `t`
nTried = no of courses tried IN period `t`
"""
struct CourseGrid
    # No of periods
    T :: ncInt
    # No of courses that can be taken in each period
    nCourseV :: Vector{ncInt}
    # For each period: possible no of courses at start of that period
    nGridV :: Vector{Vector{ncInt}}
end

## Possible choices for `n` at `t`
# No of courses tried this period.
n_tried_grid(cg :: CourseGrid, t :: Integer) = cg.nCourseV;
n_n_tried(cg :: CourseGrid, t :: Integer) = length(cg.nCourseV);

## Possible values for no of courses taken at start of t
function n_completed_grid(cg :: CourseGrid, t :: Integer)
    @assert t <= cg.T
    return cg.nGridV[t]
end

function n_n_completed(cg :: CourseGrid, t :: Integer)
    return length(cg.nGridV[t])
end

# Constructor when no of courses is the same for all periods
# Especially efficient when each entry in `nCourseV` is a multiple of one number
function make_course_grid(T :: Integer, nCourseV :: Vector{T1}) where T1 <: Integer
    nGridV = make_ncompleted_grid(T, nCourseV);
    return CourseGrid(T, nCourseV, nGridV)
end

function make_test_course_grid(T :: Integer)
    return make_course_grid(T, [5, 8, 10]);
end

# ----------------  Retrieve

"""
    $(SIGNATURES)

Convert no of courses to index values on the grid.
If any `nV` are not on the grid: index is set to 0.
Always returns a `Vector`, but the method that is called with scalar `nV` does return a scalar.
Performance critical.
"""
function ncompleted_to_indices(cg :: CourseGrid, t :: Integer, nV)
    # This is much faster than the hand-written approach below
    if dbgHigh
        @assert issorted(nV)  "Sorted nV required for `find_indices`"
    end
    idxOutV = find_indices(nV, n_completed_grid(cg, t));
    return idxOutV
end

# input and output are scalar
ncompleted_to_indices(cg :: CourseGrid, t :: Integer, n :: Integer) =
    find_index(n, n_completed_grid(cg, t));

"""
    $(SIGNATURES)

Grid: possible no of courses completed at start of each period.
Builds up the grids that are stored inside the object.
"""
function make_ncompleted_grid(T :: Integer, nCourseV :: Vector{T1}) where T1 <: Integer

    gridV = Vector{Vector{ncInt}}(undef, T);
    # No courses attempted at start of first period
    gridV[1] = zeros(ncInt, 1);
    for t = 2 : T
        tGridV = Vector{ncInt}();
        prevGridV = gridV[t-1];
        for i1 = 1 : length(prevGridV)
            append!(tGridV, prevGridV[i1] .+ nCourseV);
        end
        gridV[t] = sort(unique(tGridV));
    end

    return gridV
end

"""
    $(SIGNATURES)

Law of motion for no of courses attempted.
Inputs are indices into course grids (which become the states of the student
problem).
This is for `idx_nTakenV` courses taken in `t` yielding grid points for `t+1`.
Performance critical.

# Arguments
- nCompleted
    No of courses taken at start of t.
"""
update_ncourses(cg :: CourseGrid, t :: Integer, nCompleted :: Integer, idx_nTriedV) =
    ncompleted_to_indices(cg, t + 1, nCompleted .+ cg.nCourseV[idx_nTriedV]);

# Same, for all possible choices of `nCompleted` at `t`
update_ncourses(cg :: CourseGrid, t :: Integer, nCompleted :: Integer) =
    ncompleted_to_indices(cg, t+1, nCompleted .+ cg.nCourseV);

# function update_ncourses(cg :: CourseGrid, t :: Integer, nCompleted :: Integer)
#     return update_ncourses(cg, t, nCompleted, 1 : n_n_tried(cg, t))
# end

# For a single no of courses tried: return Integer index, not vector
update_ncourses(cg :: CourseGrid, t :: Integer, nCompleted :: Integer, idxNTried :: Integer) =
    ncompleted_to_indices(cg, t+1, nCompleted + cg.nCourseV[idxNTried]);

# -----------------

ext: jl


text:
import NonUniformRandomVariateGeneration.sampleSortedUniforms!

# sample nans categorical r.v.s where Pr(X=i+offset) ∝ p[i]
# required: length(Fs) = length(p), unifs[start+1:start+nans] are valid
@inline function sampleSerial!(vs::Vector{Int64}, rng::RNG, p::V,
  nans::Int64, unifs::Vector{Float64}, Fs::V, start::Int64 = 0,
  offset::Int64 = 0) where V<:AbstractVector{Float64}
  # generateSortedUniformsBSerial!(unifs, start, nans, rng)
  sampleSortedUniforms!(unifs, start, nans, rng)
  cumsum!(Fs, p)
  @inbounds maxFs = Fs[length(p)]
  j = 1
  for i = 1:nans
    @inbounds while maxFs * unifs[start + i] > Fs[j]
      j += 1
    end
    @inbounds vs[start + i] = j + offset
  end
end


text:
include("helpers/listing.jl")

function main(window)
    push!(window.assets, "widgets")
    push!(window.assets, "codemirror")
    push!(window.assets, "tex")

    vbox(
        title(3, "Graphics"),

        vskip(1em),
        h2("Vector Graphics"),
        md"Any [Compose.jl](http://composejl.org/) graphic can be used inside Escher if you install the `ComposeDiff` package. ",
        vskip(1em),
        listing("""
        using Compose
        using ComposeDiff
        duck = compose(context(),
          (context(), circle(0.45cx, 0.15cy, 0.02cx), fill("white")),
          (context(), circle(0.55cx, 0.15cy, 0.02cx), fill("white")),
          (context(), ellipse(0.5cx, 0.25cy, 0.04cx, 0.02cx), fill("orange")),
          (context(), circle(0.5cx,0.2cy,0.12cx)),
          (context(), circle(0.5cx, 0.6cy, 0.2cx)),
          fill("#eeeeee")
          )"""
        ),
        md"To resize a Compose graphic, you can use the `drawing` function.",
        listing("""
        drawing(4Compose.inch, 3Compose.inch, duck)"""),

        vskip(1em),
        h2("Plots"),
        md"Gadfly plots are rendered in Escher.",
        vskip(1em),
        listing("""
        using Gadfly
        randomwalk = plot(x=[1:100],y=cumsum(randn(100)), Geom.line)
        """),

        md"If you want to render the graphic as a png file instead of using the default `Patchwork` renderer, you can specify the backend in `drawing`.",
        listing("""
        drawing(PNG(4Compose.inch, 3Compose.inch), randomwalk)"""),

    ) |> pad(2em)
end


text:
@testset "utils" begin

    # Test permutation matrix
    p = [5, 3, 4, 1, 2, 6]
    x = randn(length(p))
    @test x[p] == Chordal.build_perm_matrix(p)*x

    # Test random SDP
    n = 100
    c, F, G, xstar, D = Chordal.generate_random_sdp(n)
    Fx = sum(F[i]*xstar[i] for i in 1:n) + G
    @test (
        all(eigvals(Matrix(Fx)) .>= -sqrt(eps())) &&
        all(eigvals(D) .>= -sqrt(eps())) &&
        all([0 ≈ c[i] - tr(F[i]*D) for i in 1:n]) &&
        tr(Fx*D) ≈ 0
    )
end


text:
"""
    reconstruction_matrix(filter, domain, M, N)

Approximate `R` of size `N × M`, the inverse of the discrete filtering matrix `W` of size `M × N`.
"""
function reconstruction_matrix end

function reconstruction_matrix(::IdentityFilter, domain::AbstractIntervalDomain, M, N; kwargs...)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    interpolation_matrix(ξ, x)
end

function reconstruction_matrix(
    filter::TopHatFilter,
    domain::ClosedIntervalDomain,
    M,
    N;
    degmax = 10,
    λ = 0,
)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    h = filter.width
    R = spzeros(N + 1, M + 1)
    for n = 1:(N+1)
        # Reconstruction point
        ξₙ = ξ[n]

        # Find m such that ξₙ is reachable from xₘ
        dists = abs.(x .- ξₙ)
        m = dists .< h.(x)
        Mₙ = sum(m)
        xₘ = x[m]'
        hₘ = h.(xₘ)
        aₘ = max.(xₘ - hₘ, domain.left)
        bₘ = min.(xₘ + hₘ, domain.right)
        Δhₘ = bₘ - aₘ

        # Polynomial degree (Taylor series order)
        Q = min(degmax + 1, max(1, Mₙ == 2 ? 2 : 3 ≤ Mₙ < 9 ? 3 : floor(Int, √Mₙ)))
        q = 1:Q

        # Vandermonde matrix
        Zₙ = @. 1 / (Δhₘ * factorial(q)) * ((bₘ - ξₙ)^q - (aₘ - ξₙ)^q)

        # Moments
        vₙ = fill(0.0, Q)
        vₙ[1] = 1.0

        # Fit weights
        rₙ = ridge(Zₙ, vₙ)
        # rₙ = Zₙ \ vₙ

        # Store row of weights
        R[n, m] .= rₙ[:] ./ sum(rₙ)
    end
    dropzeros!(R)

    R
end

# function reconstruction_matrix(
#     filter::TopHatFilter,
#     domain::PeriodicIntervalDomain,
#     M,
#     N;
#     degmax = 10,
#     λ = 0,
# )
#     x = discretize(domain, M)
#     ξ = discretize(domain, N)
#     L = (domain.right - domain.left)
#     h = filter.width
#
#     # Get reconstruction weights for each point
#     R = spzeros(N, M)
#     for n = 1:N
#         # Reconstruction point
#         ξₙ = ξ[n]
#
#         # Move x by L * (shifts - 2) to get closer to ξₙ
#         dists = @. abs(ξₙ - x - [-L 0 L])
#         mininds = argmin(dists, dims = 2)
#         shifts = [mininds[m].I[2] for m in eachindex(mininds)]
#
#         # Find m such that ξₙ is reachable from xₘ
#         m = dists[mininds][:] .< h.(x)
#         Mₙ = sum(m)
#
#         # Vandermonde matrix
#         xₘ = x[m]'
#         hₘ = h.(xₘ)
#         sₘ = L * (shifts[m]' .- 2)
#
#         # Polynomial degree (Taylor series order)
#         Q = min(degmax + 1, max(1, Mₙ == 2 ? 2 : 3 ≤ Mₙ < 9 ? 3 : floor(Int, √Mₙ)))
#         q = 1:Q
#
#         # Vandermonde matrix
#         Zₙ = @. 1 / (2hₘ * factorial(q)) * ((xₘ + sₘ + hₘ - ξₙ)^q - (xₘ + sₘ - hₘ - ξₙ)^q)
#
#         # Moments
#         vₙ = fill(0.0, Q)
#         vₙ[1] = 1.0
#
#         # Fit weights
#         rₙ = ridge(Zₙ, vₙ)
#         # rₙ = Zₙ \ vₙ
#
#         # Store row of weights
#         R[n, m] .= rₙ[:] ./ sum(rₙ)
#     end
#     dropzeros!(R)
#
#     R
# end

function reconstruction_matrix(
    filter::TopHatFilter,
    domain::PeriodicIntervalDomain,
    M,
    N;
    degmax = 10,
    λ = 0,
)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    L = (domain.right - domain.left)
    h = filter.width
    R = spzeros(N, M)
    for n = 1:N
        # Reconstruction point
        ξₙ = ξ[n]

        # Move x by L * (shifts - 2) to get closer to ξₙ
        dists = @. abs(ξₙ - x - [-L 0 L])
        mininds = argmin(dists, dims = 2)
        shifts = [mininds[m].I[2] for m in eachindex(mininds)]

        # Find m such that ξₙ is reachable from xₘ and include one point outside
        m = dists[mininds][:] .< 1.5h.(x)
        if !any(m)
            # No coarse point is inside the interval. Set closest weight to one
            m = argmin(dists; dims = (1,2))[1].I[1]
            R[n, m] = 1
            continue
        end
        # m = mapreduce(s -> circshift(m, s), .|, [-1, 0, 1])
        hₘ = h.(x[m]')
        sₘ = L * (shifts[m]' .- 2)
        xₘ = x[m]' + sₘ
        Mₙ = length(xₘ)
        # bₘ = min.(xₘ + hₘ, domain.right + 0.001L)
        # Δhₘ = bₘ - aₘ
        ivals = xₘ .± hₘ
        # ivals = Interval.(aₘ, bₘ)
        # aₘ = max.(xₘ - hₘ, domain.left - 0.001L)

        # Polynomial basis
        Q = min(degmax + 1, Mₙ)
        # Q = min(degmax + 1, max(1, Mₙ == 2 ? 2 : 3 ≤ Mₙ < 9 ? 3 : floor(Int, √Mₙ)))
        ϕ = chebyshevt.(0:Q-1, [minimum(xₘ - 1.5hₘ)..maximum(xₘ + 1.5hₘ)])
        # ϕ = chebyshevt.(0:Q-1, [domain.left..domain.right])
        ϕ_int = integrate.(ϕ)
        # Build reconstruction system
        Zₙ = spzeros(Q, Mₙ)
        vₙ = zeros(Q)
        for q = 1:Q
            for (i, ival) ∈ enumerate(ivals)
                Zₙ[q, i] =
                    (ϕ_int[q](ival.right) - ϕ_int[q](ival.left)) / (ival.right - ival.left)
            end
            vₙ[q] = ϕ[q](ξₙ)
        end

        # Fit reconstruction weights
        rₙ = ridge(Zₙ, vₙ, λ)
        # rₙ = Zₙ \ vₙ
        # rₙ = fit(LassoRegression(λ; fit_intercept = false), Zₙ, vₙ)

        # Store row of weights
        R[n, m] .= rₙ[:] ./ sum(rₙ)
    end
    dropzeros!(R)

    R
end

function reconstruction_matrix(
    filter::ConvolutionalFilter,
    domain::ClosedIntervalDomain,
    M,
    N;
    degmax = 10,
    λ = 0,
)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    Δx = x[2] - x[1]
    Ival_domain = (domain.left - Δx / 1000)..(domain.right + Δx / 1000)
    h = filter.width
    G = filter.kernel
    R = spzeros(N + 1, M + 1)
    for n = 1:(N+1)
        # Reconstruction point
        ξₙ = ξ[n]

        # Find m such that ξₙ is reachable from xₘ
        dists = abs.(x .- ξₙ)
        m = dists .< h.(x)
        xₘ = x[m]
        hₘ = h.(xₘ)
        aₘ = max.(xₘ - hₘ, Ival_domain.left)
        bₘ = min.(xₘ + hₘ, Ival_domain.right)
        # Δhₘ = bₘ - aₘ
        # ivals = xₘ .± hₘ
        ivals = Interval.(aₘ, bₘ)

        # Polynomial basis
        Mₙ = sum(m)
        Q = min(degmax + 1, Mₙ)
        # Q = min(degmax + 1, max(1, Mₙ == 2 ? 2 : 3 ≤ Mₙ < 9 ? 3 : floor(Int, √Mₙ)))
        # ϕ = chebyshevt.(0:Q-1, [minimum(xₘ - hₘ)..maximum(xₘ + hₘ)])
        ϕ = chebyshevt.(0:Q-1, [minimum(aₘ)..maximum(bₘ)])
        # ϕ = chebyshevt.(0:Q-1, [Ival_domain])

        # Polynomials evaluated at ξₙ
        vₙ = [ϕ(ξₙ) for ϕ ∈ ϕ]

        # Moment matrix
        Zₙ = spzeros(Q, Mₙ)
        for mm = 1:Mₙ
            ival = ivals[mm]
            ϕₘₘ = Fun.(ϕ, [ival])
            kern = Fun(x -> G(x - xₘ[mm]), ival)
            Zₙ[:, mm] = sum.(kern .* ϕₘₘ) / sum(kern)
        end

        # Fit weights
        rₙ = ridge(Zₙ, vₙ, λ)
        # rₙ = Zₙ \ vₙ

        # Store row of weights
        R[n, m] .= rₙ[:] ./ sum(rₙ)
    end
    dropzeros!(R)

    R
end

function reconstruction_matrix(
    filter::ConvolutionalFilter,
    domain::PeriodicIntervalDomain,
    M,
    N;
    degmax = 10,
    λ = 0,
)
    x = discretize(domain, M)
    ξ = discretize(domain, N)
    Δx = x[2] - x[1]
    L = (domain.right - domain.left)
    h = filter.width
    G = filter.kernel
    R = spzeros(N, M)
    for n = 1:N
        # Reconstruction point
        ξₙ = ξ[n]

        # Move x by L * (shifts - 2) to get closer to ξₙ
        dists = @. abs(ξₙ - x - [-L 0 L])
        mininds = argmin(dists, dims = 2)
        shifts = [mininds[m].I[2] for m in eachindex(mininds)]

        # Find m such that ξₙ is reachable from xₘ
        m = dists[mininds][:] .< h.(x)
        sₘ = L * (shifts[m]' .- 2)
        xₘ = x[m]' + sₘ
        hₘ = h.(xₘ)
        aₘ = max.(xₘ - hₘ, domain.left - Δx / 1000)
        bₘ = min.(xₘ + hₘ, domain.right + Δx / 1000)
        # Δhₘ = bₘ - aₘ
        ivals = xₘ .± hₘ
        # ivals = Interval.(aₘ, bₘ)
        Mₙ = sum(m)

        # Polynomial degree
        Q = min(degmax + 1, Mₙ)
        # Q = min(degmax + 1, max(1, Mₙ == 2 ? 2 : 3 ≤ Mₙ < 9 ? 3 : floor(Int, √Mₙ)))
        ϕ = chebyshevt.(0:Q-1, [minimum(xₘ - hₘ)..maximum(xₘ + hₘ)])
        # ϕ = chebyshevt.(0:Q-1, [domain.left..domain.right])

        # Polynomials evaluated at ξₙ
        vₙ = [ϕ(ξₙ) for ϕ ∈ ϕ]

        # Polynomial moments
        Zₙ = spzeros(Q, Mₙ)
        for mm = 1:Mₙ
            ival = ivals[mm]
            ϕₘₘ = Fun.(ϕ, [ival])
            kern = Fun(x -> G(x - xₘ[mm]), ival)
            kern /= sum(kern)
            Zₙ[:, mm] = sum.(kern .* ϕₘₘ)
        end

        # Fit weights
        rₙ = ridge(Zₙ, vₙ, λ)
        # rₙ = Zₙ \ vₙ

        # Store row of weights
        R[n, m] .= rₙ[:] ./ sum(rₙ)
    end
    dropzeros!(R)

    R
end


text:
#******************************************************************************#
#                                                                              #
#                                                         :::      ::::::::    #
#    shell_colors.jl                                    :+:      :+:    :+:    #
#                                                     +:+ +:+         +:+      #
#    By: mngomane <mngomane@student.42.fr>          +#+  +:+       +#+         #
#                                                 +#+#+#+#+#+   +#+            #
#    Created: 2015/02/12 07:30:42 by mngomane          #+#    #+#              #
#    Updated: 2015/02/12 07:30:42 by mngomane         ###   ########.fr        #
#                                                                              #
#******************************************************************************#

C_PURPLE = "\033[35;1m"
C_RED = "\033[31;1m"
C_GREEN = "\033[32;1m"
C_YELLOW = "\033[33;1m"
C_END = "\033[0m"

max_stars_repo_name: mngomane/projects

max_issues_repo_name: mngomane/projects

max_forks_repo_name: mngomane/projects


text:
function test_content(f::Base.Callable, fname)
  try
    open(fname) do io
      f(read(io, String))
    end
  finally
    rm(fname, force = true)
  end
end

@testset "Send" begin
  logfile = tempname()
  server = "smtp://127.0.0.1:1025"
  addr = "<julian@julialang.org>"
  mock = joinpath(dirname(@__FILE__), "mock.py")

  cmd = `python3.7 $mock $logfile`
  smtpsink = run(pipeline(cmd, stderr=stdout), wait = false)
  sleep(.5)  # wait for fake smtp server ready

  try
    let  # send with body::IOBuffer
      body = IOBuffer("body::IOBuffer test")
      send(server, [addr], addr, body)
      test_content(logfile) do s
        @test occursin("body::IOBuffer test", s)
      end
    end

    let  # send with body::IOStream
      mktemp() do path, io
        write(io, "body::IOStream test")
        seekstart(io)

        send(server, [addr], addr, io)
        test_content(logfile) do s
          @test occursin("body::IOStream test", s)
        end
      end
    end

    let  # AUTH PLAIN
      opts = SendOptions(username = "foo@example.org", passwd = "bar")
      body = IOBuffer("AUTH PLAIN test")
      send(server, [addr], addr, body, opts)
      test_content(logfile) do s
        @test occursin("AUTH PLAIN test", s)
      end
    end

    let
      opts = SendOptions(username = "foo@example.org", passwd = "invalid")
      body = IOBuffer("invalid password")
      @test_throws Exception send(server, [addr], addr, body, opts)
    end

    let  # multiple RCPT TO
      body = IOBuffer("multiple rcpt")
      rcpts = ["<foo@example.org>", "<bar@example.org>", "<baz@example.org>"]
      send(server, rcpts, addr, body)

      test_content(logfile) do s
        @test occursin("multiple rcpt", s)
        @test occursin("X-RCPT: foo@example.org", s)
        @test occursin("X-RCPT: bar@example.org", s)
        @test occursin("X-RCPT: baz@example.org", s)
      end
    end

    let  # non-blocking send
      body = IOBuffer("non-blocking send")
      task = @async send(server, [addr], addr, body)
      wait(task)
      test_content(logfile) do s
        @test occursin("non-blocking send", s)
      end
    end

    let  # SendOptions.verbose no error
      opts = SendOptions(verbose = true, username = "foo@example.org", passwd = "bar")
      body = IOBuffer("SendOptions.verbose")
      send(server, [addr], addr, body, opts)
      test_content(logfile) do s
        @test occursin("SendOptions.verbose", s)
      end
    end

  finally
    kill(smtpsink)
    rm(logfile, force = true)
  end
end  # @testset "Send"


text:
module Utils

"""
    has_header(headers::Dict{String, String}, key::String, value::String) :: Bool

Check whether corresponding (key, value) pair exists in Header
"""
function has_header(headers::Dict{String, String}, key::String, value::String) :: Bool
    return headers[key] == value
end

end


text:
using Bildefin
using Test

@testset "Bildefin.jl" begin
    # Write your tests here.
end


text:
import Base: copy, copyto!, ==
import Random: rand, seed!
import RandomNumbers: AbstractRNG, gen_seed, split_uint, seed_type

"""
```julia
AbstractXorshift128 <: AbstractRNG{UInt64}
```

The base abstract type for `Xorshift128`, `Xorshift128Star` and `Xorshift128Plus`.
"""
abstract type AbstractXorshift128 <: AbstractRNG{UInt64} end

for (star, plus) in (
        (false, false),
        (false, true),
        (true, false),
    )
    rng_name = Symbol(string("Xorshift128", star ? "Star" : plus ? "Plus" : ""))
    @eval begin
        mutable struct $rng_name <: AbstractXorshift128
            x::UInt64
            y::UInt64
        end

        function $rng_name(seed::NTuple{2, UInt64}=gen_seed(UInt64, 2))
            r = $rng_name(zero(UInt64), zero(UInt64))
            seed!(r, seed)
            r
        end

        $rng_name(seed::Integer) = $rng_name(init_seed(seed, UInt64, 2))

        @inline function xorshift_next(r::$rng_name)
            t = r.x ⊻ r.x << 23
            r.x = r.y
            r.y = t ⊻ (t >> 3) ⊻ r.y ⊻ (r.y >> 24)
            $(star ? :(r.y * 2685821657736338717) :
              plus ? :(r.y + r.x) : :(r.y))
        end
    end
end

@inline seed_type(::Type{T}) where T <: AbstractXorshift128 = NTuple{2, UInt64}

function copyto!(dest::T, src::T) where T <: AbstractXorshift128
    dest.x = src.x
    dest.y = src.y
    dest
end

copy(src::T) where T <: AbstractXorshift128 = copyto!(T(), src)

==(r1::T, r2::T) where T <: AbstractXorshift128 = r1.x == r2.x && r1.y == r2.y

seed!(r::AbstractXorshift128, seed::Integer) = seed!(r, split_uint(seed % UInt128))

function seed!(r::AbstractXorshift128, seed::NTuple{2, UInt64}=gen_seed(UInt64, 2))
    all(==(0), seed) && error("0 cannot be the seed")
    r.x = seed[1]
    r.y = seed[2]
    xorshift_next(r)
    xorshift_next(r)
    r
end

@inline rand(r::AbstractXorshift128, ::Type{UInt64}) = xorshift_next(r)


text:
# Import modules
using Distributions  # probability distributions and associated functions
using StatsBase
using LinearAlgebra
# Import Calibrate-Emulate-Sample modules
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributionStorage
include(joinpath(@__DIR__, "helper_funcs.jl"))
using NPZ

###
###  Define the parameters and their priors
###

# Define the parameters that we want to learn
param_names = ["entrainment_factor", "detrainment_factor", "sorting_power",
	"tke_ed_coeff", "tke_diss_coeff", "pressure_normalmode_adv_coeff",
         "pressure_normalmode_buoy_coeff1", "pressure_normalmode_drag_coeff", "static_stab_coeff"]
n_param = length(param_names)

###
###  Retrieve true LES samples from PyCLES data
###

# This is the true value of the observables (e.g. LES ensemble mean for EDMF)
ti = [10800.0, 28800.0, 10800.0, 18000.0]
tf = [14400.0, 32400.0, 14400.0, 21600.0]
y_names = Array{String, 1}[]
push!(y_names, ["thetal_mean", "ql_mean", "qt_mean", "total_flux_h", "total_flux_qt"]) #DYCOMS_RF01
push!(y_names, ["thetal_mean", "u_mean", "v_mean", "tke_mean"]) #GABLS
push!(y_names, ["thetal_mean", "total_flux_h"]) #Nieuwstadt
push!(y_names, ["thetal_mean", "ql_mean", "qt_mean", "total_flux_h", "total_flux_qt"]) #Bomex

# Get observations
sim_names = ["DYCOMS_RF01", "GABLS", "Nieuwstadt", "Bomex"]
for i in range(0,1,length=2)
    normalize = Bool(i)
    suffix = normalize ? "_normalized" : "_wscales"

    yt = zeros(0)
    yt_var_list = []
    yt_pca_list = zeros(0)
    yt_pca_var_list = []
    P_pca_list = []

    les_dir = string("/groups/esm/ilopezgo/Output.", sim_names[1],".may20")
    sim_dir = string("Output.", sim_names[1],".00000")
    z_scm = get_profile(sim_dir, ["z_half"])
    yt_, yt_var_ = obs_LES(y_names[1], les_dir, ti[1], tf[1], z_scm = z_scm, normalize=normalize)
    append!(yt, yt_)
    push!(yt_var_list, yt_var_)
    npzwrite("dycoms_z.npy", z_scm)
    yt_pca, yt_var_pca, P_pca = obs_PCA(yt_, yt_var_)
    append!(yt_pca_list, yt_pca)
    push!(yt_pca_var_list, yt_var_pca)
    push!(P_pca_list, P_pca)
    @assert length(yt_pca) == length(yt_var_pca[1,:])

    les_dir = string("/groups/esm/ilopezgo/Output.", sim_names[2],".iles128wCov")
    sim_dir = string("Output.", sim_names[2],".00000")
    z_scm = get_profile(sim_dir, ["z_half"])
    yt_, yt_var_ = obs_LES(y_names[2], les_dir, ti[2], tf[2], z_scm = z_scm, normalize=normalize)
    append!(yt, yt_)
    push!(yt_var_list, yt_var_)
    npzwrite("gabls_z.npy", z_scm)
    yt_pca, yt_var_pca, P_pca = obs_PCA(yt_, yt_var_)
    append!(yt_pca_list, yt_pca)
    push!(yt_pca_var_list, yt_var_pca)
    push!(P_pca_list, P_pca)
    @assert length(yt_pca) == length(yt_var_pca[1,:])

    les_dir = string("/groups/esm/ilopezgo/Output.Soares.dry11")
    sim_dir = string("Output.", sim_names[3],".00000")
    z_scm = get_profile(sim_dir, ["z_half"])
    yt_, yt_var_ = obs_LES(y_names[3], les_dir, ti[3], tf[3], z_scm = z_scm, normalize=normalize)
    append!(yt, yt_)
    push!(yt_var_list, yt_var_)
    npzwrite("nieuwstadt_z.npy", z_scm)
    yt_pca, yt_var_pca, P_pca = obs_PCA(yt_, yt_var_)
    append!(yt_pca_list, yt_pca)
    push!(yt_pca_var_list, yt_var_pca)
    push!(P_pca_list, P_pca)
    @assert length(yt_pca) == length(yt_var_pca[1,:])

    les_dir = string("/groups/esm/ilopezgo/Output.Bomex.may18")
    sim_dir = string("Output.", sim_names[4],".00000")
    z_scm = get_profile(sim_dir, ["z_half"])
    yt_, yt_var_ = obs_LES(y_names[4], les_dir, ti[4], tf[4], z_scm = z_scm, normalize=normalize)
    append!(yt, yt_)
    push!(yt_var_list, yt_var_)
    npzwrite("bomex_z.npy", z_scm)
    yt_pca, yt_var_pca, P_pca = obs_PCA(yt_, yt_var_)
    append!(yt_pca_list, yt_pca)
    push!(yt_pca_var_list, yt_var_pca)
    push!(P_pca_list, P_pca)
    @assert length(yt_pca) == length(yt_var_pca[1,:])

    yt_var = zeros(length(yt), length(yt))
    vars_num = 1
    for sim_covmat in yt_var_list
        vars = length(sim_covmat[1,:])
        yt_var[vars_num:vars_num+vars-1, vars_num:vars_num+vars-1] = sim_covmat
        vars_num = vars_num+vars
    end

    yt_var_pca = zeros(length(yt_pca_list), length(yt_pca_list))
    vars_num = 1
    for sim_covmat in yt_pca_var_list
        vars = length(sim_covmat[1,:])
        yt_var_pca[vars_num:vars_num+vars-1, vars_num:vars_num+vars-1] = sim_covmat
        vars_num = vars_num+vars
    end

    # Write covariances to file
    npzwrite("obs_cov"*suffix*".npy", yt_var)
    npzwrite("obs_cov_pca"*suffix*".npy", yt_var_pca)
end


text:
delete_missing(x) = [collect(skipmissing(x[:,i])) for i in 1:ncol(x)]

function bin_dissociations(data, n_bins; bin_type = "force", edges=nothing)
    F = data[1]
    t = data[2]
    # even force bins
    if edges != nothing
        bins,bin_ids = partition(F,edges)
        binned_t = [similar(t,0) for i in bins]

    elseif bin_type == "force"
        f_min = minimum(F)
        f_max = maximum(F)

        edges = range(f_min,stop=f_max,length=n_bins+1)
        bins,bin_ids = partition(F,edges)
        binned_t = [similar(t,0) for i in bins]

    elseif bin_type == "equal number"
        edges_ind = Int64.(round.(range(1,stop=length(F),length=n_bins+1)))
        edges = F[edges_ind]
        # edges = [0,9,12,15,19,30]
        bins,bin_ids = partition(F,edges)
        binned_t = [similar(t,0) for i in bins]
    end
    # h = fit(Histogram, F)

    for (i,bin) in enumerate(bin_ids)
        push!(binned_t[bin],t[i])
    end

    return bins, binned_t, bin_ids
end

function partition(data, intervals)
    ranges = intervals[1:end-1] .=> intervals[2:end]
    bins = [similar(data, 0) for _ in 1:length(ranges)]
    ids = similar(data,0)
    for x in data
        for (i, (a, b)) in pairs(ranges)
            if i != length(bins) && a <= x < b
                push!(bins[i], x)
                push!(ids,i)
                break

            elseif i == length(bins) && a <= x <= b
                push!(bins[i], x)
                push!(ids,i)
                break
        end
        end
    end
    return bins, Int64.(ids)
end

function resort(f,t)
    t_sorted = sort.(t)
    sort_id = sortperm.(t)
    f_sorted = [f[i][sort_id[i]] for i in 1:length(sort_id)]
    return f_sorted, t_sorted
end

survival_p(t) = collect((length(t) :-1:1) / length(t))

p2lnp(p) = map((x) -> log.(x), p)


text:
push!(LOAD_PATH, "$(homedir())/.julia/dev")

using Test
using SafeTestsets# , TestSetExtensions

@testset verbose = true "SymbolicCRN.jl" begin
    @safetestset "Cones.jl" begin include("Cones.jl") end;

    @safetestset "ConvexParameters.jl" begin include("ConvexParameters.jl") end;

    @safetestset "FindRoots.jl" begin end;

    @safetestset "HomotopyContinuation.jl" begin end;

    @safetestset "NemoUtils.jl" begin include("NemoUtils.jl") end;

    @safetestset "PosNegVertices.jl" begin end;

    @safetestset "Stability.jl" begin end;

    @safetestset "Utils.jl" begin end;

end


text:
using Test
using WorldBankData

# checks the iso 2-letter country codes for updates
@test WorldBankData.download_countries(verbose=true)[!, :iso2c] == WorldBankData.all_countries


text:
#Code written by Skipper af Rosenborg and Esben Bager
using StatsBase
using DataFrames
using CSV

#Load training size
trainingSizeInput = parse(Int64, ARGS[1])
#trainingSizeInput = 12
println(typeof(trainingSizeInput))

#set path
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/ML/Lasso_Test"
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/ML"
cd(path)
@everywhere include("ParallelModelGeneration.jl")
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Function to run LASSO and output result files
function runLassos(VIX, raw, expTrans, timeTrans, TA, trainingSize)
	#Set path and load data

	#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
	path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"

	#= industry must be one of the following
	NoDur
	Durbl
	Manuf
	Enrgy
	HiTec
	Telcm
	Shops
	Hlth
	Utils
	Other
	=#
	industry = ARGS[2]
	#industry = "NoDur"
	mainData = loadIndexDataLOGReturn(industry, path)
	mainDataArr = Array(mainData[:,1:end-2])

	path = "/zhome/9f/d/88706/SpecialeCode/"
	#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/"
	dateAndReseccion = Array(mainData[:,end-1:end])
	mainDataArr = Array(mainData[:,1:end-2])

	#Extract names
	colNames = names(mainData)
	if VIX == 1
		if timeTrans == 1
			featureNames = expandColNamesTimeToString(colNames[1:end-3])
		else
			featureNames = colNames[1:end-3]
		end
	else
		if timeTrans == 1
			featureNames = expandColNamesTimeToString(colNames[1:end-4])
		else
			featureNames = colNames[1:end-4]
		end
	end
	s = split(expandedColNamesToString(featureNames,expTrans, TA),",")
	s = Array{String}(s)

	nRows = size(mainDataArr)[1]
	nCols = size(mainDataArr)[2]

	#Pipeline data and include different transformations
	fileName = path*"Results/IndexData/LassoTest/"*industry*"/"*string(trainingSize)*"-1/"*string(trainingSize)*"_"

	mainYarr = mainDataArr[:, nCols:nCols]

	mainXarr = mainDataArr[:,1:nCols-1]
	if raw == 1
		fileName = fileName*"Raw"
		mainXarr = mainXarr[:,1:10]
	else
		if VIX == 1
			fileName = fileName*"VIX_Macro"
		else
			fileName = fileName*"Macro"
			mainXarr = mainXarr[:,1:end-1]
		end
	end

	# Transform with time elements
	if timeTrans == 1
		fileName = fileName*"Time"
		mainXarr = expandWithTime3612(mainXarr)
	end

	# Transform #
	if expTrans == 1
		fileName = fileName*"Exp"
		mainXarr = expandWithTransformations(mainXarr)
	end

	if TA == 1
		fileName = fileName*"TA"
		mainXarr = expandWithMAandMomentum(mainXarr, mainYarr, (nCols-1))
	end

	# Standardize
	standX = zScoreByColumn(mainXarr)
	standY = mainYarr

	allData = hcat(standX, standY)
	bCols = size(standX)[2]
	nRows = size(standX)[1]
	println(" \nSolving Convex\n")

	fileName = fileName*"_"

	println(fileName)

	#Set gamma and initilise arrays
	nGammas = 10
	predictions = 1
	testRuns = nRows-trainingSize-predictions

	ISR = SharedArray{Float64}(testRuns, nGammas)
	Indi = SharedArray{Float64}(testRuns, nGammas)
	#gammaArray = log.(logspace(0, SSTO/2, nGammas))
	predictedArr = SharedArray{Float64}(testRuns, nGammas)
	realArr = SharedArray{Float64}(testRuns, 1)
	bSolvedArr = SharedArray{Float64}(testRuns, bCols)
	bSolMatrix = SharedArray{Float64}(testRuns, bCols, nGammas)

	xTrainInput = [allData[r:(trainingSize+r-1), 1:bCols] for r = 1:testRuns]
	YtrainInput = [allData[r:(trainingSize+r-1), bCols+1] for r = 1:testRuns]
	XpredInput  = [allData[(trainingSize+r+1-1), 1:bCols] for r = 1:testRuns]
	YpredInput  = [allData[(trainingSize+r+1-1), bCols+1] for r = 1:testRuns]

	gammaArray = zeros(testRuns,nGammas)
	for r = 1:testRuns
		maxVal = findmax(abs.(xTrainInput[r]'*YtrainInput[r]))[1]
		gammaArray[r,:] = Array(linspace(0,maxVal,nGammas))
 	end

	for i=1:nGammas
		gammaInput = gammaArray[:,i]

		#Loop over gamma and solve LASSO in parallel
		results = pmap(generatSolveAndProcess, xTrainInput, YtrainInput, XpredInput, YpredInput, gammaInput)
		results
		#Log results
		for j in 1:testRuns
			ISR[j,i] = results[j][1]
			Indi[j,i] = results[j][2]
			predictedArr[j,i] = results[j][3]
			realArr[j,1] = YpredInput[j][1]
			bSolMatrix[j,:,i] = results[j][4]
		end
	end

	#Output bmatrix files
	for i=1:nGammas
		writedlm(fileName*"bMatrix"*string(trainingSizeInput)*"_"*string(i)*".csv",bSolMatrix[:,:,i],",")
	end

	#Output files
	dateAndReseccionOutput = dateAndReseccion[trainingSize+1:trainingSize+testRuns,:]
	toInput = Array{String}(1,2)
	toInput[1,1] = "Date"
	toInput[1,2] = "Recession"
	dateAndReseccionOutput = vcat(toInput,dateAndReseccionOutput)
	runCounter = collect(0:testRuns)

	combinedArray = vcat((round.(gammaArray[1,:],3))', ISR)
	ISRcomb = hcat(runCounter,dateAndReseccionOutput, combinedArray)
	writedlm(fileName*"ISRsquared.CSV", ISRcomb,",")

	combinedArray = vcat((round.(gammaArray[1,:],3))', Indi)
	Indicomb = hcat(runCounter, dateAndReseccionOutput, combinedArray)
	writedlm(fileName*"Indi.CSV", Indicomb,",")

	combinedArray = vcat("Real value", realArr)
	realArr = hcat(runCounter, dateAndReseccionOutput, combinedArray)
	writedlm(fileName*"real.CSV", realArr,",")

	combinedArray = vcat((round.(gammaArray[1,:],3))', predictedArr)
	predictedArr = hcat(runCounter, dateAndReseccionOutput, combinedArray)
	writedlm(fileName*"predicted.CSV", predictedArr,",")

	#writedlm(fileName*"bSolved.CSV", bSolvedArr,",")
end

#Loops to process all dataset transformations
#Testing raw dataset
for timeTrans = 0:1
	for expTrans = 0:1
		for TA = 0:1
			runLassos(0, 1, timeTrans, expTrans, TA, trainingSizeInput)
		end
	end
end

#Testing macro dataset, including and excluding VIX
for VIX = 0:1
	for timeTrans = 0:1
		for expTrans = 0:1
			for TA = 0:1
				runLassos(VIX, 0, timeTrans, expTrans, TA, trainingSizeInput)
			end
		end
	end
end

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/ML"
#cd(path)
println("Finished")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
#! /usr/bin/julia

# Rosetta Code, URL encoding

function sortselected{S<:Real,T<:Integer}(a::AbstractArray{S,1},
                                          s::AbstractArray{T,1})
    sel = unique(sort(s))
    if sel[1] < 1 || length(a) < sel[end]
        throw(ArgumentError("Tried to select outside of input array."))
    end
    b = collect(copy(a))
    b[sel] = sort(b[sel])
    return b
end

a = [7, 6, 5, 4, 3, 2, 1, 0]
sel = [7, 2, 8]
b = sortselected(a, sel)

print("Original Array:  ", a)
println(" sorted on  ", sel)
println("Sorted Array:    ", b)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# Meta-Learning
# ==============
#
# This is a shortened version of https://www.domluna.me/meta-learning/.
#
# > Meta-learning is an algorithm which learns to learn.
#
# We will use meta-learning to train a neural net which
# can be fintuned to accurately models a sine wave after
# training for few iterations on a handful of x, y coordinates.
#
# References:
#
# [Model-Agnostic Meta-Learning for Fast Adapdation of Deep Networks](https://arxiv.org/abs/1703.03400).
# [Reptile: A Scalable Meta-Learning Algorithm](https://blog.openai.com/reptile/)

using Flux
using Printf
using Plots
using Distributions: Uniform, Normal
using Statistics: mean
using Base.Iterators: partition
using Random: randperm, seed!

include("utils.jl")
include("linear.jl")

# Sine Waves
# -----------
#
# We'll be training the neural network on sine waves in the
# x coordinate of range $[-5, 5]$. The sine wave can be shifted left
# or right by the `phase_shift` and shrunk or streched by the `amplitude`.

struct SineWave
    amplitude::Float32
    phase_shift::Float32
end
SineWave() = SineWave(rand(Uniform(0.1, 5)), rand(Uniform(0, 2pi)))

(s::SineWave)(x::AbstractArray) = s.amplitude .* sin.(x .+ s.phase_shift)

function Base.show(io::IO, s::SineWave)
    print(io, "SineWave(amplitude = ", s.amplitude, ", phase shift = ", s.phase_shift, ")")
end

# Now let's make some waves ...

x = LinRange(-5, 5, 100)
#-
wave = SineWave(4, 1)
#-
plot(x, wave(x))

# Nice, time to learn to learn!

# `tanh` activations seem to produce smoother lines than `relu`.

# First-Order MAML (FOMAML)
# --------------------------
#
# MAML can be succinctly described as:
#
# $$
# \phi_i = \theta - \alpha \bigtriangledown_{\theta} Loss(D_{t_i}^{train}; \theta) \\
# \theta \leftarrow \theta - \beta \sum_i \bigtriangledown_{\theta} Loss(D_{t_i}^{test}; \phi_i)
# $$
#
# The first line shows the inner update, which optimizes $\theta$ towards a solution for the task
# training set $D_{t_i}^{train}$ producing $\phi_i$. The following line is the meta update which
# aims to generalize to new data. This involves evaluating all $\phi_i$ on the test sets $D_{t_i}^{test}$
# and accumulating the resulting gradients. $\alpha$ and $\beta$ are learning rates.
#
# The difference between MAML and FOMAML (first-order MAML) is the inner gradient, shown is red
# is ignored during backpropagation:
#
# $$
# \theta = \theta - \beta \sum_i \bigtriangledown_{\theta} Loss( \theta - \alpha {\color{red} \bigtriangledown_{\theta} Loss(\theta, D_{t_i}^{train})}, D_{t_i}^{test})
# $$
#
# We'll focus on FOMAML since it's less computationally expensive and achieves similar performance
# to MAML in practice and has a simpler implementation.
#
# FOMAML generalizes by further adjusting parameters based on how performance on a validation
# or test set (used interchangeably) during the meta update. This is best shown visually:
#
# ![fomaml gradient updates](./fomaml_grad.png)
#
# $\theta$ is the starting point, $\phi_{3}$ is the parameter values after 3 gradient updates on a task.
# Notice before the meta update the parameters shift in the direction of the new task's solution (red arrow)
# but after the meta update they change direction (blue arrow). This illustrates how the meta update
# adjusts the gradient for task generalization.

function fomaml(model; meta_opt=Descent(0.01), inner_opt=Descent(0.02), epochs=30_000,
              n_tasks=3, train_batch_size=10, eval_batch_size=10, eval_interval=1000)

    weights = params(model)
    dist = Uniform(-5, 5)
    testx = Float32.(range(-5, stop=5, length=50))

    for i in 1:epochs
        prev_weights = deepcopy(Flux.data.(weights))

        for _ in 1:n_tasks
            task = SineWave()

            x = Float32.(rand(dist, train_batch_size))
            y = task(x)
            grad = Flux.Tracker.gradient(() -> Flux.mse(model(x'), y'), weights)

            for w in weights
                w.data .-= Flux.Optimise.apply!(inner_opt, w.data, grad[w].data)
            end

            testy = task(testx)
            grad = Flux.Tracker.gradient(() -> Flux.mse(model(testx'), testy'), weights)

            # reset weights and accumulate gradients
            for (w1, w2) in zip(weights, prev_weights)
                w1.data .= w2
                w1.grad .+= grad[w1].data
            end

        end

        Flux.Optimise._update_params!(meta_opt, weights)

        if i % eval_interval == 0
            @printf("Iteration %d, evaluating model on random task...\n", i)
            evalx = Float32.(rand(dist, eval_batch_size))
            eval_model(model, evalx, testx, SineWave())
        end

    end
end

# Training a neural net using FOMAML.

seed!(0)
fomaml_model = Chain(Linear(1, 64, tanh), Linear(64, 64, tanh), Linear(64, 1))
fomaml(fomaml_model, meta_opt=Descent(0.01), inner_opt=Descent(0.02), epochs=50_000, n_tasks=3)

# Evaluate the model on 10 (x, y) sampled uniformly from $[-5, 5]$

x = rand(Uniform(-5, 5), 10)
testx = LinRange(-5, 5, 50)
#-
fomaml_data = eval_model(fomaml_model, x, testx, wave, updates=32, opt=Descent(0.02))
#-
p = plot_eval_data(fomaml_data, "FOMAML - SGD Optimizer")
#-
plot(p)

# Reptile
# --------------------------
#
# Reptile optimizes to find a point (parameter representation) in manifold space which is
# closest in euclidean distance to a point in each task's manifold of optimal solutions.
# To achieve this we minimize the expected value for all tasks $t$ of $(\theta - \phi_{*}^{t})^2$
# where $\theta$ is the model's parameters and $\phi_{*}^{t}$ are the optimal parameters for task $t$.
#
# $$
# E_{t}[\frac{1}{2}(\theta - \phi_{*}^{t})^2]
# $$
#
# In each iteration of Reptile we sample a task and update $\theta$ using SGD:
#
# $$
# \theta \leftarrow \theta - \alpha \bigtriangledown_{\theta} \frac{1}{2}(\theta - \phi_{*}^{t})^2 \\
# \theta \leftarrow \theta - \alpha(\theta - \phi_{*}^{t})
# $$
#
# In practice an approximation of $\phi_{*}^{t}$ is used since it's not feasible to compute. The approximation
# is $\phi$ after $i$ gradient steps $\phi_{i}^{t}$.
#
# ![Reptile gradient updates](./reptile_grad.png)
#
# This a Reptile update after training for 3 gradient steps on task data. Note with Reptile there's
# no train and test data, just data. The direction of the gradient update (blue arrow) is directly
# in the direction towards $\phi_i$.  It's kind of crazy that this actually works. Section 5 of
# the [Reptile paper](https://arxiv.org/abs/1803.02999) has an analysis showing the gradients of
# MAML, FOMMAL and Reptile are similar within constants.

function reptile(model; meta_opt=Descent(0.1), inner_opt=Descent(0.02), epochs=30_000,
                 train_batch_size=10, eval_batch_size=10, eval_interval=1000)
    weights = params(model)
    dist = Uniform(-5, 5)
    testx = Float32.(range(-5, stop=5, length=50))
    x = testx

    for i in 1:epochs
        prev_weights = deepcopy(Flux.data.(weights))
        task = SineWave()

        # Train on task for k steps on the dataset
        y = task(x)
        for idx in partition(randperm(length(x)), train_batch_size)
            l = Flux.mse(model(x[idx]'), y[idx]')
            Flux.back!(l)
            Flux.Optimise._update_params!(inner_opt, weights)
        end

        # Reptile update
        for (w1, w2) in zip(weights, prev_weights)
            gw = Flux.Optimise.apply!(meta_opt, w2, w1.data - w2)
            @. w1.data = w2 + gw
        end

        if i % eval_interval == 0
            @printf("Iteration %d, evaluating model on random task...\n", i)
            evalx = Float32.(rand(dist, eval_batch_size))
            eval_model(model, evalx, testx, SineWave())
        end

    end
end

# Training a neural net using Reptile.

seed!(0)
reptile_model = Chain(Linear(1, 64, tanh), Linear(64, 64, tanh), Linear(64, 1))
reptile(reptile_model, meta_opt=Descent(0.1), inner_opt=Descent(0.02), epochs=50_000)

# Evaluate the model on 10 x coordinates sampled uniformly from $[-5, 5]$

reptile_data = eval_model(reptile_model, x, testx, wave, updates=32, opt=Descent(0.02))
#-
p = plot_eval_data(reptile_data, "Reptile - SGD Optimizer")
#-
plot(p)

# Testing Robustness
# -------------------

# To test if the FOMAML and Reptile representations *learned to learn quickly with minimal data* we'll
# finetune on 5 datapoints for 10 update steps. The x values are sampled from a uniform distribution
# of $[0, 5]$, the right half of the sine wave. Can the entire wave be learned?

x = rand(Uniform(0, 5), 5)

# First FOMAML:

fomaml_data = eval_model(fomaml_model, x, testx, wave, updates=10, opt=Descent(0.02))
#-
p = plot_eval_data(fomaml_data, "FOMAML - 5 samples, 10 updates")
#-
plot(p)

# Now Reptile:

reptile_data = eval_model(reptile_model, x, testx, wave, updates=10, opt=Descent(0.02))
#-
p = plot_eval_data(reptile_data, "Reptile - 5 samples, 10 updates")
#-
plot(p)

# Lastly, let's see how quickly a representation is learned for a new task. We'll use the same
# 5 element sample as before and train for 10 gradient updates.

plot([fomaml_data.test_losses, reptile_data.test_losses],
     label=["FOMAML", "Reptile"],
     xlabel="Updates", ylabel="Loss")

# Cool! Both FOMAML and Reptile learn a very useful representation within a few updates.

max_stars_repo_name: c-p-murphy/model-zoo

max_issues_repo_name: c-p-murphy/model-zoo

max_forks_repo_name: c-p-murphy/model-zoo


text:
"""
    Pendulum{R}

The most simple canonical nonlinear system, consisting of a simple pendulum controlled by
torque at the base. Has 2 state and 1 control.

# Constructor
    Pendulum(; kwargs...)

with keyword arguments
* `mass` - mass of the pendulum, in kg (default = 1.0)
* `length` - length of the pendulum, in m (default = 0.5)
* `b` - friction coefficient (default=0.1)
* `I` - inertia, in kg⋅m² (default=0.25)
* `g` - gravity, in kg/m² (default=9.81)
"""
@with_kw mutable struct Pendulum{T} <: AbstractModel
    mass::T = 1.
    length::T = 0.5
    b::T = 0.1
    lc::T = 0.5
    I::T = 0.25
    g::T = 9.81
end

RobotDynamics.state_dim(::Pendulum) = 2
RobotDynamics.control_dim(::Pendulum) = 1

function dynamics(p::Pendulum, x, u)
    m = p.mass * p.lc * p.lc
    @SVector [x[2],
              u[1]/m - p.g*sin(x[1])/p.lc - p.b*x[2]/m]
end
@inline Base.position(::Pendulum, x) = @SVector zeros(3)
orientation(::Pendulum, x) = expm((pi-x[1])*@SVector [1,0,0.])


text:
export entropyM

"""
function entropyM(param,w;pwr)

Compute the entropy of the probability p(M|w,d) for given w

Input:

	1. param - StrainReconParam
	2. w     - frequency vector

Keyword argument:

    3. pwr   - exponent in the likelihood function (default = param.pwr)

Output:

	1. E     - entropy
"""

function entropyM(param::StrainReconParam,w::Vector;pwr::Float64=param.pwr)
	n  = length(w)
	d  = param.d
	Rp = diag(param.R.^pwr)
	p  = param.p
	m  = param.m

	B  = getB(param)
	Ba = reshape(B*w,p-1,:)'
	Bs = reshape(param.lambda * vec(sum(B,2)), p-1,:)'

	blockEntropies = zeros(m)

	for k=1:m
		blockidx = (k-1)*(p-1)+1 : k*(p-1)
		Jk = - vec(sum(Rp[blockidx]'.*abs.(Ba.-d[blockidx]').^pwr,2)/pwr - sum(Bs,2))
		Jk-= maximum(Jk)
		Prob = exp.(Jk)
		Prob/=sum(Prob)

		nonzero = Prob .> 0

		blockEntropies[k] = abs(dot(Prob[nonzero], log2.(Prob[nonzero])))
        end

	entropy = sum(blockEntropies)

	return entropy
end


text:
function hexwordsum(string)
    sum( split(string,' ') ) do s
        s = replace(s, 'S' => '5' )
        s = replace(s, 'O' => '0' )
        if !isempty(s) && all(isxdigit,s)
            parse(Int,s,base=16)
        else
            0
        end
    end
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
###
# FB Oracle
###
# At the moment only supports cutting planes

export FBOracle

type FBOracle <: AbstractOracle
    mfs::Vector{Float64}
    mbs::Vector{Float64}
    sigfs::Vector{Float64}
    sigbs::Vector{Float64}
    log_eps::Float64

    # Cutting plane algorithm
    cut_tol::Float64  ##defaults to 1e-6

    # Other options
    debug_printcut::Bool
end

FBOracle(mfs, mbs, sigfs, sigbs, eps; TOL=1e-6) =
    FBOracle(mfs, mbs, sigfs, sigbs, log(1/eps), TOL, false)

#Preferred constructors
function FBOracle(data, eps, alpha1, alpha2; CUT_TOL=1e-6, numBoots=10_000)
    N, d  = size(data)
    mfs   = zeros(Float64, d)
    mbs   = zeros(Float64, d)
    sigfs = zeros(Float64, d)
    sigbs = zeros(Float64, d)

    for i = 1:d
        mbs[i], mfs[i] = calcMeansT(data[:, i], alpha1/d)
        sigfs[i], sigbs[i] = calcSigsBoot(data[:, i], alpha2/d, numBoots)
    end
    FBOracle(mfs, mbs, sigfs, sigbs, eps, TOL=CUT_TOL)
end
FBOracle(data, eps, alpha; CUT_TOL=1e-6, numBoots=10_000) =
    FBOracle(data, eps, alpha/2, alpha/2, CUT_TOL=CUT_TOL, numBoots=numBoots)

suppFcn(xs, w::FBOracle, cut_sense) =
    suppFcnFB(xs, w.mfs, w.mbs, w.sigfs, w.sigbs, w.log_eps, cut_sense)

#log_eps = log(1/eps_)
#returns zstar, ustar
function suppFcnFB(xs, mfs, mbs, sigfs, sigbs, log_eps, cut_sense=:Max)
    sign_flip = 1
    if cut_sense == :Min
        xs = -xs
        sign_flip = -1
    end
    lam = 0.0
    for i = 1:length(xs)
        if xs[i] >= 0
            lam += sigfs[i]^2 * xs[i]^2
        else
            lam += sigbs[i]^2 * xs[i]^2
        end
    end
    lam = sqrt(lam / (2 * log_eps))
    ustar = zeros(length(xs))
    for i = 1:length(xs)
        if xs[i] >= 0
            ustar[i] = mfs[i] + xs[i] * sigfs[i]^2/lam
        else
            ustar[i] = mbs[i] + xs[i] * sigbs[i]^2/lam
        end
    end
    zstar = dot(xs, ustar) * sign_flip
    return zstar, ustar
end

function setup(w::FBOracle, rm::Model, prefs)
    # Extract preferences we care about
    w.debug_printcut = get(prefs, :debug_printcut, false)
    w.cut_tol        = get(prefs, :cut_tol, w.cut_tol)

    rd = JuMPeR.getRobust(rm)
    @assert (rd.numUncs == length(w.mfs)) "Num Uncertainties $(rd.numUncs) doesn't match columns in data $(size(w.mfs))"
    @assert (length(w.mfs) == length(w.mbs) == length(w.sigfs) == length(w.sigbs)) "Lengths of means and fb devs dont match uncertainties"

    #ignore any additional constraints on uncertainties for now
    @assert (length(rd.uncertaintyset) == 0) "Auxiliary constraints on uncertainties not yet supported"
end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 04-2019
# Define operations for 2D spaces
#--------------------------------------------------------------------

import Base: maximum, minimum, range, size, ndims

ndims(PS::T)   where {T<:ProductSpace{S1, S2}} where {S1, S2} = 2
minimum(PS::T) where {T<:ProductSpace{S1, S2}} where {S1, S2} = (minimum(PS.S1), minimum(PS.S2))
maximum(PS::T) where {T<:ProductSpace{S1, S2}} where {S1, S2} = (maximum(PS.S1), maximum(PS.S2))
order(PS::T)   where {T<:ProductSpace{S1, S2}} where {S1, S2} = (order(PS.S1), order(PS.S2))
size(PS::T)    where {T<:ProductSpace{S1, S2}} where {S1, S2} = (size(PS.S1),   size(PS.S2))
range(PS::T)   where {T<:ProductSpace{S1, S2}} where {S1, S2} = (minimum(PS),   maximum(PS))

function Base. kron(A::Operator{S1}, B::Operator{S2})::Operator{ProductSpace{S2, S1}} where {S2, S1}
    return reshape(ProductSpace(B.space, A.space), kron(reshape(A), reshape(B)))
end

function Base. identity(PS::ProductSpace{S1, S2})::Operator{ProductSpace{S1, S2}} where {S1, S2}
    return kron(identity(PS.S2), identity(PS.S1))
end

function incomingboundary(PS::ProductSpace{S1, S2})::Operator{ProductSpace{S1,  S2}} where {S1, S2}
    return kron(identity(PS.S2), incomingboundary(PS.S1)) ⊕ kron(incomingboundary(PS.S2), identity(PS.S1))
end

function outgoingboundary(PS::ProductSpace{S1, S2})::Operator{ProductSpace{S1,  S2}} where {S1, S2}
    return kron(identity(PS.S2), outgoingboundary(PS.S1)) ⊕ kron(outgoingboundary(PS.S2), identity(PS.S1))
end

function derivative(PS::ProductSpace{S1, S2})::NTuple{2, Operator{ProductSpace{S1,  S2}}} where {S1, S2}
    return (kron(identity(PS.S2), derivative(PS.S1)), kron(derivative(PS.S2), identity(PS.S1)))
end

function integral(PS::ProductSpace{S1, S2})::Operator{ProductSpace{S1,  S2}}  where {S1, S2}
    return kron(integral(PS.S2), integral(PS.S1))
end

function Field(PS::ProductSpace{S1, S2}, umap::Function)::Field{ProductSpace{S1, S2}} where {S1, S2 <: Cardinal{Tag, N, T}} where {Tag, N, T}
    value = zeros(T, size(PS))
    for index in CartesianIndices(value)
        value[index] = umap(collocation(PS.S1, index.I[1]),
                            collocation(PS.S2, index.I[2]))
    end
    return Field(PS, value)
end


text:
aux_field(::BernoulliLikelihood{<:LogisticLink}, Ω) = getproperty(Ω, :ω)

function init_aux_variables(rng::AbstractRNG, ::BernoulliLikelihood{<:LogisticLink}, n::Int)
    return TupleVector(; ω=rand(rng, PolyaGamma(1, 0.0), n))
end

function init_aux_posterior(T::DataType, ::BernoulliLikelihood{<:LogisticLink}, n::Int)
    return For(TupleVector(; c=zeros(T, n))) do φ
        NTDist(PolyaGamma(1, φ.c))
    end
end

function aux_full_conditional(::BernoulliLikelihood{<:LogisticLink}, ::Any, f::Real)
    return NTDist(PolyaGamma(1, abs(f)))
end

function aux_posterior!(
    qΩ,
    ::BernoulliLikelihood{<:LogisticLink},
    ::AbstractVector,
    qf::AbstractVector{<:Normal},
)
    map!(sqrt ∘ second_moment, only(qΩ.inds).c, qf)
    return qΩ
end

function auglik_potential(::BernoulliLikelihood{<:LogisticLink}, ::Any, y::AbstractVector)
    return (sign.(y .- 0.5) / 2,)
end

function auglik_precision(::BernoulliLikelihood{<:LogisticLink}, Ω, ::AbstractVector)
    return (Ω.ω,)
end

function expected_auglik_potential(
    lik::BernoulliLikelihood{<:LogisticLink}, qΩ, y::AbstractVector
)
    return auglik_potential(lik, qΩ, y)
end

function expected_auglik_precision(
    ::BernoulliLikelihood{<:LogisticLink}, qΩ, ::AbstractVector
)
    return (tvmean(qΩ).ω,)
end

function logtilt(::BernoulliLikelihood{<:LogisticLink}, ω::Real, y::Real, f::Real)
    return -log(2) + (sign(y - 0.5) * f - abs2(f) * ω) / 2
end

function aux_prior(lik::BernoulliLikelihood{<:LogisticLink}, y::AbstractVector)
    return For(length(y)) do _
        NTDist(aux_prior(lik))
    end
end

aux_prior(::BernoulliLikelihood{<:LogisticLink}) = PolyaGamma(1, 0.0)

function expected_logtilt(
    ::BernoulliLikelihood{<:LogisticLink}, qωᵢ::NTDist{<:PolyaGamma}, yᵢ::Real, qfᵢ::Normal
)
    m = mean(qfᵢ)
    θ = ntmean(qωᵢ)
    return -log(2) + (sign(yᵢ - 0.5) * m - (abs2(m) + var(qfᵢ)) * θ.ω) / 2
end


text:
############################################################################################
# Types to distinguish between Gibbs and Metropolis updating
"Referencing struct: no updating via reference index"
struct Marginal <: ParticleReferencing end
"Referencing struct: last element in particle history always reference index"
struct Conditional <: ParticleReferencing end
"Referencing struct: last element in particle history weighted via reference index"
struct Ancestral <: ParticleReferencing end #PG with Ancestor Sampling

############################################################################################
"""
$(SIGNATURES)
Assign reference, depending on whether reference is tracked in Model.

# Examples
```julia
```

"""
function get_reference(referencing, model::ModelWrapper, data::D, tagged::Tagged) where {D}
    ## Check if only 1 parameter is tracked that has equal number of samples as data
    ArgCheck.@argcheck length(tagged.parameter) == 1 "Only 1 parameter can be tracked for Particle Filter."
    referenceᵖ = getindex(model.val, keys(tagged.parameter)[begin])
    return referenceᵖ
end
function get_reference(referencing, objective::Objective)
    return get_reference(referencing, objective.model, objective.data, objective.tagged)
end

############################################################################################
"""
$(SIGNATURES)
Inplace calculate ancestor weights.

# Examples
```julia
```

"""
function ancestralweight!(
    weights::BaytesCore.ParameterWeights,
    valₜ::Union{P,AbstractArray{P}},
    kernel::ParticleKernel,
    val::AbstractArray,
    iter::Integer,
) where {P}
    #!NOTE: buffer first overwritten in ℓtransition! and then assigned again on LHS
    weights.buffer .=
        weights.ℓweights .+ ℓtransition!(weights.buffer, valₜ, kernel, val, iter)
    weights.buffer .-= logsumexp(weights.buffer)
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Sample ancestor reference for particle history.

# Examples
```julia
```

"""
function sample_ancestor(
    _rng::Random.AbstractRNG,
    weights::BaytesCore.ParameterWeights,
    valₜ::Union{P,AbstractArray{P}},
    kernel::ParticleKernel,
    val::AbstractArray,
    iter::Integer,
) where {P}
    #!NOTE: buffer will be first overwritten on right side to inplace calculate new weights, and then assigned inplace on left side
    ancestralweight!(weights, valₜ, kernel, val, iter)
    #!NOTE Need to exp() particles.weights.buffer due to ancestralweights! calculation in log space
    weights.buffer .= exp.(weights.buffer)
    return BaytesCore.randcat(_rng, weights.buffer)
end

############################################################################################
"""
$(SIGNATURES)
Resampling function for particles, dispatched on ParticleReferencing subtype.

# Examples
```julia
```

"""
function ancestors!(
    _rng::Random.AbstractRNG,
    referencing::ParticleReferencing,
    resampling::R,
    ancestor::AbstractArray, #referencing::Marginal
    iter::Integer,
    Nparticles::Integer,
    weights::AbstractVector{F},
) where {R<:BaytesCore.ResamplingMethod,F<:AbstractFloat}
    #!NOTE: Weights.buffer already exp.(weightsₙ) during BaytesCore.issmaller( BaytesCore.computeESS(particles.weights), X) in previous step
    resample!(
        _rng,
        resampling,
        ancestor,
        iter,
        weights,
        isa(referencing, Marginal) ? Nparticles : Nparticles - 1,
    )
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Set last particle at current iteration.

# Examples
```julia
```

"""
function set_reference!(
    referencing::Marginal,
    val::AbstractVector{P},
    reference::AbstractVector{P},
    iter::Integer,
) where {P}
    #!NOTE: Nothing to be done in Marginal case
    return nothing
end
function set_reference!(
    referencing::R, val::AbstractVector{P}, reference::AbstractVector{P}, iter::Integer
) where {R<:Union{Conditional,Ancestral},P}
    #!NOTE: Reference particle at t is the same for PG and PGAS
    val[end] = reference[iter]
    return nothing
end
function set_reference!(
    referencing, val::AbstractMatrix{P}, reference::AbstractVector{P}, iter::Integer
) where {P}
    return set_reference!(referencing, view(val, :, iter), reference, iter)
end

############################################################################################
"""
$(SIGNATURES)
Compute new reference index from current reference index.

# Examples
```julia
```

"""
function get_reference!(
    _rng::Random.AbstractRNG,
    referencing::Marginal,
    ancestor::AbstractVector{I},
    particles::F,
    tune::ParticleFilterTune,
    reference::Union{P,AbstractArray{P}}
) where {I<:Integer, F<:AbstractParticles, P}
    return nothing
end

function get_reference!(
    _rng::Random.AbstractRNG,
    referencing::Conditional,
    ancestor::AbstractVector{I},
    particles::F,
    tune::ParticleFilterTune,
    reference::Union{P,AbstractArray{P}}
) where {I<:Integer, F<:AbstractParticles, P}
    ancestor[end] = length(ancestor)
    return nothing
end

function get_reference!(
    _rng::Random.AbstractRNG,
    referencing::Ancestral,
    ancestor::AbstractVector{I},
    particles::F,
    tune::ParticleFilterTune,
    reference::Union{P,AbstractArray{P}}
) where {I<:Integer, F<:AbstractParticles, P}
    referenceₜ = BaytesCore.grab(reference, tune.iter.current, tune.config.particle)
    #!NOTE: Weights.buffer already exp.(weightsₙ) during BaytesCore.issmaller( BaytesCore.computeESS(particles.weights), X) in previous step
    ancestor[end] = sample_ancestor(
        _rng,
        particles.weights,
        referenceₜ,
        particles.kernel,
        particles.val,
        tune.iter.current
    )
    return nothing
end

#!NOTE: If buffer is a matrix, reference will be assigned at iter-1 within function, but input here is iter as we grab data for forward looking ancestral step.
function get_reference!(
    _rng::Random.AbstractRNG,
    referencing::R,
    ancestor::AbstractMatrix{I},
    particles::F,
    tune::ParticleFilterTune,
    reference::Union{P,AbstractArray{P}}
) where {R<:ParticleReferencing, I<:Integer, F<:AbstractParticles, P}
    return get_reference!(
        _rng,
        referencing,
        view(ancestor, :, tune.iter.current - 1),
        particles,
        tune,
        reference
    )
end
function get_reference!(
    _rng::Random.AbstractRNG,
    particles::F,
    tune::ParticleFilterTune,
    reference::Union{P,AbstractArray{P}}
) where {F<:AbstractParticles, P}
    return get_reference!(_rng, tune.referencing, particles.ancestor, particles, tune, reference)
end

############################################################################################
#export
export
    Marginal,
    Conditional,
    Ancestral,
    get_reference,
    get_reference!,
    set_reference!,
    sample_ancestor


text:
include("abstract_model2fit_peak.jl")

include("gaussian.jl")

include("motif.jl")


text:
function train!(NN::NeuralNetwork, n_epochs::Int)
    for epoch = 1:n_epochs
        NN.epoch += 1
        push!(NN.J, 0)

        mini_batches = random_mini_batches(NN)
        for mini_batch in mini_batches
            propagate_forward!(NN, mini_batch)
            compute_cost!(NN, mini_batch)
            propagate_back!(NN, mini_batch)
            update_params!(NN)
        end
    end
end

function predict(params::Parameters, X::Array{T,2}) where {T<:AbstractFloat}
    X = eltype(params.W[1]).(X)
    gs = params.activation_functions
    L = length(params.layer_sizes)

    A = gs[1](params.W[1] * X .+ params.b[1])
    for l = 2:L
        A = gs[l](params.W[l] * A .+ params.b[l])
    end

    return A
end

function predict(NN::NeuralNetwork, X::Array{T,2}) where {T<:AbstractFloat}
    predict(NN.params, X)
end

function random_mini_batches(NN::NeuralNetwork)
    m = size(NN.X, 2)
    K = ceil(m / NN.hparams.mini_batch_size)
    permutation = randperm(NN.rng, m)
    mini_batches = MiniBatch[]

    for k = 1:K-1
        a = Int((k - 1) * NN.hparams.mini_batch_size + 1)
        b = Int(k * NN.hparams.mini_batch_size)
        push!(
            mini_batches,
            MiniBatch(NN.X[:, permutation][:, a:b], NN.Y[:, permutation][:, a:b]),
        )
    end

    a = Int((K - 1) * NN.hparams.mini_batch_size + 1)
    push!(
        mini_batches,
        MiniBatch(NN.X[:, permutation][:, a:end], NN.Y[:, permutation][:, a:end]),
    )

    return mini_batches
end

function one_cold(Y)
    Y_out = zeros(Float32, size(Y, 2))
    for i = 1:size(Y, 2)
        Y_out[i] = findmax(Y[:, i])[2]
    end
    return Y_out
end

function one_hot(Y, n_classes)
    Y_out = zeros(Float32, n_classes, size(Y, 2))
    for i = 1:size(Y, 2)
        Y_out[Y[i]+1, i] = 1
    end
    return Y_out
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashwindowcloselistener

"""
    dashwindowcloselistener(;kwargs...)

A DashWindowCloseListener component.
ExampleComponent is an example component.
It takes a property, `label`, and
displays it.
It renders an input with the property `value`
which is editable by the user.
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `status` (String; optional): The status of the window. Mounted or closing.
"""
function dashwindowcloselistener(; kwargs...)
        available_props = Symbol[:id, :status]
        wild_props = Symbol[]
        return Component("dashwindowcloselistener", "DashWindowCloseListener", "dwcl", available_props, wild_props; kwargs...)
end

max_stars_repo_name: raymon-ai/dwcl

max_issues_repo_name: raymon-ai/dwcl

max_forks_repo_name: raymon-ai/dwcl


text:
Pkg.add("Images")
Pkg.add("DataFrames")
using Images
using DataFrames

#typeData could be either "train" or "test.
#labelsInfo should contain the IDs of each image to be read
#The images in the trainResized and testResized data files
#are 20x20 pixels, so imageSize is set to 400.
#path should be set to the location of the data files.

function read_data(typeData, labelsInfo, imageSize, path)
 #Intialize x matrix
 x = zeros(size(labelsInfo, 1), imageSize)

 for (index, idImage) in enumerate(labelsInfo["ID"])
  #Read image file
  nameFile = "$(path)/$(typeData)Resized/$(idImage).Bmp"
  img = imread(nameFile)

  #Convert img to float values
  temp = float32sc(img)

  #Convert color images to gray images
  #by taking the average of the color scales.
  if ndims(temp) == 3
   temp = mean(temp.data, 1)
  end

  #Transform image matrix to a vector and store
  #it in data matrix
  x[index, :] = reshape(temp, 1, imageSize)
 end
 return x
end

imageSize = 400 # 20 x 20 pixel

#Set location of data files, folders
path = ...

#Read information about training data , IDs.
labelsInfoTrain = readtable("$(path)/trainLabels.csv")

#Read training matrix
xTrain = read_data("train", labelsInfoTrain, imageSize, path)

#Read information about test data ( IDs ).
labelsInfoTest = readtable("$(path)/sampleSubmission.csv")

#Read test matrix
xTest = read_data("test", labelsInfoTest, imageSize, path)

#Get only first character of string (convert from string to character).
#Apply the function to each element of the column "Class"
yTrain = map(x -> x[1], labelsInfoTrain["Class"])

#Convert from character to integer
yTrain = int(yTrain)

xTrain = xTrain'
xTest = xTest'

addprocs(3)

@everywhere function euclidean_distance(a, b)
 distance = 0.0
 for index in 1:size(a, 1)
  distance += (a[index]-b[index]) * (a[index]-b[index])
 end
 return distance
end

@everywhere function get_k_nearest_neighbors(x, i, k)
 nRows, nCols = size(x)
 imageI = Array(Float32, nRows)
 for index in 1:nRows
  imageI[index] = x[index, i]
 end
 imageJ = Array(Float32, nRows)
 distances = Array(Float32, nCols)
 for j in 1:nCols
  for index in 1:nRows
   imageJ[index] = x[index, j]
  end
  distances[j] = euclidean_distance(imageI, imageJ)
 end
 sortedNeighbors = sortperm(distances)
 kNearestNeighbors = sortedNeighbors[2:k+1]
 return kNearestNeighbors
end

@everywhere function assign_label(x, y, k, i)
 kNearestNeighbors = get_k_nearest_neighbors(x, i, k)
 counts = Dict{Int, Int}()
 highestCount = 0
 mostPopularLabel = 0
 for n in kNearestNeighbors
  labelOfN = y[n]
  if !haskey(counts, labelOfN)
   counts[labelOfN] = 0
  end
  counts[labelOfN] += 1
  if counts[labelOfN] > highestCount
   highestCount = counts[labelOfN]
   mostPopularLabel = labelOfN
  end
 end
 return mostPopularLabel
end

tic()
k = 1
sumValues = @parallel (+) for i in 1:size(xTrain, 2)
 assign_label(xTrain, yTrain, k, i) == yTrain[i, 1]
end
loofCvAccuracy = sumValues / size(xTrain, 2)
println("The LOOF-CV accuracy of 1NN is $(loofCvAccuracy)")
toc()

#Similar to function assign_label.
#Only changes are commented
@everywhere function assign_label_each_k(x, y, maxK, i)
 kNearestNeighbors = get_k_nearest_neighbors(x, i, maxK)

 #The next array will keep the labels for each value of k
 labelsK = zeros(Int, 1, maxK)

 counts = Dict{Int, Int}()
 highestCount = 0
 mostPopularLabel = 0

 #We need to keep track of the current value of k
 for (k, n) in enumerate(kNearestNeighbors)
  labelOfN = y[n]
  if !haskey(counts, labelOfN)
   counts[labelOfN] = 0
  end
  counts[labelOfN] += 1
  if counts[labelOfN] > highestCount
   highestCount = counts[labelOfN]
   mostPopularLabel = labelOfN
  end
  #Save current most popular label
  labelsK[k] = mostPopularLabel
 end
 #Return vector of labels for each k
 return labelsK
end

tic()
maxK = 20 #Any value can be chosen
yPredictionsK = @parallel (vcat) for i in 1:size(xTrain, 2)
 assign_label_each_k(xTrain, yTrain, maxK, i)
end
for k in 1:maxK
 accuracyK = mean(yTrain .== yPredictionsK[:, k])
 println("The LOOF-CV accuracy of $(k)-NN is $(accuracyK)")
end
toc()

@everywhere function get_k_nearest_neighbors(xTrain, imageI, k)
 nRows, nCols = size(xTrain)
 imageJ = Array(Float32, nRows)
 distances = Array(Float32, nCols)
 for j in 1:nCols
  for index in 1:nRows
   imageJ[index] = xTrain[index, j]
  end
  distances[j] = euclidean_distance(imageI, imageJ)
 end
 sortedNeighbors = sortperm(distances)
 kNearestNeighbors = sortedNeighbors[1:k]
 return kNearestNeighbors
end

@everywhere function assign_label(xTrain, yTrain, k, imageI)
 kNearestNeighbors = get_k_nearest_neighbors(xTrain, imageI, k)
 counts = Dict{Int, Int}()
 highestCount = 0
 mostPopularLabel = 0
 for n in kNearestNeighbors
  labelOfN = yTrain[n]
  if !haskey(counts, labelOfN)
   counts[labelOfN] = 0
  end
  counts[labelOfN] += 1 #add one to the count
  if counts[labelOfN] > highestCount
   highestCount = counts[labelOfN]
   mostPopularLabel = labelOfN
  end
 end
 return mostPopularLabel
end

println("Running kNN on test data")
tic()
k = 3 # The CV accuracy shows this value to be the best.
yPredictions = @parallel (vcat) for i in 1:size(xTest, 2)
 nRows = size(xTrain, 1)
 imageI = Array(Float32, nRows)
 for index in 1:nRows
  imageI[index] = xTest[index, i]
 end
 assign_label(xTrain, yTrain, k, imageI)
end
toc()

#Convert integer predictions to character
labelsInfoTest["Class"] = char(yPredictions)

#Save predictions
writetable("$(path)/juliaKNNSubmission.csv", labelsInfoTest, separator=',', header=true)

println("Submission file saved in $(path)/juliaKNNSubmission.csv")

max_stars_repo_name: JimmyFW/stop-making-sense-julia

max_issues_repo_name: JimmyFW/stop-making-sense-julia

max_forks_repo_name: JimmyFW/stop-making-sense-julia


text:
using Lumira
using Gadfly

plot = plot([sin,cos], 5, 25)

exportToSVG(plot)

setId("trigo")
setName("HelloWorld1")
createLumiraExtension()


text:
# Testing ObjectId

# struct TestObj4 <: ModelObject
#     objId :: ObjectId
# end

# function single_id_test()
#     @testset "SingleId" begin
#         id1 = SingleId(:id1, [1, 2])
#         show(id1);
#         @test id1.index == [1,2]
#         @test has_index(id1)

#         id11 = SingleId(:id1, [1, 2]);
#         @test isequal(id1, id11)

#         id2 = SingleId(:id2, 3)
#         @test id2.index == [3]

#         id3 = SingleId(:id3);
#         @test isempty(id3.index)
#         @test !has_index(id3)
#         @test isequal(id3, SingleId(:id3))

#         @test isequal([id1, id2], [id1, id2])
#         @test !isequal([id1, id1], [id2, id1])
#         @test !isequal([id1, id2], [id1])

#         id4 = SingleId(:id4);
#         s4 = ModelParams.make_string(id4);
#         @test isequal(s4, "id4")
#         id4a = ModelParams.make_single_id(s4);
#         @test isequal(id4, id4a)

#         id5 = SingleId(:id5, [4, 2]);
#         s5 = ModelParams.make_string(id5);
#         id5a = ModelParams.make_single_id(s5);
#         @test isequal(id5, id5a)

#         id6 = SingleId(:id6, 4)
#         s6 = ModelParams.make_string(id6);
#         id6a = ModelParams.make_single_id(s6);
#         @test isequal(id6, id6a)
#     end
# end

# function object_id_test()
#     @testset "ObjectId" begin
#         # Simplest case
#         id1 = SingleId(:id1);
#         o1 = ObjectId(id1);
#         @test !ModelParams.has_parent(o1)

#         # Roundtrip between objectId and String
#         s1 = make_string(id1);
#         @test isequal(s1, "id1")
#         o1a = make_object_id(s1)
#         @test isequal(o1, o1a)

#         # Index, no parents
#         o2 = ObjectId(:id2, 2)
#         show(o2);
#         @test ModelParams.own_index(o2) == [2]
#         @test ModelParams.own_name(o2) == :id2

#         # Has id1 as parent
#         o3 = ObjectId(:id3, 2, o1);
#         show(o3)
#         p3 = ModelParams.get_parent_id(o3);
#         @test ModelParams.is_parent_of(p3, o3)
#         s3 = ModelParams.make_string(o3);
#         @test isequal(s3, "id1 > id3[2]")
#         o3a = make_object_id(s3)
#         @test isequal(o3, o3a)

#         # Make child id
#         o4 = ObjectId(:id4, p3);
#         obj4 = TestObj4(o4);
#         childId = ModelParams.make_child_id(obj4, :child);
#         @test isequal(ModelParams.own_id(childId),  SingleId(:child))
#         @test isequal(o4, ModelParams.get_parent_id(childId))
#         @test isequal(ModelParams.own_name(obj4), :id4)

#         # Check `isequal` when "depth" of `ObjectId`s is different
#         @test !isequal(obj4.objId, childId)
#     end
# end

# @testset "ObjectId" begin
#     single_id_test()
#     object_id_test()
# end

# -----------


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__

d_str = "DAG(u ~ x + w, x ~ v + z, y ~ v, s1 ~ u, w ~ y, k ~ y, s2 ~ w + s3 + k, s3 ~ s1)"

dag = DAG("test_desc", d_str)

to_dagitty(dag)

#=
fname = joinpath(ProjDir, "test_descendents_01.dot")
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)
=#

bs = basis_set(dag)
as = adjustment_sets(dag, :u, :s2)

@testset "Test_descendents_01" begin

  @test dag.a[:x, :u] == 0
  @test dag.e[:x, :u] == 1
  @test length(bs) == 83
  @test bs[70] == [:y, :s1, :u]
  @test bs[71] == [:y, :s1, :v, :u]
  @test bs[81] == [:z, :k, :y]
  @test bs[82] == [:z, :s2, :k, :w, :s3]
  @test bs[83] == [:z, :s3, :s1]

end

@testset "Adjustment_sets" begin

  @test as[1] == [:x, :w]
  @test as[2] == [:v, :w]
  @test as[3] == [:y, :w]
  @test as[4] == [:w, :k]

end