text:
#    Bessel functions of the first kind of order zero and one
#                       besselj0, besselj1
#
#    Calculation of besselj0 is done in three branches using polynomial approximations
#
#    Branch 1: x <= 5.0
#              besselj0 = (x^2 - r1^2)*(x^2 - r2^2)*P3(x^2) / Q8(x^2)
#    where r1 and r2 are zeros of J0
#    and P3 and Q8 are a 3 and 8 degree polynomial respectively
#    Polynomial coefficients are from [1] which is based on [2]
#    For tiny arugments the power series expansion is used.
#
#    Branch 2: 5.0 < x < 25.0
#              besselj0 = sqrt(2/(pi*x))*(cos(x - pi/4)*R7(x) - sin(x - pi/4)*R8(x))
#    Hankel's asymptotic expansion is used
#    where R7 and R8 are rational functions (Pn(x)/Qn(x)) of degree 7 and 8 respectively
#    See section 4 of [3] for more details and [1] for coefficients of polynomials
#
#   Branch 3: x >= 25.0
#              besselj0 = sqrt(2/(pi*x))*beta(x)*(cos(x - pi/4 - alpha(x))
#   See modified expansions given in [3]. Exact coefficients are used
#
#   Calculation of besselj1 is done in a similar way as besselj0.
#   See [3] for details on similarities.
#
# [1] https://github.com/deepmind/torch-cephes
# [2] Cephes Math Library Release 2.8:  June, 2000 by Stephen L. Moshier
# [3] Harrison, John. "Fast and accurate Bessel function computation."
#     2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
#

"""
    besselj0(x::T) where T <: Union{Float32, Float64}

Bessel function of the first kind of order zero, ``J_0(x)``.
"""
function besselj0(x::Float64)
    T = Float64
    x = abs(x)
    isinf(x) && return zero(x)

    if x <= 5
        z = x * x
        if x < 1.0e-5
            return 1.0 - z / 4.0
        end
        DR1 = 5.78318596294678452118e0
        DR2 = 3.04712623436620863991e1
        p = (z - DR1) * (z - DR2)
        p = p * evalpoly(z, RP_j0(T)) / evalpoly(z, RQ_j0(T))
        return p
    elseif x < 25.0
        w = 5.0 / x
        q = 25.0 / (x * x)

        p = evalpoly(q, PP_j0(T)) / evalpoly(q, PQ_j0(T))
        q = evalpoly(q, QP_j0(T)) / evalpoly(q, QQ_j0(T))
        xn = x - PIO4(T)
        sc = sincos(xn)
        p = p * sc[2] - w * q * sc[1]
        return p * SQ2OPI(T) / sqrt(x)
    else
        if x < 120.0
            p = (one(T), -1/16, 53/512, -4447/8192, 3066403/524288, -896631415/8388608, 796754802993/268435456, -500528959023471/4294967296)
            q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296, 24713030909/46137344, -7780757249041/436207616)
        else
            p = (one(T), -1/16, 53/512, -4447/8192)
            q = (-1/8, 25/384, -1073/5120, 375733/229376)
        end
        xinv = inv(x)
        x2 = xinv*xinv

        p = evalpoly(x2, p)
        a = SQ2OPI(T) * sqrt(xinv) * p
        xn = muladd(xinv, evalpoly(x2, q), - PIO4(T))

        # the following computes b = cos(x + xn) more accurately
        # see src/misc.jl
        b = cos_sum(x, xn)
        return a * b
    end
end
function besselj0(x::Float32)
    T = Float32
    x = abs(x)
    isinf(x) && return zero(x)

    if x <= 2.0f0
        z = x * x
        if x < 1.0f-3
            return 1.0f0 - 0.25f0 * z
        end
        DR1 = 5.78318596294678452118f0
        p = (z - DR1) * evalpoly(z, JP_j0(T))
        return p
    else
        q = inv(x)
        w = sqrt(q)
        p = w * evalpoly(q, MO_j0(T))
        w = q * q
        xn = q * evalpoly(w, PH_j0(T)) - PIO4(Float32)
        p = p * cos(xn + x)
        return p
    end
end

"""
    besselj1(x::T) where T <: Union{Float32, Float64}

Bessel function of the first kind of order one, ``J_1(x)``.
"""
function besselj1(x::Float64)
    T = Float64
    x = abs(x)
    isinf(x) && return zero(x)

    if x <= 5.0
        z = x * x
        w = evalpoly(z, RP_j1(T)) / evalpoly(z, RQ_j1(T))
        w = w * x * (z - 1.46819706421238932572e1) * (z - 4.92184563216946036703e1)
        return w
    elseif x < 25.0
        w = 5.0 / x
        z = w * w
        p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
        q = evalpoly(z, QP_j1(T)) / evalpoly(z, QQ_j1(T))
        xn = x - THPIO4(T)
        sc = sincos(xn)
        p = p * sc[2] - w * q * sc[1]
        return p * SQ2OPI(T) / sqrt(x)
    else
        if x < 120.0
            p = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296)
            q = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144, -30413055339/46137344, 9228545313147/436207616)
        else
            p = (one(T), 3/16, -99/512, 6597/8192)
            q = (3/8, -21/128, 1899/5120, -543483/229376)
        end
        xinv = inv(x)
        x2 = xinv*xinv

        p = evalpoly(x2, p)
        a = SQ2OPI(T) * sqrt(xinv) * p
        xn = muladd(xinv, evalpoly(x2, q), - 3 * PIO4(T))

        # the following computes b = cos(x + xn) more accurately
        # see src/misc.jl
        b = cos_sum(x, xn)
        return a * b
    end
end
function besselj1(x::Float32)
    x = abs(x)
    isinf(x) && return zero(x)

    if x <= 2.0f0
        z = x * x
        Z1 = 1.46819706421238932572f1
        p = (z - Z1) * x * evalpoly(z, JP32)
        return p
    else
        q = inv(x)
        w = sqrt(q)
        p = w * evalpoly(q, MO132)
        w = q * q
        xn = q * evalpoly(w, PH132) - THPIO4(Float32)
        p = p * cos(xn + x)
        return p
    end
end


text:
#Benchmarking of RBC and FVGQ variants
using DifferenceEquations, BenchmarkTools
using DelimitedFiles, Distributions, Zygote, LinearAlgebra

# for benchmarking construction itself
function make_problem_1(A, B, C, u0, noise, observables, D; kwargs...)
    prob = LinearStateSpaceProblem(A, B, u0, (0, size(observables, 2)); C, observables_noise = D,
                                   noise, observables, kwargs...)
    return prob.A[1, 1] + prob.B[1, 1]
end
function make_problem_kalman(A, B, C, u0_prior, observables, D; kwargs...)
    prob = LinearStateSpaceProblem(A, B, u0_prior, (0, size(observables, 2)); C, u0_prior,
                                   observables_noise = D, noise = nothing, observables, kwargs...)
    return prob.A[1, 1] + prob.B[1, 1]
end

# for benchmarking likelihoods
function joint_likelihood_1(A, B, C, u0, noise, observables, D; kwargs...)
    prob = LinearStateSpaceProblem(A, B, u0, (0, size(observables, 2)); C, observables_noise = D,
                                   noise, observables, kwargs...)
    return solve(prob).logpdf
end
function kalman_likelihood(A, B, C, u0_prior, observables, D; kwargs...)
    prob = LinearStateSpaceProblem(A, B, u0_prior, (0, size(observables, 2)); C, u0_prior,
                                   observables_noise = D, noise = nothing, observables, kwargs...)
    return solve(prob).logpdf
end

function simulate_model_no_noise_1(A, B, C, u0, observables, D; kwargs...)
    prob = LinearStateSpaceProblem(A, B, u0, (0, size(observables, 2)); C, observables_noise = D,
                                   observables, kwargs...)
    sol = solve(prob)
    return sol.retcode
end

function simulate_model_no_observations_1(A, B, C, u0, T; kwargs...)
    prob = LinearStateSpaceProblem(A, B, u0, (0, T); C, kwargs...)
    sol = solve(prob)
    return sol.retcode
end

# Matrices from RBC
const A_rbc = [0.9568351489231076 6.209371005755285;
               3.0153731819288737e-18 0.20000000000000007]
const B_rbc = reshape([0.0; -0.01], 2, 1) # make sure B is a matrix
const C_rbc = [0.09579643002426148 0.6746869652592109; 1.0 0.0]
const D_rbc = MvNormal(Diagonal(abs2.([0.1, 0.1])))
const u0_rbc = zeros(2)
const u0_prior_rbc = MvNormal(diagm(ones(length(u0_rbc))))

const observables_rbc = readdlm(joinpath(pkgdir(DifferenceEquations),
                                         "test/data/RBC_observables.csv"), ',')' |> collect
const noise_rbc = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/RBC_noise.csv"), ',')' |>
                  collect
const T_rbc = size(observables_rbc, 2)
# Matrices from FVGQ
# Load FVGQ data for checks
const A_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_A.csv"), ',')
const B_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_B.csv"), ',')
const C_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_C.csv"), ',')
const D_FVGQ = MvNormal(Diagonal(abs2.(ones(6) * 1e-3)))

const observables_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations),
                                          "test/data/FVGQ20_observables.csv"), ',')' |> collect

const noise_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_noise.csv"),
                           ',')' |> collect
const u0_FVGQ = zeros(size(A_FVGQ, 1))
const u0_prior_FVGQ = MvNormal(diagm(ones(length(u0_FVGQ))))
const T_FVGQ = size(observables_FVGQ, 2)
# executing gradients once to avoid compilation time in benchmarking

make_problem_1(A_rbc, B_rbc, C_rbc, u0_rbc, noise_rbc, observables_rbc, D_rbc)
gradient((args...) -> make_problem_1(args..., observables_rbc, D_rbc), A_rbc, B_rbc, C_rbc, u0_rbc,
         noise_rbc)

make_problem_kalman(A_rbc, B_rbc, C_rbc, u0_prior_rbc, observables_rbc, D_rbc)
gradient((args...) -> make_problem_kalman(args..., observables_rbc, D_rbc), A_rbc, B_rbc, C_rbc,
         u0_prior_rbc)

kalman_likelihood(A_rbc, B_rbc, C_rbc, u0_prior_rbc, observables_rbc, D_rbc)
gradient((args...) -> kalman_likelihood(args..., observables_rbc, D_rbc), A_rbc, B_rbc, C_rbc,
         u0_prior_rbc)
joint_likelihood_1(A_rbc, B_rbc, C_rbc, u0_rbc, noise_rbc, observables_rbc, D_rbc)
gradient((args...) -> joint_likelihood_1(args..., observables_rbc, D_rbc), A_rbc, B_rbc, C_rbc,
         u0_rbc, noise_rbc)

kalman_likelihood(A_FVGQ, B_FVGQ, C_FVGQ, u0_prior_FVGQ, observables_FVGQ, D_FVGQ)
gradient((args...) -> kalman_likelihood(args..., observables_FVGQ, D_FVGQ), A_FVGQ, B_FVGQ, C_FVGQ,
         u0_prior_FVGQ)
joint_likelihood_1(A_FVGQ, B_FVGQ, C_FVGQ, u0_FVGQ, noise_FVGQ, observables_FVGQ, D_FVGQ)
gradient((args...) -> joint_likelihood_1(args..., observables_FVGQ, D_FVGQ), A_FVGQ, B_FVGQ, C_FVGQ,
         u0_FVGQ, noise_FVGQ)

####### Becnmarks

const LINEAR = BenchmarkGroup()

const LINEAR["rbc"] = BenchmarkGroup()

const LINEAR["rbc"]["make_problem_1"] = @benchmarkable make_problem_1($A_rbc, $B_rbc, $C_rbc,
                                                                      $u0_rbc, $noise_rbc,
                                                                      $observables_rbc, $D_rbc)
const LINEAR["rbc"]["make_problem_1_gradient"] = @benchmarkable gradient((args...) -> make_problem_1(args...,
                                                                                $observables_rbc,
                                                                                $D_rbc),
                                                                         $A_rbc, $B_rbc, $C_rbc,
                                                                         $u0_rbc, $noise_rbc)

const LINEAR["rbc"]["simulate_model_no_noise_1"] = @benchmarkable simulate_model_no_noise_1($A_rbc,
                                                                                $B_rbc,
                                                                                $C_rbc,
                                                                                $u0_rbc,
                                                                                $observables_rbc,
                                                                                $D_rbc)
const LINEAR["rbc"]["simulate_model_no_observations_1"] = @benchmarkable simulate_model_no_observations_1($A_rbc,
                                                                                $B_rbc,
                                                                                $C_rbc,
                                                                                $u0_rbc,
                                                                                T_rbc)
const LINEAR["rbc"]["joint_1"] = @benchmarkable joint_likelihood_1($A_rbc, $B_rbc, $C_rbc, $u0_rbc,
                                                                   $noise_rbc, $observables_rbc,
                                                                   $D_rbc)
const LINEAR["rbc"]["joint_1_gradient"] = @benchmarkable gradient((args...) -> joint_likelihood_1(args...,
                                                                                $observables_rbc,
                                                                                $D_rbc),
                                                                  $A_rbc, $B_rbc, $C_rbc, $u0_rbc,
                                                                  $noise_rbc)
const LINEAR["rbc"]["make_problem_kalman"] = @benchmarkable make_problem_kalman($A_rbc, $B_rbc,
                                                                                $C_rbc,
                                                                                $u0_prior_rbc,
                                                                                $observables_rbc,
                                                                                $D_rbc)
const LINEAR["rbc"]["make_problem_kalman_gradient"] = @benchmarkable gradient((args...) -> make_problem_kalman(args...,
                                                                                $observables_rbc,
                                                                                $D_rbc),
                                                                              $A_rbc, $B_rbc,
                                                                              $C_rbc, $u0_prior_rbc)
const LINEAR["rbc"]["kalman"] = @benchmarkable kalman_likelihood($A_rbc, $B_rbc, $C_rbc,
                                                                 $u0_prior_rbc, $observables_rbc,
                                                                 $D_rbc)
const LINEAR["rbc"]["kalman_gradient"] = @benchmarkable gradient((args...) -> kalman_likelihood(args...,
                                                                                $observables_rbc,
                                                                                $D_rbc),
                                                                 $A_rbc, $B_rbc, $C_rbc,
                                                                 $u0_prior_rbc)

# FVGQ
const LINEAR["FVGQ"] = BenchmarkGroup()
const LINEAR["FVGQ"]["make_problem_1"] = @benchmarkable make_problem_1($A_FVGQ, $B_FVGQ, $C_FVGQ,
                                                                       $u0_FVGQ, $noise_FVGQ,
                                                                       $observables_FVGQ, $D_FVGQ)
const LINEAR["FVGQ"]["make_problem_1_gradient"] = @benchmarkable gradient((args...) -> make_problem_1(args...,
                                                                                $observables_FVGQ,
                                                                                $D_FVGQ),
                                                                          $A_FVGQ, $B_FVGQ, $C_FVGQ,
                                                                          $u0_FVGQ, $noise_FVGQ)
const LINEAR["FVGQ"]["simulate_model_no_noise_1"] = @benchmarkable simulate_model_no_noise_1($A_FVGQ,
                                                                                $B_FVGQ,
                                                                                $C_FVGQ,
                                                                                $u0_FVGQ,
                                                                                $observables_FVGQ,
                                                                                $D_FVGQ)
const LINEAR["FVGQ"]["simulate_model_no_observations_1"] = @benchmarkable simulate_model_no_observations_1($A_FVGQ,
                                                                                $B_FVGQ,
                                                                                $C_FVGQ,
                                                                                $u0_FVGQ,
                                                                                $T_FVGQ)
const LINEAR["FVGQ"]["joint_1"] = @benchmarkable joint_likelihood_1($A_FVGQ, $B_FVGQ, $C_FVGQ,
                                                                    $u0_FVGQ, $noise_FVGQ,
                                                                    $observables_FVGQ, $D_FVGQ)
const LINEAR["FVGQ"]["joint_1_gradient"] = @benchmarkable gradient((args...) -> joint_likelihood_1(args...,
                                                                                $observables_FVGQ,
                                                                                $D_FVGQ),
                                                                   $A_FVGQ, $B_FVGQ, $C_FVGQ,
                                                                   $u0_FVGQ, $noise_FVGQ)
const LINEAR["FVGQ"]["make_problem_kalman"] = @benchmarkable make_problem_kalman($A_FVGQ, $B_FVGQ,
                                                                                $C_FVGQ,
                                                                                $u0_prior_FVGQ,
                                                                                $observables_FVGQ,
                                                                                $D_FVGQ)
const LINEAR["FVGQ"]["make_problem_kalman_gradient"] = @benchmarkable gradient((args...) -> make_problem_kalman(args...,
                                                                                $observables_FVGQ,
                                                                                $D_FVGQ),
                                                                               $A_FVGQ, $B_FVGQ,
                                                                               $C_FVGQ,
                                                                               $u0_prior_FVGQ)
const LINEAR["FVGQ"]["kalman"] = @benchmarkable kalman_likelihood($A_FVGQ, $B_FVGQ, $C_FVGQ,
                                                                  $u0_prior_FVGQ, $observables_FVGQ,
                                                                  $D_FVGQ)
const LINEAR["FVGQ"]["kalman_gradient"] = @benchmarkable gradient((args...) -> kalman_likelihood(args...,
                                                                                $observables_FVGQ,
                                                                                $D_FVGQ),
                                                                  $A_FVGQ, $B_FVGQ, $C_FVGQ,
                                                                  $u0_prior_FVGQ)

# return for the test suite
LINEAR

ext: jl


text:
## eer.jl  Equal Error Rate and friends
## (c) 2014--2015 David A. van Leeuwen
##
## This code is licenced under the MIT software license

## eerch: computes the eer using the convex hull method
"""
`eerch()` computes the Equal Error Rate using the Convex Hull interpretation.  This computes the error
rate at which the convex hull of the ROC curve crosses the `y=x` diagonal in the ROC.

This value is less sensitive to trial sampling at low number of trials, and has the interpretation of
`the maximum while vary prior of the minimum decision costs'---a useful quantity in decision cost
function analysis of a two class classification system.

Arguments:

 - `r::Roc`: an object of type `Roc`, the result of `roc()`

 - `tar, non`: Vectors of target and non-target scores

 - `pfa, pmiss, ch`: Vectors of the false positive and false negative rate, and an array indicating
  the membership of the `(pfa, pmiss)` point on the convex hull.  These points are found by `roc()`.
"""
eerch(r::Roc) = eerch(r.pfa, r.pmiss, r.ch)
eerch(tar::Vector, non::Vector) = eerch(roc(tar, non))

## compute the EER given corresponding pfa and pmiss array, but use the
## convex hull points only.
function eerch(pfa::Vector{T}, pmiss::Vector{T}, ch::BitVector) where T<:AbstractFloat
    @assert length(pfa) == length(pmiss) == length(ch)
    ## find the index on the ch where the crossing of pfa=pmiss occurs
    chi = findall(ch)
    i = chi[1]
    direction = sign(pfa[i] - pmiss[i]) # >0 if first pfa > pmiss
    ## does it pay off to do a binary search?  It should, really.
    li = i
    for loop_variable in 1:length(chi)
        i = chi[loop_variable]
        if sign(pfa[i] - pmiss[i]) != direction
            break
        end
        li = i
    end
    ## compute the crossing
    (ax,bx) = pfa[[li,i]]
    (ay,by) = pmiss[[li,i]]
    return crossing(ax, ay, bx, by)
end

## compute crossing with y=x points (ax,ay) and (bx,by)
function crossing(ax, ay, bx, by)
    den = ax-ay-bx+by
    if den == 0
        return (ax + bx + ay + by) / 4
    else
        return ax + (ax-ay)*(bx-ax) / den
    end
end

## just a simple EER approximation, optimized for memory and speed
"""
`eer(tar, non)` computes a simple approximation to the Equal Error
Rate, the intersection of the ROC with the line `y=x`, i.e., the error
rate at which the false positive rate and the false negative rate are
approximately equal.

For a more consistent interpretation and implementation of the EER, please consider `eerch()`,
the Equal Error Rate - Convex Hull computation.
"""
eer(tar::Vector{T}, non::Vector{T}) where {T<:Real} = eer_sorted(sort(tar), sort(non))
eer(r::Roc) = eerch(r.pfa, r.pmiss, trues(length(r.pfa)))

## if tar and non are sorted, we may be doing even faster...
function eer_sorted(tar::Vector{T}, non::Vector{T}) where T<:Real
    ntar = length(tar)
    nnon = length(non)
    ntar > 0 && nnon > 0 || error("target and non-target arrays must not be empty")
    Δfa = 1. / nnon
    Δmiss = 1. / ntar
    pmiss = 0.
    ni = binsearch(tar[1], non, lower=true) # index in non-targets of lowest target score
    pfa = 1. - ni * Δfa
    ti = 1
    ni += 1
    lpfa, lpmiss = pfa, pmiss
    while ni ≤ nnon && ti ≤ ntar
        if tar[ti] < non[ni]
            pmiss += Δmiss
            ti += 1
        elseif tar[ti] > non[ni]
            pfa -= Δfa
            ni += 1
        else
            pmiss += Δmiss
            ti += 1
            pfa -= Δfa
            ni += 1
        end
        if pfa < pmiss
            break
        end
        lpfa, lpmiss = pfa, pmiss
    end
    return crossing(pfa, pmiss, lpfa, lpmiss)
end

## Returns the index to the largest value in the sorted array `a` ≤ `x` if lower==false
## If lower==true, the value must be strictly < `x`
function binsearch(x::Real, a::Vector{T}; lower=false, check=false) where T<:Real
    check && (issorted(a) || error("Array needs to be sorted"))
    mi = 1
    ma = length(a)
    if x < a[mi] || lower && x == a[mi]
        return 0
    elseif x > a[ma] || !lower && x == a[ma]
        return ma
    end
    while ma - mi > 1
        h = mi + div(ma-mi,2)
        if x > a[h] || !lower && x == a[h]
            mi = h
        else
            ma = h
        end
    end
    return mi
end

## EER, given the sort order of `[tar, non]`
## The idea is that you don't need to sort multiple times if you have
## multiple selections.  This cannot be accurate for cases where
## `tar` and `non` have the same value
function eer_so(so::Vector{T}, ntar::T, selection::BitVector) where T<:Integer
    length(so) == length(selection) || error("selection not the same length as sort order")
    eer(so[findin(so, find(selection))], ntar)
end

## This computes the actual approximate eer from sortorder and number of targets and nontargets
## The sort order "so" may be a subset of an original sort order, "thres" should be the
## original number of targets.
function eer_so(so::Vector{T}, thres::T, ntar::T=sum(so .<= thres)) where T<:Integer
    nnon = length(so) - ntar
    Δfa = 1. / nnon
    Δmiss = 1. / ntar
    pfa = 1.
    pmiss = 0.
    lpfa, lpmiss = pfa, pmiss
    for s in so
        if s <= thres            # target score
            pmiss += Δmiss
        else
            pfa -= Δfa
        end
        if pfa < pmiss
            break
        end
        lpfa, lpmiss = pfa, pmiss
    end
    return crossing(pfa, pmiss, lpfa, lpmiss)
end

## experimental routine for timing and optimization...
function eer_experimental(tar::Vector{T}, non::Vector{T}; method=:naive, fast=true) where T
    ntar = length(tar)
    nnon = length(non)
    @time scores = vcat(non,tar)      # 1
    @time so = sortperm(scores)       # 1
    @time truth = so .> nnon          # 1/64
    @time if !fast
        pfa = 1.0 - cumsum(!truth)/nnon
        pmiss = cumsum(truth)/ntar
    else                        # this takes less memory, but that's all I can say for it.
        Δ = 1/nnon
        pfa = similar(scores)   # 1
        s = one(T)
        for i=1:length(scores)
            s -= Δ*!truth[i]
            pfa[i] = s
        end
        Δ = 1/ntar
        pmiss = similar(scores) # 1
        s = zero(T)
        for i=1:length(scores)
            s += Δ*truth[i]
            pmiss[i] = s
        end
    end
    @time if method==:naive
        i = indmin(abs(pfa-pmiss))
        res = (pfa[i] + pmiss[i]) / 2
    end
    return res, (pfa, pmiss, scores[so]) # 1
end


text:
# This file was generated, do not modify it. # hide
using MLJ, RDatasets, PrettyPrinting
MLJ.color_off() # hide
@load DecisionTreeClassifier pkg=DecisionTree

carseats = dataset("ISLR", "Carseats")

first(carseats, 3) |> pretty

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
using CSV, DataFrames, Statistics, GLM
df = CSV.read("Data/slides3data.csv"; missingstrings=["NA"])
size(df)

describe(df)

N = size(df,1)
β = [1.65,.4,.06,-.0002]
σ = .4;
df.exper = df.age .- ( 18*(1 .- df.collgrad) .+ 22*df.collgrad )
df.lwsim = β[1] .+ β[2]*df.collgrad .+ β[3]*df.exper .+ β[4]*df.exper.^2 .+ σ*randn(N)
df.lw    = log.(df.wage)

describe(df;cols=[:lw,:lwsim])

β̂ = lm(@formula(lw ~ collgrad + exper + exper^2), df[df.employed.==1,:])
df.elwage = predict(β̂, df) # generates expected log wage for all observations
r2(β̂)                               # reports R2
sqrt(deviance(β̂)/dof_residual(β̂))  # reports root mean squared error

α̂ = glm(@formula(collgrad ~ parent_college + efc), df, Binomial(), LogitLink())
γ̂ = glm(@formula(employed ~ elwage + numkids), df, Binomial(), LogitLink())

df_cfl     = deepcopy(df)
df_cfl.efc = df.efc .- 1         # change value of efc to be $1,000 less
df.basesch = predict(α̂, df)     # generates expected log wage for all observations
df.cflsch  = predict(α̂, df_cfl) # generates expected log wage for all observations

describe(df;cols=[:basesch,:cflsch])


text:
@doc """
## Volume Image
This composite type contains volume image information

#### Fields

* `data`: contains the recorded data
* `x`, `y`, `z`, `t` Arrays containing spatial and time information
* `method` AbstractString of method used to compute tomography
* `info`: additional information in dictionary

#### `processing` Fields
The following standard names are used when saving data to the info dictionary.
* `Regularisation`: Regularisation used in tomography
* `NormalisationConstant`: Value used to normalise image to maximum of 1
* `FileName`: Name of file

""" ->
type VolumeImage
    data::Array{AbstractFloat, 4}
    units::AbstractString
    x::Vector{quantity(AbstractFloat, Meter)}
    y::Vector{quantity(AbstractFloat, Meter)}
    z::Vector{quantity(AbstractFloat, Meter)}
    t::Vector{quantity(AbstractFloat, Second)}
    method::AbstractString
    info::Dict
    coord_system::AbstractString

    function VolumeImage{F<:AbstractFloat, S<:AbstractString}(data::Array{F, 4}, units::S,
                         x::Vector{F}, y::Vector{F}, z::Vector{F}, t::Vector{F},
                         method::S, info::Dict, coord_system::S)

        @assert size(data, 1) == length(x)
        @assert size(data, 2) == length(y)
        @assert size(data, 3) == length(z)
        @assert size(data, 4) == length(t)

        new(data, units, x, y, z, t, method, info, coord_system)
    end

    function VolumeImage{F<:AbstractFloat, S<:AbstractString, Met<:quantity(AbstractFloat, Meter),
                         Sec<:quantity(AbstractFloat, Second)}(data::Array{F, 4}, units::S,
                         x::Vector{Met}, y::Vector{Met}, z::Vector{Met}, t::Vector{Sec},
                         method::S, info::Dict, coord_system::S)

        @assert size(data, 1) == length(x)
        @assert size(data, 2) == length(y)
        @assert size(data, 3) == length(z)
        @assert size(data, 4) == length(t)

        new(data, units, x, y, z, t, method, info, coord_system)
    end

    function VolumeImage{F<:AbstractFloat, S<:AbstractString}(data::Vector{F}, units::S,
                         x::Vector{F}, y::Vector{F}, z::Vector{F}, t::Vector{F},
                         method::S, info::Dict, coord_system::S)

        @assert length(x) == length(data)
        @assert length(y) == length(data)
        @assert length(z) == length(data)
        @assert length(t) == length(data)

        newX = sort(unique(x))
        newY = sort(unique(y))
        newZ = sort(unique(z))
        newT = sort(unique(t))

        L = zeros(typeof(data[1]), length(newX), length(newY), length(newZ), length(newT))

        for idx in 1:length(data)
            idxX = findfirst(newX, x[idx])
            idxY = findfirst(newY, y[idx])
            idxZ = findfirst(newZ, z[idx])
            idxT = findfirst(newT, t[idx])
            L[idxX, idxY, idxZ, idxT] = data[idx]
        end

        new(L, units, newX, newY, newZ, newT, method, info, coord_system)
    end

end

#
# Basic operations
# ----------------
#

import Base: show

function Base.show(io::IO, vi::VolumeImage)

    println(io, "VolumeImage of method $(vi.method) and units $(vi.units)")
    println(io, "  Spanning x: $(vi.x[1]) : $(vi.x[end])")
    println(io, "  Spanning y: $(vi.y[1]) : $(vi.y[end])")
    println(io, "  Spanning z: $(vi.z[1]) : $(vi.z[end])")
    println(io, "  Spanning t: $(vi.t[1]) : $(vi.t[end])")

    if haskey(vi.info, "Regularisation")
        println(io, "  Regularisation: $(vi.info["Regularisation"])")
    end
    if haskey(vi.info, "NormalisationConstant")
        println(io, "  Image has been normalised")
    end
end


text:
using MacroTools

include("scope.jl")
include("block.jl")

lines(s) = split(s, "\n")

codemodule(code, pos) =
  @as _ code scopes(_, pos) filter(s->s.kind==:module, _) map(s->s.name, _) join(_, ".")

precursor(s::AbstractString, i) = join(collect(s)[1:min(i-1, end)])
postcursor(s::AbstractString, i) = join(collect(s)[i:end])

function getblockcursor(code, line, c)
  code, bounds = getblock(code, line)
  code, bounds, cursor(c.line-bounds[1]+1, c.column)
end

getblockcursor(code, cursor) = getblockcursor(code, cursor.line, cursor)

charundercursor(code, cursor) = get(collect(lines(code)[cursor.line]), cursor.column, ' ')

function matchorempty(args...)
  result = match(args...)
  result == nothing ? "" : result.match
end

function getqualifiedname(str::AbstractString, index::Integer)
  pre = precursor(str, index)
  post = postcursor(str, index)

  pre = matchorempty(Regex("(?:$(identifier.pattern)\\.)*(?:$(identifier.pattern))\\.?\$"), pre)

  beginning = pre == "" || last(pre) == '.'
  post = matchorempty(Regex("^$(beginning ? identifier.pattern : identifier_inner.pattern*"*")"), post)

  if beginning && post == ""
    return pre
  else
    return pre * post
  end
end

# could be more efficient
getqualifiedname(str::AbstractString, cursor) = getqualifiedname(lines(str)[cursor.line], cursor.column)

function isdefinition(code::AbstractString)
  try
    code = parse(code)
    return isexpr(code, :function) || (isexpr(code, :(=)) && isexpr(code.args[1], :call))
  catch e
    return false
  end
end

:(a() = b).args[1].head


text:
module ManifoldMeasures

using KeywordCalls,
    LinearAlgebra, Manifolds, MeasureTheory, Random, SpecialFunctions, StatsFuns
using Distributions: Distributions
using Manifolds: ℝ

include("specialfunctions.jl")
include("primitive.jl")
include("stiefel/primitive.jl")
include("stiefel/matrixlangevin.jl")
include("stiefel/matrixbingham.jl")
include("stiefel/macg.jl")

# primitive
export Hausdorff, Haar

# stiefel
export MatrixLangevin, MatrixVonMisesFisher, MatrixBingham, MatrixAngularCentralGaussian

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -3.3773826343471955 <= q <= 0.8255956851677421)

                     add_NL_constraint(m, :(gelu(-0.5817410001457377 + -0.38051045291676777*gelu(-0.7752766980375418 + 0.20887750352529233*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.828517562717189*gelu(-0.31567568822764036 + 0.26872838410166544*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.13996198025628548*gelu(0.6085123540494801 + -0.012952465394535206*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) + gelu(-0.08252400134510385 + 0.3442818429504264*gelu(-0.7752766980375418 + 0.20887750352529233*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.7813134096358496*gelu(-0.31567568822764036 + 0.26872838410166544*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.5888247957791601*gelu(0.6085123540494801 + -0.012952465394535206*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) + gelu(-0.2140493400207517 + 0.6401317985056014*gelu(-0.7752766980375418 + 0.20887750352529233*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.4453350709432602*gelu(-0.31567568822764036 + 0.26872838410166544*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + 0.9148591417614278*gelu(0.6085123540494801 + -0.012952465394535206*gelu(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*gelu(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*gelu(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using BinaryBuilder

platforms = [
  BinaryProvider.Linux(:i686, :glibc),
  BinaryProvider.Linux(:x86_64, :glibc),
  BinaryProvider.Linux(:aarch64, :glibc),
  BinaryProvider.Linux(:armv7l, :glibc),
  BinaryProvider.Linux(:powerpc64le, :glibc),
  BinaryProvider.MacOS(),
  BinaryProvider.Windows(:i686),
  BinaryProvider.Windows(:x86_64)
]

sources = [
  "https://ftp.gnu.org/pub/gnu/libiconv/libiconv-1.15.tar.gz" =>
"ccf536620a45458d26ba83887a983b96827001e92a13847b45e4925cc8913178",
]

script = raw"""
cd $WORKSPACE/srcdir
cd libiconv-1.15/
./configure --prefix=/ --host=$target
make -j80
make install

"""

products = prefix -> [
	LibraryProduct(prefix,"libcharset"),
	LibraryProduct(prefix,"libiconv"),
	ExecutableProduct(prefix,"iconv")
]

autobuild(pwd(), "libiconv", platforms, sources, script, products)


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/functions/nonlinear/univariate.jl
# Defines switching functions for forward and reverse univariate functions.
#############################################################################

exprs = Expr[]
for i = 1:length(univariate_operators)
    op = univariate_operators[i]
    ex = :(return $op(x))
    push!(exprs, ex)
end

function binaryswitch(ids, exprs)
    if length(exprs) <= 3
        out = Expr(:if, Expr(:call, :(==), :operator_id, ids[1]), exprs[1])
        if length(exprs) > 1
            push!(out.args, binaryswitch(ids[2:end], exprs[2:end]))
        end
        return out
    else
        mid = length(exprs) >>> 1
        return Expr(:if, Expr(:call, :(<=), :operator_id, ids[mid]),
            binaryswitch(ids[1:mid], exprs[1:mid]),
            binaryswitch(ids[mid+1:end], exprs[mid+1:end]))
    end
end

switchexpr = binaryswitch(1:length(exprs), exprs)
@eval @inline function eval_univariate_set(operator_id,x::T) where T
    $switchexpr
    error("No match for operator_id")
end

#=
null_interval_rev(y::Interval{T}, x::Interval{T}) where T = y, x
const DEFINED_REVERSES = [:plus_rev, :minus_rev, :inv_rev,
                          :mul_rev, :div_rev, power_rev,
                          :sign_rev, max_rev, min_rev,
                          :sqr_rev, sqrt_rev, abs_rev,
                          :exp_rev, exp2_rev, exp10_rev, expm1_rev,
                          :log_rev, log2_rev, log10_rev, log1p_rev,
                          :sin_rev, cos_rev, tan_rev,
                          :asin_rev, acos_rev, atan_rev,
                          :sinh_rev, cosh_rev, tanh_rev,
                          :asinh_rev, acosh_rev, atanh_rev]
                          =#

# code for defining reverse operators
univariate_operators_rev = [:plus_rev, :minus_rev]
for i in 3:length(univariate_operators)
    string = String(univariate_operators[i])*"_rev"
    push!(univariate_operators_rev, Symbol(string))
end

univariate_reverse_operator_to_id = Dict{Symbol,Int}()
for i = 1:length(univariate_operators)
    univariate_reverse_operator_to_id[univariate_operators_rev[i]] = i
end

exprs_rev = Expr[]
for i = 1:length(univariate_operators_rev)
    op = univariate_operators_rev[i]
    ex = :(return $op(y,x))
    push!(exprs_rev, ex)
end

switchexpr_rev = binaryswitch(1:length(exprs_rev), exprs_rev)
@eval @inline function eval_univariate_set_reverse(operator_id, y::T, x::T) where T
    $switchexpr_rev
    y, x
end

function make_uv_gate(symbol_list)
   n = Expr(:||)
   for (i, sym) in enumerate(symbol_list)
      if i !== 1
         n = :(x == $(univariate_operator_to_id[sym]) || $n)
      else
        n = Expr(:||)
        n = :(x == $(univariate_operator_to_id[sym]))
      end
   end
   n
end

single_tp_ops = [:tan, :cot, :asin, :acos, :atan, :sinh, :tanh,
                 :erf, :erfinv, :erfc, :erfcinv, :erfi, :erfcx, :tand, :cotd,
                 :acosd, :atand, :asind, :asinh, :atanh, :asech]
single_tp_gate = make_uv_gate(single_tp_ops)
@eval @inline single_tp(x::Int) = $single_tp_gate

double_tp_ops = [:sin, :cos, :sind, :cosd, :sech, :dawson]
double_tp_gate = make_uv_gate(double_tp_ops)
@eval @inline double_tp(x::Int) = $double_tp_gate

function make_tp_gate_1(symbol_list, x::Symbol, z::Symbol)
    n = Expr(:block)
    for (i, sym) in enumerate(symbol_list)
        a = :(op == $(univariate_operator_to_id[sym]) && (return (McCormick.$sym)(x), 0.0, 0.0))
        push!(n.args, a)
    end
    push!(n.args, :(error("No operator")))

    nd = Expr(:block)
    for (i, sym) in enumerate(symbol_list)
        sym_kernel = Symbol(String(sym)*"_kernel")
        a = :(op == $(univariate_operator_to_id[sym]) && (return (McCormick.$sym_kernel)(x, z.Intv, tp1, tp2)))
        push!(nd.args, a)
    end
    push!(nd.args, :(error("No operator")))

    qn = quote
            if flag
                $n
            else
                $nd
            end
        end
    qn
end

function make_tp_gate_2(symbol_list, x::Symbol, z::Symbol)
    n = Expr(:block)
    for (i, sym) in enumerate(symbol_list)
        a = :(op == $(univariate_operator_to_id[sym]) && (return (McCormick.$sym)(x), 0.0, 0.0, 0.0, 0.0))
        push!(n.args, a)
    end
    push!(n.args, :(error("No operator")))

    nd = Expr(:block)
    for (i, sym) in enumerate(symbol_list)
        sym_kernel = Symbol(String(sym)*"_kernel")
        a = :(op == $(univariate_operator_to_id[sym]) && (return (McCormick.$sym_kernel)(x, z.Intv, tp1, tp2, tp3, tp4)))
        push!(nd.args, a)
    end
    push!(nd.args, :(error("No operator")))

    qn = quote
            if flag
                $n
            else
                $nd
            end
        end
    qn
end

single_tp_expr = make_tp_gate_1(single_tp_ops, :x, :z)
@eval @inline function single_tp_set(op::Int64, x::MC, z::MC, tp1::Float64,
                               tp2::Float64, flag::Bool)
    $single_tp_expr
end
double_tp_expr = make_tp_gate_2(double_tp_ops, :x, :z)
@eval @inline function double_tp_set(op::Int64, x::MC, z::MC, tp1::Float64,
                               tp2::Float64, tp3::Float64, tp4::Float64, flag::Bool)
    $double_tp_expr
end


text:
using Ipopt

"""
$(SIGNATURES)
Extract important information from the Ipopt output file(s)
"""
function parse_ipopt_summary(file=joinpath(dirname(pathof(TrajectoryOptimization)),"..","logs","ipopt.out"))
    props = Dict()
    obj = Vector{Float64}()
    c_max = Vector{Float64}()
    iter_lines = false  # Flag true if it's parsing the iteration summary lines

    if !isfile(file)
        println("Can't find Ipopt output file.\n Solve statistics will be unavailable.")
        return props
    end

    function stash_prop(ln::String,prop::String,prop_name::Symbol=prop,vartype::Type=Float64)
        if occursin(prop,ln)
            loc = findfirst(prop,ln)
            if vartype <: Real
                val = convert(vartype,parse(Float64,split(ln)[end]))
                props[prop_name] = val
                return true
            end
        end
        return false
    end

    function store_itervals(ln::String)
        if iter_lines
            vals = split(ln)
            if length(vals) == 10 && vals[1] != "iter" && vals[1] != "Restoration" && vals[2] != "iteration"
                push!(obj, parse(Float64,vals[2]))
                push!(c_max, parse(Float64,vals[3]))
            end
        end
    end

    open(file) do f
        for ln in eachline(f)
            stash_prop(ln,"Number of Iterations..",:iterations,Int64) ? iter_lines = false : nothing
            stash_prop(ln,"Total CPU secs in IPOPT (w/o function evaluations)",:self_time)
            stash_prop(ln,"Total CPU secs in NLP function evaluations",:function_time)
            stash_prop(ln,"Number of objective function evaluations",:objective_calls,Int64)
            length(ln) > 0 && split(ln)[1] == "iter" && iter_lines == false ? iter_lines = true : nothing
            store_itervals(ln)
        end
    end
    props[:cost] = obj
    props[:c_max] = c_max
    return props
end

function write_ipopt_options()
    root_dir = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
    if !isdir(joinpath(root_dir,"logs"))
        mkdir(joinpath(root_dir,"logs"))
    end
    outfile=joinpath(root_dir,"logs","ipopt.out")
    optfile=joinpath(root_dir,"ipopt.opt")

    if !isfile(optfile)
        println("Writing Ipopt options file to $optfile...")
        f = open(optfile,"w")
        println(f,"# IPOPT Options for TrajectoryOptimization.jl\n")
        println(f,"# Use Quasi-Newton methods to avoid the need for the Hessian")
        println(f,"hessian_approximation limited-memory\n")
        println(f,"# Output file")
        println(f,"file_print_level 5")
        println(f,"output_file"*" "*"\""*"$(outfile)"*"\"")
        close(f)
    end
end


text:
"""
Operator fileformat

Create the link between the ICAO operator code, the operator name, and
the country of origin

Example:
#;OPERATOR;1;446;20181206;20190102;6478;EAR_P
AAA;ANSETT AUSTRALIA HOLDINGS LTD;ANSETT;AU;AUSTRALIA;19000101;21001231
AAB;ABELAG AVIATION;ABG;BE;BELGIUM;19000101;21001231
AAC;ARMY AIR CORPS;ARMYAIR;GB;UNITED KINGDOM;19000101;21001231
AAD;AERO AVIATION CENTRE LTD.;SUNRISE;CA;CANADA;19000101;21001231
AAE;EXPRESS AIR, INC.  (PHOENIX, AZ);ARIZONA;US;UNITED STATES;19000101;21001231
AAF;AIGLE AZUR;AIGLE AZUR;FR;FRANCE;19000101;21001231
AAG;AIR ATLANTIQUE;ATLANTIC;GB;UNITED KINGDOM;20150101;21001231
"""
module Operator

export read

using Dates

struct OperatorICAO
    name::String
    callsign::String
    countryid::String
    country::String
    startdate::Date
    enddate::Date
    function OperatorICAO(name, callsign, countryid, country, startdate, enddate)
        sd = Date(startdate, "YYYYmmdd")
        ed = Date(enddate, "YYYYmmdd")
        new(name, callsign, countryid, country, sd, ed)
    end
end

function read(file)
    operators = Dict{String, OperatorICAO}()
    for line in eachline(file)
        line_elements = split(line, ';')
        if length(line_elements) == 7
            operators[line_elements[1]] = OperatorICAO(line_elements[2],
            line_elements[3], line_elements[4], line_elements[5],
            line_elements[6], line_elements[7])
        end
    end
    return operators
end

end # module


text:
using Patchwork
using Base.Test

import Patchwork: count, TextNode, Elem

ul(x...) = Elem(:xhtml, :ul, x)
li(x...) = Elem(:xhtml, :li, x)

@testset "Testing Nodes" begin
    @testset "testing counts" begin
        @test count(TextNode("a")) == 0
        @test count(ul()) == 0
        @test count(li("x")) == 1
        @test count(ul(li("x"))) == 2
        @test count(ul(li("x"), li("x"))) == 4
    end
end

@testset "Testing nested node creation" begin
  node = Elem(:div, TextNode("a"))
  @test count(node) == 1
end


text:
module _DefaultSalsaStorage

import ..Salsa
using ..Salsa:
    Runtime, AbstractSalsaStorage, memoized_lookup, get_user_function, collect_trace
using ..Salsa: DependencyKey, DerivedKey, InputKey, RuntimeWithStorage,
    _TopLevelRuntimeWithStorage, _TracingRuntimeWithStorage
using Base.Threads: Atomic, atomic_add!, atomic_sub!
using Base: @lock

import ..Salsa.Debug: @debug_mode, @dbg_log_trace

# --------------------
# TODO:
# - Investigate @nospecialize on user arguments for compiler performance?
#    - We tried this, but saw perf regressions, so we maybe don't understand the
#      specializations that are going on.
# --------------------

const Revision = Int

struct InputValue{T}
    value::T
    changed_at::Revision
    # Allow converting from abitrary values (e.g. creating a Vector{String} from []).
    InputValue{T}(v, changed_at) where {T} = new{T}(v, changed_at)
end
InputValue(v::T, changed_at) where {T} = InputValue{T}(v, changed_at)

_changed_at(v::InputValue)::Revision = v.changed_at

# This struct is mutable so that we can edit the Revisions in-place without having to
# reallocate a new DerivedValue. This is a performance optimization to avoid allocations.
# Since this struct contains a Vector, it will not be isbits, so it will be heap-allocated
# anyway (though this may change in the future, with Julia's stack allocation patch in 1.5;
# consider re-evaluating this decision then).
mutable struct DerivedValue{T}
    value::T
    # A list of all the computations that were accessed directly from this derived function
    # (not all recursive dependencies) when computing this derived value.
    dependencies::Vector{DependencyKey}
    # These Revisions are used to determine whether we need to re-compute this value or not.
    # We need both of them in order to correctly implement the Early-Exit Optimization.
    changed_at::Revision
    verified_at::Revision
end

_changed_at(v::DerivedValue)::Revision = v.changed_at

const InputMapType = IdDict{InputKey,Dict}
const DerivedFunctionMapType = IdDict{Type{<:DerivedKey},Dict}

mutable struct DefaultStorage <: AbstractSalsaStorage
    # The entire Salsa storage is protected by this lock. All accesses and
    # modifications to the storage should be made within this lock.
    lock::Base.ReentrantLock

    # This is bumped every time the storage changes at all, in `set_input!` and
    # `delete_input!`. Whenever this is bumped, we have to check every derived funciton the
    # first time it is called to see if it is still valid.
    # NOTE: This is guaranteed not to change while any derived functions are currently
    # running anywhere on any threads, which we ensure by tracking the count of active
    # derived functions via `derived_functions_active`, below.
    current_revision::Int64

    # We use one big dictionary for the inputs, storing them all together, to reduce
    # allocating a new dictionary for every input.
    # TODO: Do more performance investigation for the tradeoff between sharing a dictionary
    # vs having separate dicts per method. It seems like the current decision is exactly
    # opposite of what would be best: Since DerivedValues are not isbits, they will always
    # be heap allocated, so there's no reason to strongly type their dict. But inputs can
    # be isbits, so it's probably worth specailizing them.
    inputs_map::Dict{InputKey,InputValue}
    derived_function_maps::DerivedFunctionMapType

    # Tracks whether there are any derived functions currently active. It is an error to
    # modify any inputs while derived functions are active, on the current Task or any Task.
    derived_functions_active::Atomic{Int}

    function DefaultStorage()
        new(Base.ReentrantLock(), 0, InputMapType(), DerivedFunctionMapType(), Atomic{Int}(0))
    end
end

const DefaultRuntime = Salsa.Runtime{Salsa.EmptyContext,DefaultStorage}

function Base.show(io::IO, storage::DefaultStorage)
    current_revision = @lock storage.lock begin
        storage.current_revision
    end
    print(io, "Salsa.DefaultStorage($current_revision, ...)")
end

# NOTE: This implements the dynamic behavior for Salsa Components, allowing users to define
# input/derived function dynamically, by attaching new Dicts for them to the storage at
# runtime.
function get_map_for_key(
    storage::DefaultStorage, ::KT, ::Type{RT}
) where {TT, KT<:DerivedKey{<:Any, TT}, RT}
    @lock storage.lock begin
        return get!(storage.derived_function_maps, KT) do
            # PERFORMANCE NOTE: Only construct key inside this do-block to
            # ensure expensive constructor only called once, the first time.

            # TODO: Use the macro's returntype to strongly type the value.
            #       We'll have to generate this function from within the macro, like we used to
            #       in the existing open-source Salsa.
            # NOTE: Except actually after https://github.com/RelationalAI-oss/Salsa.jl/issues/11
            #       maybe we won't do this anymore, and we'll just use one big dictionary!
            Dict{TT,DerivedValue{RT}}()
        # NOTE: Somehow, julia has trouble deducing this return value!
        end::Dict{TT,DerivedValue{RT}}  # This type assertion reduces allocations by 2!!
    end
end
function get_map_for_key(storage::DefaultStorage, ::InputKey)
    # We use one big dictionary for the inputs, storing them all together, so any input key
    # would return the same value here. :)
    return storage.inputs_map
end

function Salsa._previous_output_internal(
    runtime::Salsa._TracingRuntimeWithStorage{DefaultStorage},
    key::DerivedKey,
)
    storage = Salsa.storage(runtime)
    derived_key, args = key, key.args

    previous_output = nothing

    @lock storage.lock begin
        cache = get_map_for_key(storage, derived_key)
        if haskey(cache, args)
            previous_output = getindex(cache, args)
        end
    end

    return previous_output
end

function Salsa._memoized_lookup_internal(
    runtime::Salsa._TracingRuntimeWithStorage{DefaultStorage},
    key::DerivedKey{F,TT},
)  where {F,TT}
    storage = Salsa.storage(runtime)
    try  # For storage.derived_functions_active
        atomic_add!(storage.derived_functions_active, 1)

        local existing_value, value
        found_existing = false
        should_run_user_func = true

        derived_key, args = key, key.args

        user_func = get_user_function(runtime, derived_key)
        # Always just box all the results in the same structure for max type stability,
        # since we don't gain anything by knowing the type here anyway. We just have
        # to rely on julia to deduce the return type at the very end.
        RT = Any

        # NOTE: We currently make no attempts to prevent two Tasks from simultaneously
        # computing the same derived function for the same key. For cheap derived functions
        # this may be a more optimal behavior than the overhead caused by coordination.
        # However for expensive functions, this is not ideal. We take this approach for now
        # because it is simpler.

        # The locking here is a bit involved, to prevent the lock from being held during
        # user computation.
        # We minimally lock only around the accesses to the fields of storage.
        # NOTE: It is okay to release and then reacquire the lock inside this function
        # while computing the result because we are guaranteed that inputs cannot be
        # modified while this function is running, so we do not need to fear concurrency
        # violations, e.g. overwriting newer values with outdated results.
        lock_held::Bool = false
        local cache
        trace = Salsa.trace(runtime) :: Salsa.TraceOfDependencyKeys
        try
            lock(storage.lock)
            lock_held = true

            cache = get_map_for_key(storage, derived_key, RT)

            if haskey(cache, args)
                # Optimization: Skip tracing dependencies when check still_valid on values
                #   - There's no reason to be tracing the Salsa functions during
                #     the `still_valid` check, since we're not going to use them. We _do_
                #     still want to keep the stack trace though for cycle detection and
                #     error messages.
                #   - So we set this toggle on the Trace object itself, to allow us to skip
                #     recording the deps for this phase.
                trace.should_trace = false

                existing_value = getindex(cache, args)::DerivedValue{RT}
                found_existing = true
                unlock(storage.lock)
                lock_held = false

                # storage.current_revision should be monotonically increasing.
                @debug_mode @assert (storage.current_revision >= existing_value.verified_at)

                # NOTE: There is no race condition possible here, despite that the storage
                # isn't locked, because all code that might bump `current_revision`
                # (`set_input!` and `delete_input!`) first asserts that there are no active
                # derived functions (`derived_functions_active == 0`). Meaning that this
                # value will be stable across the lifetime of this function.
                if existing_value.verified_at == storage.current_revision
                    value = existing_value
                    should_run_user_func = false
                    # NOTE: still_valid() will recursively call memoized_lookup, potentially
                    #       recomputing all our recursive dependencies.
                elseif still_valid(runtime, existing_value)
                    # Update the verified_at field, but otherwise use the existing value.
                    # NOTE: As above, current_revision is safe to read during this function
                    # without a lock, due to asserts on derived_functions_active.
                    existing_value.verified_at = storage.current_revision
                    value = existing_value
                    should_run_user_func = false
                end
            end
        finally
            if lock_held
                unlock(storage.lock)
                lock_held = false
            end
            trace.should_trace = true
        end

        # At this point (value == nothing) if (and only if) the args are not
        # in the cache, OR if they are in the cache, but they are no longer valid.
        if should_run_user_func
            @dbg_log_trace @info "invoking $key"
            if found_existing
                # Dependency array swap Optimization:
                #   If we've already got an `existing_value` object here, we can avoid an
                #   allocation and a copy by _swapping_ the `trace`'s `ordered_dependencies`
                #   with `existing_value.dependencies`, so that the deps are written
                #   in-place directly into their final destination! :)
                trace = Salsa.get_trace(runtime.immediate_dependencies_id)
                # Temporarily swap the dependency vectors while running user_func so the
                # deps are recorded in-place. Note that we must swap them back at the end.
                existing_value.dependencies, trace.ordered_deps =
                    trace.ordered_deps, existing_value.dependencies
                try
                    v = user_func(runtime, key.args...)
                finally
                    # Swap back the dependency vectors so the vector isn't modified by
                    # future traces.
                    existing_value.dependencies, trace.ordered_deps =
                        trace.ordered_deps, existing_value.dependencies
                end
            else
                v = user_func(runtime, key.args...)
            end
            # NOTE: We use `isequal` for the Early Exit Optimization, since values are
            # required to be purely immutable (but not necessarily julia `immutable
            # structs`).
            @dbg_log_trace @info "Returning from $key."
            if found_existing && isequal(existing_value.value, v)
                # Early Exit Optimization Part 2: (for Part 1 see `set_input!`, below).
                # If a derived function computes the exact same value, we can terminate
                # early and "backdate" the changed_at field to say this value has _not_
                # changed.
                # NOTE: As above, current_revision is safe to read during this function
                #       without a lock, due to asserts on derived_functions_active.
                existing_value.verified_at = storage.current_revision
                # Note that just because it computed the same value, it doesn't mean it
                # computed it in the same way, so we need to update the list of
                # dependencies as well.
                # HOWEVER, we can actually skip this, thanks to the swap-optimization above.
                # existing_value.dependencies = collect_trace(runtime)

                # We keep the old computed `.value` rather than the new value to help catch
                # bugs with users' over-permissive `isequal()` functions earlier.
                value = existing_value
            elseif found_existing
                # Reuse as much of the existing_value's structure to avoid allocations
                existing_value.value = v
                # We skip this thanks to the swap-optimization, above.
                # existing_value.dependencies = collect_trace(runtime)
                existing_value.changed_at = storage.current_revision
                existing_value.verified_at = storage.current_revision
                value = existing_value
            else
                @dbg_log_trace @info "Computed new derived value for $key."
                # The user function computed a new value, which we must now store.
                # NOTE: We set the computed RT here, which might be more abstract than the
                # actual type of the value, `v`.
                # The other option is to change the dicts to be Dict{K, DerivedValue{<:RT}},
                # but that causes extra allocations.
                value = DerivedValue{RT}(
                    v,
                    collect_trace(runtime),
                    storage.current_revision,
                    storage.current_revision,
                )
                @lock storage.lock begin
                    cache[args] = value
                end
            end # existing_value
        end # if value === nothing

        return value
    finally
        atomic_sub!(storage.derived_functions_active, 1)
    end
end # _memoized_lookup_internal

function Salsa._unwrap_salsa_value(
    runtime::RuntimeWithStorage{DefaultStorage},
    # Note: the presence of a `where T` on this function causes allocations.
    v #=::Union{DerivedValue,InputValue} =#
)
    return v.value
end

# A `value` is still valid if none of its dependencies have changed.
function still_valid(runtime, value)
    for depkey in value.dependencies
        dep_changed_at = key_changed_at(runtime, depkey)
        if dep_changed_at > value.verified_at
            return false
        end
    end # for
    return true
end

function key_changed_at(runtime, key::DependencyKey)
    return _changed_at(memoized_lookup(runtime, key))
end

# =============================================================================

# --- Inputs --------------------------------------------------------------------------

function Salsa._unmemoized_input_lookup_internal(
    runtime::Salsa._TracingRuntimeWithStorage{DefaultStorage},
    key::InputKey,
)
    storage = Salsa.storage(runtime)
    cache = get_map_for_key(storage, key)
    @lock storage.lock begin
        return get(cache, key, nothing)
    end
end

function Salsa._memoized_lookup_internal(
    runtime::Salsa._TracingRuntimeWithStorage{DefaultStorage},
    key::InputKey,
)
    storage = Salsa.storage(runtime)
    cache = get_map_for_key(storage, key)
    @lock storage.lock begin
        return cache[key]
    end
end

function Salsa.set_input!(
    runtime::_TopLevelRuntimeWithStorage{DefaultStorage},
    key::InputKey,
    value,
)
    storage = Salsa.storage(runtime)

    @lock storage.lock begin
        cache = get_map_for_key(storage, key)

        if haskey(cache, key) && _value_isequal_to_cached(cache[key], value)
            # Early Exit Optimization Part 1: Don't dirty anything if setting exactly the
            # same value for an input.
            return
        end

        # It is an error to modify any inputs while derived functions are active, even
        # concurrently on other threads.
        @assert storage.derived_functions_active[] == 0

        @dbg_log_trace @info "Setting input $key => $value"
        storage.current_revision += 1

        cache[key] = InputValue(value, storage.current_revision)
        return nothing
    end
end
# This function barrier exists to allow specializing the `.value` on the type of
# the cached InputValue. (It doesn't seem to have any effect on performance though?)
function _value_isequal_to_cached(cached::InputValue, value)
    # NOTE: We use `isequal` for the Early Exit Optimization, since we compare values _by
    # value_, not by identity. (That is, `[] == []`, despite `[] !== []`.) And we prefer
    # isequal over `==` since we want to preserve float diffs, just like a Dict would.
    return isequal(cached.value, value)
end

function Salsa.delete_input!(
    runtime::_TopLevelRuntimeWithStorage{DefaultStorage},
    key::InputKey,
)
    @dbg_log_trace @info "Deleting input $key"
    storage = Salsa.storage(runtime)
    cache = get_map_for_key(storage, key)

    @lock storage.lock begin
        # It is an error to modify any inputs while derived functions are active, even
        # concurrently on other threads.
        @assert storage.derived_functions_active[] == 0

        storage.current_revision += 1
        delete!(cache, key)
        return nothing
    end
end

function Salsa.new_epoch!(runtime::Salsa.RuntimeWithStorage{DefaultStorage})
end

#= Debug Inspection =#

# Returns two dicts (inputs: K=>V and derived_vals: K=>(V, K[deps...]))
function Salsa.Inspect.all_inputs_and_derived_vals(st::DefaultStorage)
    inputs = st.inputs_map
    derived_vals = Dict(
        key_t(k) => (v.value, v.dependencies)
        for (key_t,m) in st.derived_function_maps
        for (k,v) in m
    )
    return inputs, derived_vals
end

end  # module


text:
using Redpitaya
using PyPlot

rp = RedPitaya("192.168.1.26")

freq = 25000
dec = optimalDecimation(dec,freq)

numPeriods = 4
freqR = roundFreq(rp,dec,freq)
numSampPerPeriod = numSamplesPerPeriod(rp,dec,freqR)
numSamp = numSampPerPeriod*numPeriods

println("Frequency = $freqR Hz")
println("Number Sampling Points per Period: $numSampPerPeriod")

# start sending
send(rp,"GEN:RST")
sendAnalogSignal(rp,1,"SINE",freqR,0.4)

# receive data
u1 = receiveAnalogSignal(rp, 1, 0, numSamp, dec=dec, delay=0.8)

figure(1)
clf()
subplot(2,1,1)
plot(u1)
subplot(2,1,2)
semilogy(abs(rfft(u1))[1:numPeriods*4],"o-b",lw=2)

#Lets check if the data is periodic
figure(2)
clf()
plot(1:length(u1),u1,"x-r",lw=2)
plot(length(u1)+1:2*length(u1),u1,"x-b",lw=2)


text:
function replace_with_arrays(t::T, d::Dict) where {T <: SymbolicObject}
    sc = get_scalars(t)
    pydict = Dict([ k.__pyobject__ => v for (k,v) in d])
    pyexp = t.__pyobject__
    for s in sc
        pyexp = pyexp.subs(s,replace_with_arrays(scalar_exprs[s],d))
    end
    return sympy_type_convert(pyexp)
end
function replace_with_arrays(t::T, d::Dict) where {T <: Tensor}
    sc = get_scalars(t)
    pydict = Dict([ k.__pyobject__ => v for (k,v) in d])
    pyexp = t.__pyobject__
    for s in sc
        pyexp = pyexp.subs(s,replace_with_arrays(scalar_exprs[s],d))
    end
    pyexp = pyexp.replace_with_arrays(pydict)
    return sympy_type_convert(pyexp)
end

#function replace_with_arrays(t::T, d::Dict,TIT::TensorIndexType) where {T <: Tensor}
#    sc = get_scalars(t)
#    for s in sc
#        d[s] = replace_with_arrays(scalar_exprs[s],d)
#    end
#    heads = [ k.head  for (k,v) in d if typeof(k)==IndexedTensor ]
#    headkeys = Dict(k.head=> k  for (k,v) in d if typeof(k)==IndexedTensor )
#    ans = []
#    for t in terms(t)
#        if typeof(t) == IndexedTensor
#            if t.head in heads
#                if headkeys[t.head].is_up == t.is_up
#                    push!(ans,heads[t.head])
#                elseif TIT.metric in heads
#                    if all( x.is_up for headkeys[t.head].indices
#                        # if indices don't match, use metric to convert them
#                        push!(ans,d[headkeys[TIT.metric]]*)
#
#                else
#                    error("Couldn't replace $t, only had $(headkeys[t.head]) and no metric to convert them")
#                end
#            elseif t.head == TIT.delta
#                push!(ans,Matrix{I,dim,dim})
#            end
#        elseif typeof(t) == Sym
#            for s in sc
#                t = subs(t,s => d[s])
#            end
#            push!(ans,t)
#        else
#            push!(replace_with_arrays(t,d))
#        end
#    end
#    return prod(ans)
#end


text:
function dfs_test()
    A = load_matrix_network("dfs_example")
    dfs(MatrixNetwork(A),1)

    n = 10
    O = ones(Int64,n-1)
    Z = zeros(Int64,n)
    A = Tridiagonal(O,Z,O)
    B = sparse(full(A))
    (d,dt,ft,pred) = dfs(B,1)
    if !(d == collect(0:n-1) && dt == collect(0:n-1) && ft == collect(2*n-1:-1:n) &&
                                    pred == collect(0:n-1))
        error("dfs test failed")
    end
    return true
end