text:
module AcceleratedDubins

using Optim

export
    path_len, path_time, sample_path, speed_by_radius, speed_profile, fastest_path, retrieve_path, optimized_path, possible_tuples, radii_samples_lin, radii_samples_exp

###############################################################################
############################ TYPE DEFINITIONS #################################
###############################################################################

@enum DubinsPathTypeR2 LSL LSR RSL RSR
@enum SegmentType L_SEG S_SEG R_SEG

DIRDATA = Dict{Int,Vector{SegmentType}}(
    Int(LSL) => [L_SEG, S_SEG, L_SEG],
    Int(LSR) => [L_SEG, S_SEG, R_SEG],
    Int(RSL) => [R_SEG, S_SEG, L_SEG],
    Int(RSR) => [R_SEG, S_SEG, R_SEG],
)

const PATH_OK = 0       # no error, path found
const PATH_ERROR = 1    # path not possible

"""
The data structure that holds the full dubins path with multiple radii.

# Data fields
- `origin::Vector{Float64}`: the initial configuration [x, y, a], angle in radians
- `lengths::Vector{Float64}`: lengths of each segment
- `r::Vector{Float64}`: turn radii of each segment, where radius of straight path is Inf
- `type::DubinsPathTypeR2`: the Dubins path type given by the @enum DubinsPathTypeR2
"""
mutable struct DubinsPathR2
   origin::Vector{Float64}
   lengths::Vector{Float64}
   r::Vector{Float64}
   type::DubinsPathTypeR2
end

DubinsPathR2() = DubinsPathR2(zeros(3), zeros(3), [1, Inf, 1], LSL)

###############################################################################
########################### HELPER FUNCTIONS ##################################
###############################################################################

"""
    radii_samples_lin(r_min::Float64, r_max::Float64, num_samples::Int64, precision::Int64=100)
Calculate radii samples with linear spacing.
...
# Arguments
- `r_min::Float64`: minimum turning radius
- `r_max::Float64`: maximum turning radius
- `num_samples::Int64`: number of samples
...
return - samples::Vector{Float64}

See also: [`radii_samples_exp`](@ref)
"""
function radii_samples_lin(r_min::Float64, r_max::Float64, num_samples::Int64, precision::Int64=100)
    if precision == 100 # 16 should be reasonable maximum precision for Float64
                        # thus the default precision for rounding is this high
        return collect(r_min:(r_max-r_min)/(num_samples-1):r_max)
    end
    return [round(i, digits=precision) for i in r_min:(r_max-r_min)/(num_samples-1):r_max]
end

"""
    radii_samples_exp(r_min::Float64, r_max::Float64, num_samples::Int64, precision::Int64=100)
Calculate radii samples with exponential spacing.
...
# Arguments
- `r_min::Float64`: minimum turning radius
- `r_max::Float64`: maximum turning radius
- `num_samples::Int64`: number of samples
- `precision::Int64`: for rounding the values. Default: Inf (no rounding)
...
return - samples::Vector{Float64}

See also: [`radii_samples_lin`](@ref)
"""
function radii_samples_exp(r_min::Float64, r_max::Float64, num_samples::Int64, precision::Int64=100)
    if precision == 100
        return [exp(i) for i in log(r_min):(log(r_max) - log(r_min))/(num_samples-1):log(r_max)]
    end
    return [round(exp(i), digits=precision) for i in log(r_min):(log(r_max) - log(r_min))/(num_samples-1):log(r_max)]
end

" Compute DubinsR2 path length. Argument: `path::DubinsPathR2, return: `len::Float64`"
function path_len(path)
    return (path == nothing) ? Inf : sum(path.lengths)
end

"""
    is_path_faster(path_lengths, radii::Vector{Float64}, best_time::Float64, params)

Comparing function to use between two maneuvres.
...
# Arguments
- `path_lengths::Vector{Float64}`: vector of segment lengths
- `radii::Vector{Float64}`: vector of segment radii
- `best_time::Float64`: time to compare with
- `params::Vector{Float64}`: Vector of speed parameters: [velocity in minimal radius, maximal vehicle velocity, acceleration, deceleration]
...
return - ans::Boolean, time::Float64
"""
function is_path_faster(path_lengths, radii::Vector{Float64}, best_time::Float64, params)
    if path_lengths == nothing
        return false, nothing
    end
    path = DubinsPathR2()
    path.type = LSL # for time computing, only posible maneuvres currently are CSC
    path.lengths = path_lengths
    path.r = radii
    time = path_time(path, params)
    if time < best_time
        return true, time
    end
    return false, time
end

"""
    speed_by_radius(radius::Number, vr::Number, vs::Number)
Caculation using plane turn

# Arguments
- `radius::Number`: Radius to get speed for
- `vr::Number`: Maximum vehicle speed on radius 1
- `vs::Number`: Maximum vehicle speed on straight segment

# Return
- `v::Float64`: Maximum speed for input radius
"""
function speed_by_radius(radius::Number, vr::Number, vs::Number)
    return min(sqrt(radius)*vr, vs)
end

"""
    speed_by_radius(radius::Number)
Caculation using plane turn

# Arguments
- `radius::Number`: Radius to get speed for

# Return
- `v::Float64`: Maximum speed for input radius
"""
function speed_by_radius(radius::Number)
    v = sqrt(radius*tan(pi/3)*9.80665) #g
    return v
end

"""
    path_time(path::DubinsPathR2, params::Vector{Number})
Get path time of single path.

# Arguments
- `path::DubinsPathR2`: Path to compute time from
- `params::Vector{Float64}`: Vector of speed parameters: [velocity in minimal radius, maximal vehicle velocity, acceleration, deceleration]

# Return
- `time::Float64`: Path time
"""
function path_time(path::DubinsPathR2, params::Vector{Float64})
    if path == nothing
        return Inf
    end
    t, _ = speed_profile(path, params[1], params[2], params[3], params[4])
    return t[end]
end

"""
    point_by_angle(origin::Vector{Float64}, ang::Number, len::Number)
Get point that is in certain distance by certain angle from origin point

# Arguments
- `origin::Vector{Float64}`: starting angle in radians from point [0, 0]
- `ang::Number`: starting angle in radians from point [0, 0]
- `len::Number`: starting angle in radians from point [0, 0]

# Return
- `[x2,y2]::Vector{Float64}`: desired point
"""
function point_by_angle(origin::Vector{Float64}, ang::Number, len::Number)
    x, y = origin[1], origin[2]
    y2 = y + len * sin(ang)
    x2 = x + len * cos(ang)
    return [x2, y2]
end

"""
    function possible_tuples(radii::Vector{Float64})
Get all possible combinations of input radii

# Arguments
- `radii::Vector{Float64}`: array of possible radii

# Return
- `tup::Vector{Float64}`: pairs of radii
"""
function possible_tuples(radii::Vector{Float64})
    tup::Vector{Vector{Float64}} = []
    for i in 1  :length(radii)
        for j in i:length(radii)
            push!(tup, [radii[i], radii[j]])
            if i != j
                push!(tup, [radii[j], radii[i]])
            end
        end
    end
    return tup
end

"""
    rotate_points(start::Vector{Float64}, stop::Vector{Float64})
Rotate given input points to format start -> [0, 0, alfa], stop -> [dist, 0, beta]

# Arguments
- `start::Vector{Float64}`: start point of the vehicle, [x1, y1, theta1]
- `stop::Vector{Float64}`: final point of the vehicle, [x2, y2, theta2]

# Return
- `alfa::Float64`: new starting angle
- `beta::Float64`: new ending angle
- `dist::Float64`: distance between start and stopex
"""
function rotate_points(start::Vector{Float64}, stop::Vector{Float64})
    diff = stop-start
    global_rotation = atan(diff[2], diff[1])
    alfa = start[3] - global_rotation
    beta = stop[3] - global_rotation
    dist = sqrt(diff[1]^2 + diff[2]^2)
    return alfa, beta, dist
end

###############################################################################
######################### MAIN PATH FUNCTIONS #################################
###############################################################################

"""
    retrieve_path(start::Float64, stop::Float64, r1::Float64, r2::Float64, params::Vector{Float64})
Function to return best DubinsR2 path with first and second radius as specified.

# Arguments
- `start::Vector{Float64}`: starting angle in radians from point [0, 0]
- `stop::Vector{Float64}`: starting angle in radians from point [0, 0]
- `radii::Vector{Float64}`: [radius of starting arc, radius of ending arc], must be of length 2
- `params::Array{Float64}`: Vector of speed parameters: [velocity in minimal radius, maximal vehicle velocity, acceleration, deceleration]

# Return
- `path::DubinsPathR2`: fastest DubinsR2 path
- `err::const`: error code of the path calculation
"""
function retrieve_path(start::Vector{Float64}, stop::Vector{Float64}, radii::Vector{Float64}, params::Vector{Float64})
    @assert length(radii) == 2
    alfa, beta, dist = rotate_points(start, stop)

    path, err = best_maneuver(alfa, beta, dist, radii[1], radii[2], params)

    if (err == PATH_ERROR)
        return nothing, PATH_ERROR
    end
    path.origin = start
    return path, PATH_OK
end

"""
    fastest_path(start::Vector{Float64}, stop::Vector{Float64}, radii::Vector{Float64}, params::Array{Float64})
based on the input radii, tests each combination and select best DubinsR2 path.

# Arguments
- `start::Vector{Float64}`: starting angle in radians from point [0, 0]
- `stop::Vector{Float64}`: starting angle in radians from point [0, 0]
- `radii::Vector{Float64}`: radius of starting arc
- `params::Array{Float64}`: Vector of speed parameters: [velocity in minimal radius, maximal vehicle velocity, acceleration, deceleration]

# Return
- `path::DubinsPathR2`: fastest DubinsR2 path
- `err::const`: error code of the path calculation
"""
function fastest_path(start::Vector{Float64}, stop::Vector{Float64}, radii::Vector{Float64}, params::Array{Float64})
    alfa, beta, dist = rotate_points(start, stop)

    # get all possible combinations
    radii_pairs = possible_tuples(radii)

    best_time = Inf
    best_len = Inf
    best_path = nothing
    ret_err = PATH_ERROR
    for pair in radii_pairs
        path, err = best_maneuver(alfa, beta, dist, pair[1], pair[2], params)
        time = path_time(path, params)
        if time < best_time || (isapprox(time, best_time, atol=1e-8) && path_len(path) < best_len)
            best_time = time
            best_len = path_len(path)
            best_path = path
            ret_err = err
        end
    end

    best_path.origin = start
    return best_path, ret_err
end

"""
    optimized_path(start::Vector{Float64}, stop::Vector{Float64}, r_min::Float64, r_max::Float64, params::Array{Float64})
based on the input radii, tests each combination and select best DubinsR2 path.

# Arguments
- `start::Vector{Float64}`: starting angle in radians from point [0, 0]
- `stop::Vector{Float64}`: starting angle in radians from point [0, 0]
- `radii::Vector{Float64}`: vector containing minimum and maximum vehicle radius
- `params::Array{Float64}`: Vector of speed parameters: [velocity in minimal radius, maximal vehicle velocity, acceleration, deceleration]

# Return
- `path::DubinsPathR2`: fastest DubinsR2 path
- `err::const`: error code of the path calculation
"""
function optimized_path(start::Vector{Float64}, stop::Vector{Float64}, radii::Vector{Float64}, params::Array{Float64})
    alfa, beta, dist = rotate_points(start, stop)
    r_min = minimum(radii)
    r_max = maximum(radii)

    function optim_radii(v)
        path, _ = best_maneuver(alfa, beta, dist, v[1], v[2], params)
        return path_time(path, params)
    end

    v = [0.5*(r_min+r_max) for i in 1:2]

    v2 = [r_min+1e-12 for i in 1:2]
    lb = [r_min for i in 1:2] # lowerbound vector, smaller than v,
                                    # otherwise Optim would set it between ub and lb
    ub = [r_max for i in 1:2] # upperbound vector

    max_iter= Int(10e5)

    res1 = optimize(optim_radii, lb, ub, v, SAMIN(verbosity=0)) # returns the best path
    res2 = optimize(optim_radii, lb, ub, v2, SAMIN(verbosity=0)) # returns the best path

    result = res2
    if Optim.minimum(res1) < Optim.minimum(res2)
        result = res1
    end

    if Optim.minimum(result) == Inf
        return nothing, PATH_ERROR
    end

    best_r = Optim.minimizer(result) # get the optimal value
    @assert best_r[1] <= r_max
    @assert best_r[1] >= r_min
    @assert best_r[2] <= r_max
    @assert best_r[2] >= r_min
    path, err = best_maneuver(alfa, beta, dist, best_r[1], best_r[2], params)

    path.origin = start
    return path, err
end

"""
    best_maneuver(alfa::Float64, beta::Float64, dist::Float64, r1::Float64, r2::Float64, params::Array{Float64})
Function to return best DubinsR2 path according to is_path_better argument

# Arguments
- `alfa::Vector{Float64}`: starting angle in radians from point [0, 0]
- `beta::Float64`: ending angle in radians from point [dist, 0]
- `dist::Float64`: distance between starting and ending point
- `r1::Float64`: radius of starting arc
- `r2::Float64`: radius of ending arc
- `params::Array{Float64}`: Vector of speed parameters: [velocity in minimal radius, maximal vehicle velocity, acceleration, deceleration]

# Return
- `path::DubinsPathR2`: path
- `err`: result of operation
"""
function best_maneuver(alfa::Float64, beta::Float64, dist::Float64, r1::Float64, r2::Float64, params::Array{Float64})
    path = DubinsPathR2()
    path.origin = [0, 0, alfa]
    best = Inf

    "shortcut function to set all the path variables"
    function pset(flengths, fbest, fr1, fr2, fr3, ftype)
        path.lengths = flengths
        best = fbest
        path.r[1] = fr1
        path.r[2] = fr2
        path.r[3] = fr3
        path.type = DubinsPathTypeR2(ftype)
    end

    _, path_lengths = RSR_path(alfa, beta, dist, r1, r2)
    ans, calc = is_path_faster(path_lengths, [r1, Inf, r2], best, params)
    if ans
        pset(path_lengths, calc, r1, Inf, r2, 3)
    end

    _, path_lengths = LSL_path(alfa, beta, dist, r1, r2)
    ans, calc = is_path_faster(path_lengths, [r1, Inf, r2], best, params)
    if ans
        pset(path_lengths, calc, r1, Inf, r2, 0)
    end

    _, path_lengths = RSL_path(alfa, beta, dist, r1, r2)
    ans, calc = is_path_faster(path_lengths, [r1, Inf, r2], best, params)
    if ans
        pset(path_lengths, calc, r1, Inf, r2, 2)
    end

    _, path_lengths = LSR_path(alfa, beta, dist, r1, r2)
    ans, calc = is_path_faster(path_lengths, [r1, Inf, r2], best, params)
    if ans
        pset(path_lengths, calc, r1, Inf, r2, 1)
    end

    if best == Inf
        return nothing, PATH_ERROR
    end
    return path, PATH_OK
end

###############################################################################
#################### PATH SAMPLING FOR PLOTTING ###############################
###############################################################################

"""
    sample_arc(start::Vector{Float64}, a1::Float64, a2::Float64, r::Float64, part::Float64, dir::Int)
Sample single arc: first three arguments similar to TikZ definition
Center of the arc is at coordinates [start[1] + r*cos(a1+pi), start[2] + r*sin(a1+pi)]

# Arguments
- `start::Vector{Float64}`: vector of starting point: [x, y, _] or [x, y]
- `a1::Float64`: starting angle in radians
- `a2::Float64`: ending angle in radians
- `r::Float64`: radius of the arc
- `part::Float64`: how smoothly to sample the arc, smaller -> finer
- `dir::Int`: direction of the arc, 1 -> clockwise, -1 -> counter-clockwise

# Return
- `confx::Vector{Float64}`: x points of the path to plot
- `confy::Vector{Float64}`: y points of the path to plot
"""
function sample_arc(start::Vector{Float64}, a1::Float64, a2::Float64, r::Float64, part::Float64, dir::Int)
    ang1 = mod2pi(a1)
    ang2 = mod2pi(a2)
    len = part/(2*pi)
    S = (start[1]-r*cos(ang1), start[2]-r*sin(ang1))
    if isapprox(ang1, ang2, atol = 1e-500)
        x = S[1] + r * cos(ang1)
        y = S[2] + r * sin(ang1)
        return [x], [y]
    end

    confx::Vector{Float64} = []
    confy::Vector{Float64} = []

    if dir == 1
        len = - len
        if ang1 < ang2
            ang2 -= 2*pi
        end
    elseif dir == -1
        if ang2 < ang1
            ang2 += 2*pi
        end
    end
    for i in ang1:len:ang2
        x = S[1] + r * cos(i)
        y = S[2] + r * sin(i)
        push!(confx, x)
        push!(confy, y)
    end

    return confx, confy
end

"""
    sample_path(path::DubinsPathR2, resolution::Float64 = 0.01)
Sample single path for plotting

# Arguments
- `path::DubinsPathR2`: path to sample
- `resolution::Float64`: how smoothly to sample the arc, smaller -> finer

# Return
- `x::Vector{Float64}`: x points of the path to plot
- `y::Vector{Float64}`: y points of the path to plot
"""
function sample_path(path::DubinsPathR2, resolution::Float64 = 0.01)
    if path == nothing
        return [], []
    elseif Integer(path.type) == Integer(RSR)
        ang1 = path.origin[3]+pi/2
        ang2 = ang1 - path.lengths[1]/path.r[1]
        confx, confy = sample_arc(path.origin, ang1, ang2, path.r[1], resolution, 1)
        A = [confx[end], confy[end]]
        B = point_by_angle(A, ang2-pi/2, path.lengths[2])
        push!(confx, B[1])
        push!(confy, B[2])
        ang3 = ang2 - path.lengths[3]/path.r[3]
        confx2, confy2 = sample_arc(B, ang2, ang3, path.r[3], resolution, 1)
        append!(confx, confx2)
        append!(confy, confy2)
        return confx, confy
    end
    if Integer(path.type) == Integer(LSL)
        ang1 = path.origin[3]-pi/2
        ang2 = ang1 + path.lengths[1]/path.r[1]
        confx, confy = sample_arc(path.origin, ang1, ang2, path.r[1], resolution, -1)
        A = [confx[end], confy[end]]
        B = point_by_angle(A, (ang2+pi/2), path.lengths[2])
        push!(confx, B[1])
        push!(confy, B[2])
        ang3 = ang2 + path.lengths[3]/path.r[3]
        confx2, confy2 = sample_arc(B, ang2, ang3, path.r[3], resolution, -1)
        append!(confx, confx2)
        append!(confy, confy2)
        return confx, confy
    end
    if Integer(path.type) == Integer(LSR)
        ang1 = path.origin[3]-pi/2
        ang2 = ang1 + path.lengths[1]/path.r[1]
        confx, confy = sample_arc(path.origin, ang1, ang2, path.r[1], resolution, -1)
        A = [confx[end], confy[end]]
        B = point_by_angle(A, ang2+(pi/2), path.lengths[2])
        push!(confx, B[1])
        push!(confy, B[2])
        ang3 = ang2 - pi - path.lengths[3]/path.r[3]
        confx2, confy2 = sample_arc(B, mod2pi(ang2+pi), ang3, path.r[3], resolution, 1)
        append!(confx, confx2)
        append!(confy, confy2)
        return confx, confy
    end
    if Integer(path.type) == Integer(RSL)
        ang1 = path.origin[3]+pi/2
        ang2 = ang1 - path.lengths[1]/path.r[1]
        confx, confy = sample_arc(path.origin, ang1, ang2, path.r[1], resolution, 1)
        A = [confx[end], confy[end]]
        B = point_by_angle(A, ang2-(pi/2), path.lengths[2])
        push!(confx, B[1])
        push!(confy, B[2])
        ang3 = (ang2 + pi) + path.lengths[3]/path.r[3]
        confx2, confy2 = sample_arc(B, mod2pi(ang2+pi), ang3, path.r[3], resolution, -1)
        append!(confx, confx2)
        append!(confy, confy2)
        return confx, confy
    end
    return [], []
end

###############################################################################
###################### PATH TYPES COMPUTATION #################################
###############################################################################

function RSR_path(alfa::Float64, beta::Float64, dist::Number, radius1::Number, radius2::Number)
    S1 = (radius1*sin(alfa), -radius1*cos(alfa))
    S2 = (dist+radius2*sin(beta), -radius2*cos(beta))

    dx = S2[1] - S1[1]
    dy = S2[2] - S1[2]
    d = sqrt(dx^2 + dy^2)
    if (d^2 - (radius1 - radius2)^2) < 0
        return PATH_ERROR, nothing
    end
    l = sqrt(d^2 - (radius1 - radius2)^2)

    gamma = atan(dy/dx)
    xi = acos((radius1-radius2)/d)

    mu = mod2pi(gamma+xi)
    if (S2[1] < S1[1])
        mu += pi
    end

    out = Vector{Float64}(undef,3)
    out[1] = mod2pi(alfa + pi/2 - mu) * radius1
    out[2] = l
    out[3] = mod2pi(2*pi - (beta + pi/2 - mu)) * radius2

    return PATH_OK, out
end

function LSL_path(alfa, beta, dist, radius1, radius2)
    S1 = (-radius1*sin(alfa), radius1*cos(alfa))
    S2 = (dist-radius2*sin(beta), radius2*cos(beta))
    dx = S2[1] - S1[1]
    dy = S2[2] - S1[2]
    d = sqrt(dx^2 + dy^2)
    if (d^2 - (radius1 - radius2)^2) < 0
        return PATH_ERROR, nothing
    end
    l = sqrt(d^2 - (radius1 - radius2)^2)

    gamma = atan(dy/dx)
    xi = asin((radius1-radius2)/d)

    mu = mod2pi(gamma+xi + 3pi/2)
    if (S2[1] < S1[1])
        mu += pi
    end

    out = Vector{Float64}(undef,3)
    out[1] = mod2pi((mu - (alfa - pi/2))) * radius1
    out[2] = l
    out[3] = mod2pi((beta - pi/2 - mu)) * radius2

    return PATH_OK, out
end

function RSL_path(alfa, beta, dist, radius1, radius2)

    S1 = (radius1*sin(alfa), -radius1*cos(alfa))
    S2 = (dist-radius2*sin(beta), radius2*cos(beta))
    dx = S2[1] - S1[1]
    dy = S2[2] - S1[2]
    d = sqrt(dx^2 + dy^2)
    if (d^2 - (radius1 + radius2)^2) < 0
        return PATH_ERROR, nothing
    end
    l = sqrt(d^2 - (radius1 + radius2)^2)

    gamma = atan(dy/dx)
    xi = asin((radius1+radius2)/d)

    mu = mod2pi(gamma-xi + pi/2)

    out = Vector{Float64}(undef,3)
    out[1] = mod2pi(((alfa + pi/2) - mu)) * radius1
    out[2] = l
    out[3] = mod2pi(beta - pi/2 - (mu + pi)) * radius2

    return PATH_OK, out
end

function LSR_path(alfa, beta, dist, radius1, radius2)

    S1 = (-radius1*sin(alfa), radius1*cos(alfa))
    S2 = (dist+radius2*sin(beta), -radius2*cos(beta))

    dx = S2[1] - S1[1]
    dy = S2[2] - S1[2]
    d = sqrt(dx^2 + dy^2)
    if (d^2 - (radius1 + radius2)^2) < 0
        return PATH_ERROR, nothing
    end
    l = sqrt(d^2 - (radius1 + radius2)^2)

    gamma = atan(dy/dx)
    xi = asin((radius1+radius2)/d)

    mu = gamma + xi - pi/2

    out = Vector{Float64}(undef,3)
    out[1] = mod2pi(mu - (alfa - pi/2)) * radius1
    out[2] = l
    out[3] = mod2pi((mu + pi) - (beta + pi/2)) * radius2

    if isapprox(mod2pi((mu + pi) - (beta + pi/2)), 2pi, atol=1e-500)
        out[3] = 0.
    end

    return PATH_OK, out
end

###############################################################################
###################### SPEED PROFILE FUNCTIONS ################################
###############################################################################

"""
    speed_profile(final_path::Vector{DubinsPathR2}, params::Vector{Number})
Sample single path for plotting

# Arguments
- `final_path::DubinsPathR2`: Vector of dubins paths
- `params::Vector{Float64}`: Vector of speed parameters: [velocity in minimal radius, maximal vehicle velocity, acceleration, deceleration]

# Return
- `times::Vector{Float64}`: times when the speed changes; Inf if invalid path
- `speed_at_time::Vector{Float64}`: speed at each time; Inf if invalid path
"""
function speed_profile(path::DubinsPathR2, params::Vector{Float64})
    speed_profile(path, params[1], params[2], params[3], params[4])
end

function speed_profile(path::DubinsPathR2, v_min::Number, v_max::Number, a1::Float64, a2::Float64)
    # a1 = a_max, a2 = -a_min
    speeds = [speed_by_radius(r, v_min, v_max) for r in path.r]
    lengths = path.lengths

    times::Vector{Float64} = []
    speed_at_time::Vector{Float64} = []

    function push_vals(time, velocity)
        push!(times, time)
        push!(speed_at_time, velocity)
    end

    time = 0
    push_vals(time, speeds[1])

    time += lengths[1] / speeds[1]
    push_vals(time, speeds[1])
    # speedup in between
    sp_up = (speeds[2] - speeds[1])/a1
    sp_down = (speeds[2] - speeds[3])/a2
    sp_up_len = sp_up*(speeds[2] + speeds[1])/2
    sp_down_len = sp_down*(speeds[2] + speeds[3])/2

    if speeds[1] == speeds[3]
    #### CASE WITH SAME RADII

        if (sp_up_len + sp_down_len <= lengths[2])
        # full speedup (to maximum vs)
            time += sp_up
            push_vals(time, speeds[2])
            time += (lengths[2] - sp_down_len - sp_up_len) / speeds[2]
            push_vals(time, speeds[2])
            time += sp_down
            push_vals(time, speeds[3])
            time += lengths[3]/speeds[3]
            push_vals(time, speeds[3])

        else
        # small speedup (only to have enough time to speedup and speeddown)
        # find the intersect where it is necessary to change from speedup to speeddown
            vs = speeds[1] # start speed = end_speed
            s = lengths[2]
            # quadratic equation coeficients a, b, c, from the uniform speedup formula
            xa = a1 + a2 # xb = 0
            xc = -(vs^2)*xa - 2*a1*a2*s
            det = -4*xa*xc
            if det < 0
                return Inf, Inf
            end

            vx = sqrt(det) / (2*xa) # second root is smaller than beggining speed
            time += (vx - vs) / a1
            push_vals(time, vx)
            time += (vx - vs) / a2
            push_vals(time, vs)
            time += lengths[3] / vs
            push_vals(time, vs)
        end

    elseif speeds[1] < speeds[3]
    #### CASE WITH SPEEDUP TO THE LATTER SPEED
        # compute if straight part + latter part is long enough to speedup
        sp_up_finish = (speeds[3] - speeds[1])/a1
        sp_up_finish_len = sp_up_finish*(speeds[3] + speeds[1])/2

        # full speedup (to maximum vs)
        if (sp_up_len + sp_down_len) <= lengths[2]
            time += sp_up
            push_vals(time, speeds[2])
            time += (lengths[2] - sp_down_len - sp_up_len) / speeds[2]
            push_vals(time, speeds[2])
            time += sp_down
            push_vals(time, speeds[3])
            time += lengths[3] / speeds[3]
            push_vals(time, speeds[3])

        elseif sp_up_finish_len > lengths[2] && (sp_up_finish_len < lengths[2] + lengths[3])
            time += sp_up_finish
            push_vals(time, speeds[3])
            time += (lengths[3] + lengths[2] - sp_up_finish_len) / speeds[3]
            push_vals(time, speeds[3])

        elseif sp_up_finish_len < lengths[2] # small speedup (only to have enough time to speedup and speeddown)

            vs = speeds[1]
            ve = speeds[3]
            s = lengths[2]
            xa = a1 + a2 # xb = 0
            xc = -(vs^2)*a2 -(ve^2)*a1 - 2*a1*a2*s
            det = -1 * (4*xa*xc) #
            if det < 0
                @warn "det < 0"
                return Inf, Inf
            end

            vx = sqrt(det) / (2*xa) # second root is smaller than previous speed - unnecessary
            if vx < ve
                return Inf, Inf
            end
            time += abs(vx - vs) / a1
            push_vals(time, vx)
            time += abs(vx - ve) / a2
            push_vals(time, ve)
            time += lengths[3] / ve
            push_vals(time, ve)
        else
            return Inf, Inf
        end

    elseif speeds[1] > speeds[3]
    #### CASE WITH SPEEDDOWN TO THE LATTER SPEED
        sp_down_finish = (speeds[1] - speeds[3])/a2
        sp_down_finish_len = sp_down_finish*(speeds[3] + speeds[1])/2

        if (sp_up_len + sp_down_len) < lengths[2]
            time += lengths[1] / speeds[1]
            push_vals(time, speeds[1])
            time += sp_up
            push_vals(time, speeds[2])
            time += (lengths[2] - sp_down_len - sp_up_len) / speeds[2]
            push_vals(time, speeds[2])
            time += sp_down
            push_vals(time, speeds[3])
            time += lengths[3] / speeds[3]
            push_vals(time, speeds[3])
        elseif sp_down_finish_len > lengths[2] && sp_down_finish_len < (lengths[1] + lengths[2])
            time += (lengths[1] + lengths[2] - sp_down_finish_len) / speeds[1]
            push_vals(time, speeds[1])
            time += sp_down_finish
            push_vals(time, speeds[3])
            time += lengths[3] / speeds[3]
            push_vals(time, speeds[3])
        elseif sp_down_finish_len < lengths[2]
            vs = speeds[1]
            ve = speeds[3]
            s = lengths[2]
            xa = a1 + a2 # xb = 0
            xc = -(vs^2)*a2 -(ve^2)*a1 - 2*a1*a2*s
            det = -1 * (4*xa*xc)
            if det < 0
                return Inf, Inf
            end
            vx = sqrt(det) / (2*xa) # second root is smaller than previous speed - unnecessary
            if vx < ve
                return Inf, Inf
            end

            time += lengths[1] / speeds[1]
            push_vals(time, speeds[1])
            time += abs(vx - vs) / a1
            push_vals(time, vx)
            time += abs(vx - ve) / a2
            push_vals(time, ve)
            time += lengths[3] / speeds[3]
            push_vals(time, speeds[3])
        else
            return Inf, Inf
        end

    end
    return times, speed_at_time
end

end


text:
@testset "BatchQueryStrategies" begin
    BQS = OneClassActiveLearning.QueryStrategies.BatchQueryStrategies

    dummy_data, dummy_labels = load_data(TEST_DATA_FILE)

    classifier = SVDD.SVDDneg(dummy_data, dummy_labels)
    init_strategy = SimpleCombinedStrategy(FixedGammaStrategy(MLKernels.GaussianKernel(1)), FixedCStrategy(0.5))
    initialize!(classifier, init_strategy)

    empty_params = Dict{Symbol, Any}()
    sequential_strategy = Dict{Symbol, Any}(
        :type => :RandomPQs,
        :param => Dict{Symbol, Any}()
    )
    batch_size = 5

    function test_batch_qs(qs_type, params::Dict{Symbol, Any}, classifier, data, labels, candidate_indices)
        @test_throws Union{UndefKeywordError, MethodError, ArgumentError} initialize_qs(qs_type, classifier, data, empty_params)
        qs = initialize_qs(qs_type, classifier, data, Dict(params))
        for candidates in candidate_indices
            batch_indices = select_batch(qs, data, labels, candidates)
            @test length(batch_indices) == min(length(candidates), batch_size)
            @test issubset(batch_indices, Set(candidates))
        end
    end

    @testset "batch_qs" begin
        labels = labelmap(fill(:U, 10))
        candidate_indices = [collect(inds) for inds in [1:3, 1:10]]

        @testset "RandomBatchQs" begin
            qs_type = RandomBatchQs
            params = Dict{Symbol, Any}(
                :k => batch_size
            )
            test_batch_qs(qs_type, params, classifier, dummy_data, labels, candidate_indices)
        end

        @testset "TopK" begin
            qs_type = TopKBatchQs
            params = Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy
            )
            test_batch_qs(qs_type, params, classifier, dummy_data, labels, candidate_indices)
        end

        @testset "GappedTopKBatchQs" begin
            qs_type = GappedTopKBatchQs
            params = Dict{Symbol, Any}(
                :k => batch_size,
                :m => 0,
                :SequentialStrategy => sequential_strategy
            )
            @test_throws ArgumentError initialize_qs(qs_type, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(params..., :m => -1))
            test_batch_qs(qs_type, params, classifier, dummy_data, labels, candidate_indices)
        end

        @testset "ClusterBatchQs" begin
            qs_type = ClusterBatchQs
            params = Dict{Symbol, Any}(
                :k => batch_size
            )
            @test_throws UndefKeywordError initialize_qs(qs_type, SVDD.RandomOCClassifier(dummy_data), dummy_data, empty_params)
            test_batch_qs(qs_type, params, classifier, dummy_data, labels, candidate_indices)
        end

        @testset "ClusterTopKBatchQs" begin
            qs_type = ClusterTopKBatchQs
            params = Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy
            )
            qs = initialize_qs(qs_type, SVDD.RandomOCClassifier(dummy_data), dummy_data, params)
            test_batch_qs(qs_type, params, classifier, dummy_data, labels, candidate_indices)
        end

        @testset "EnsembleBatchQs" begin
            qs_type = EnsembleBatchQs
            params = Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy,
                :solver => TEST_SOLVER
            )
            qs = initialize_qs(EnsembleBatchQs, classifier, dummy_data, params)
            labels_ensemble = labelmap(fill(:U, TEST_DATA_NUM_OBSERVATIONS))
            for candidates in candidate_indices
                try
                    batch_indices = select_batch(qs, dummy_data, labels_ensemble, candidates)
                    @test length(batch_indices) == min(length(candidates), batch_size)
                    @test issubset(batch_indices, Set(candidates))
                catch e
                    # This strategy depends on the solver finding a solution
                    # for the submodels.
                    # The standard Ipopt solver does not always find one.
                    # If no solution is found, this strategy simply does not
                    # work and the use of another solver is recommended.
                    # In that case we skip these tests.
                    if e isa ErrorException
                        println(e)
                        println("Solver could not find a solution, skipping test.")
                        continue
                    else
                        throw(e)
                    end
                end
            end
        end

        @testset "weighted sum objective" begin
            qs_types = [IterativeBatchQs, EnumerativeBatchQs]
            params = [Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy,
                :representativeness => rep,
                :diversity => div,
                :λ_inf => 1,
                :λ_rep => 1,
                :λ_div => 1
            ) for rep in [:KDE] for div in [:AngleDiversity, :EuclideanDistance]]
            invalid_rep_params = Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy,
                :representativeness => :INVALID,
                :diversity => :AngleDiversity,
                :λ_inf => 1,
                :λ_rep => 1,
                :λ_div => 1)
            invalid_div_params = Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy,
                :representativeness => :KDE,
                :diversity => :INVALID,
                :λ_inf => 1,
                :λ_rep => 1,
                :λ_div => 1)
            for qs_type in qs_types
                @testset "$(qs_type)" begin
                    for invalid_params in [invalid_rep_params, invalid_div_params]
                        @test_throws ArgumentError initialize_qs(qs_type, classifier, dummy_data, invalid_params)
                    end
                    for filled_params in params
                        @testset "$(filled_params)" begin
                            test_batch_qs(qs_type, filled_params, classifier, dummy_data, labels, candidate_indices)
                        end
                    end
                end
            end
        end

        @testset "hierarchical" begin
            qs_types = [FilterHierarchicalBatchQs, EnumFilterHierarchicalBatchQs]
            params = Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy,
                :representativeness => :KDE,
                :diversity => :EuclideanDistance
            )
            invalid_rep_params = Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy,
                :representativeness => :INVALID,
                :diversity => :AngleDiversity)
            invalid_div_params = Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy,
                :representativeness => :KDE,
                :diversity => :INVALID)
            for qs_type in qs_types
                @testset "$(qs_type)" begin
                    for invalid_params in [invalid_rep_params, invalid_div_params]
                        @test_throws ArgumentError initialize_qs(qs_type, classifier, dummy_data, invalid_params)
                    end
                    test_batch_qs(qs_type, params, classifier, dummy_data, labels, candidate_indices)
                end
            end
        end

        @testset "FilterSimilarBatchQs" begin
            qs_type = FilterSimilarBatchQs
            params = [Dict{Symbol, Any}(
                :k => batch_size,
                :SequentialStrategy => sequential_strategy,
                :diversity => div
            ) for div in [:AngleDiversity, :EuclideanDistance]]
            @test_throws UndefVarError initialize_qs(qs_type, classifier, dummy_data, invalid_div_params)
            for filled_params in params
                @testset "$(filled_params)" begin
                    test_batch_qs(qs_type, filled_params, classifier, dummy_data, labels, candidate_indices)
                end
            end
        end
    end

    @testset "get_query_objects" begin
        qs = RandomBatchQs(k=4)
        @testset "a" begin
            data = rand(2, 6)
            pools = fill(:U, 6)
            pools[5] = :Lin
            indices = [1, 3, 5, 7, 9, 11]
            history = [[7]]
            @test_throws ArgumentError OneClassActiveLearning.get_query_objects(qs, data, fill(:Lin, 5), indices, history)
            @test OneClassActiveLearning.get_query_objects(qs, data, pools, indices, history) == [1, 3, 5, 11]
        end

        @testset "b" begin
            data = rand(2, 5)
            pools = fill(:U, 5)
            pools[5] = :Lin
            indices = [1, 2, 4, 7, 9]
            history = [[7]]
            @test_throws ArgumentError OneClassActiveLearning.get_query_objects(qs, data, fill(:Lin, 5), indices, history)
            @test OneClassActiveLearning.get_query_objects(qs, data, pools, indices, history) == [1, 2, 4]
        end
    end
end

ext: jl


text:
using PyCall
using JSON
using Pipe
using HTTP

function getSimpleResults(results,geneset="KEGG_2015")
    @pipe results[1][geneset] |> map(r->r[2:3],_)
end

function enrichall(userListId::Int, gene_set_libraries::Array{String,1} = ["KEGG_2015"])

    results = [enrich(userListId, gene_set_library) for gene_set_library in gene_set_libraries]
end

function enrich(userListId::Int, gene_set_library::String = "KEGG_2015")

    ENRICHR_URL = "http://amp.pharm.mssm.edu/Enrichr/enrich"
    query_string = "?userListId=$(string(userListId))&backgroundType=$(gene_set_library)"

    r = HTTP.request("GET", ENRICHR_URL * query_string)

    r.body |> String |> JSON.parse
end

function list(userListId::Int)
    ENRICHR_URL = "http://amp.pharm.mssm.edu/Enrichr/view?userListId=$(string(userListId))"

    r = HTTP.request("GET", ENRICHR_URL)

    r.body |> String |> JSON.parse
end

function addGeneList(gene_list::Array{String,1}; description)

    py"""
    import requests
    import http

    ENRICHR_URL = "http://amp.pharm.mssm.edu/Enrichr/addList"

    def query(genes_str,description):

        payload = {
            'list': (None, genes_str),
            'description': (None, description)
        }

        response = requests.post(ENRICHR_URL, files=payload)

        if not response.ok:
            raise Exception('Error analyzing gene list')

        return response.text
    """

    genes_str = join(gene_list,"\n")

    response = py"query"(genes_str, description)

    JSON.parse(response)["userListId"]
end

max_stars_repo_name: Fredrikp-ume/Enrich

max_issues_repo_name: Fredrikp-ume/Enrich

max_forks_repo_name: Fredrikp-ume/Enrich


text:
mutable struct ConfigData <: AbstractMetaData
    filepath::AbstractString
    data::Union{Missing, AbstractDict}
    tables::AbstractDict
    mtime::Float64
end
function ConfigData(file)
    x = ConfigData(file, missing, Dict{Symbol, Any}(), 0.)
    if !isfile(file)
        throw(SystemError(file, 2))
    end
    loaddata!(x)
    return x
end

function ismodified(x::ConfigData)
    x.mtime != mtime(filepath(x))
end

function loaddata!(x::ConfigData)
    x.data = JSON.parsefile(filepath(x); dicttype=OrderedDict{String,Any}, use_mmap=false)

    # create environment paths
    for (k, path) in x.data["environment"]
        if isabspath(path)
            fullpath = normpath(path)
        else
            fullpath = abspath(joinpath(GAMEENV["PROJECT"], path))
        end
        if !isdir(fullpath)
            @info "$fullpath is created"
            mkpath(fullpath)
        end
        GAMEENV[uppercase(k)] = fullpath
    end
    # check if xlsx file exists
    for (fname, sheetdata) in x.data["xlsxtables"]
        file = joinpath(GAMEENV["XLSX"], fname)
        if !isfile(file)
            @warn "$file does not exist, please insert a valid path in `config.json`"
        else
            jwb = XLSXTable(file, sheetdata)
            x.tables[_filename(jwb)] = jwb
        end
    end

    return x
end

function xlsxfilenames(x::ConfigData)
    keys(x.tables)
end


text:
export MLP, weightCount, release!, predict!, setWeights!

abstract MLP

weightCount(mlp::MLP) = mlp.data.weightCount

release!(mlp::MLP) = release!(mlp.data)

predict!(mlp::MLP, input) = predict!(mlp, mlp.data, input)

setWeights!(mlp::MLP, weights) = setWeights!(mlp.data, weights)

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
__precompile__()
module Quadmath

export Float128, Complex256

import Base: (*), +, -, /,  <, <=, ==, ^, convert,
          reinterpret, sign_mask, exponent_mask, exponent_one, exponent_half,
          significand_mask,
          promote_rule, widen,
          string, print, show, showcompact, parse,
          acos, acosh, asin, asinh, atan, atanh, cosh, cos,
          exp, expm1, log, log2, log10, log1p, sin, sinh, sqrt,
          tan, tanh,
          ceil, floor, trunc, round, fma,
          atan2, copysign, max, min, hypot,
          abs, imag, real, conj, angle, cis,
          eps, realmin, realmax, isinf, isnan, isfinite

import SpecialFunctions: erf, erfc,
          besselj, besselj0, besselj1, bessely, bessely0, bessely1,
          gamma, lgamma

if Sys.isapple()
    const quadoplib = "libquadmath.0"
    const libquadmath = "libquadmath.0"
elseif Sys.isunix()
    const quadoplib = "libgcc_s.so.1"
    const libquadmath = "libquadmath.so.0"
    const mpfr_float128 = joinpath(dirname(@__FILE__),
                "..", "deps", "lib", "mpfr_float128.so")
elseif Sys.iswindows()
    const quadoplib = "libgcc_s_seh-1.dll"
    const libquadmath = "libquadmath-0.dll"
end

@static if Sys.isunix()
    # we use this slightly cumbersome definition to ensure that the value is passed
    # on the xmm registers, matching the x86_64 ABI for __float128.
    const Cfloat128 = NTuple{2,VecElement{Float64}}

    struct Float128 <: AbstractFloat
        data::Cfloat128
    end

    const Complex256 = Complex{Float128}

    Base.cconvert(::Type{Cfloat128}, x::Float128) = x.data

    # reinterpret
    function reinterpret(::Type{UInt128}, x::Float128)
        hi = reinterpret(UInt64, x.data[2].value)
        lo = reinterpret(UInt64, x.data[1].value)
        UInt128(hi) << 64 | lo
    end
    function reinterpret(::Type{Float128}, x::UInt128)
        fhi = reinterpret(Float64, (x >> 64) % UInt64)
        flo = reinterpret(Float64, x % UInt64)
        Float128((VecElement(flo), VecElement(fhi)))
    end
    reinterpret(::Type{Unsigned}, x::Float128) = reinterpret(UInt128, x)

    reinterpret(::Type{Int128}, x::Float128) =
        reinterpret(Int128, reinterpret(UInt128, x))
    reinterpret(::Type{Float128}, x::Int128) =
        reinterpret(Float128, reinterpret(UInt128, x))

elseif Sys.iswindows()
    primitive type Float128 128 end
    const Cfloat128 = Float128
end

sign_mask(::Type{Float128}) =        0x8000_0000_0000_0000_0000_0000_0000_0000
exponent_mask(::Type{Float128}) =    0x7fff_0000_0000_0000_0000_0000_0000_0000
exponent_one(::Type{Float128}) =     0x3fff_0000_0000_0000_0000_0000_0000_0000
exponent_half(::Type{Float128}) =    0x3ffe_0000_0000_0000_0000_0000_0000_0000
significand_mask(::Type{Float128}) = 0x0000_ffff_ffff_ffff_ffff_ffff_ffff_ffff

fpinttype(::Type{Float128}) = UInt128

# constructors

Float128(x::Float128) = x

Float128(x::Float64) =
    Float128(ccall((:__extenddftf2, quadoplib), Cfloat128, (Cdouble,), x))
Float64(x::Float128) =
    ccall((:__trunctfdf2, quadoplib), Cdouble, (Cfloat128,), x)

Base.Int32(x::Float128) =
    ccall((:__fixtfsi, quadoplib), Int32, (Cfloat128,), x)
Float128(x::Int32) =
    Float128(ccall((:__floatsitf, quadoplib), Cfloat128, (Int32,), x))

Float128(x::UInt32) =
    Float128(ccall((:__floatunsitf, quadoplib), Cfloat128, (UInt32,), x))

Base.Int64(x::Float128) =
    ccall((:__fixtfdi, quadoplib), Int64, (Cfloat128,), x)
Float128(x::Int64) =
    Float128(ccall((:__floatditf, quadoplib), Cfloat128, (Int64,), x))

Float128(x::Rational{T}) where T = Float128(numerator(x))/Float128(denominator(x))

# conversion

convert(::Type{Float128}, x::Float64) = Float128(x)
convert(::Type{Float64}, x::Float128) = Float64(x)
convert(::Type{Int32}, x::Float128) = Int32(x)
convert(::Type{Float128}, x::Int32) = Float128(x)
convert(::Type{Float128}, x::UInt32) = Float128(x)
convert(::Type{Int64}, x::Float128) = Int32(x)
convert(::Type{Float128}, x::Int64) = Float128(x)
convert(::Type{Float128}, x::Type{Rational{T}}) where T = Float128(x)

const ROUNDING_MODE = Cint[0] # TODO: CHECK!!!!

Float128(x::BigFloat) =
    Float128(ccall((:mpfr_get_float128, mpfr_float128), Cfloat128, (Ref{BigFloat},Int32), x, ROUNDING_MODE[end]))

function BigFloat(x::Float128)
    z = BigFloat()
    res = ccall((:mpfr_set_float128, mpfr_float128), Int32, (Ref{BigFloat}, Cfloat128, Int32), z, x, ROUNDING_MODE[end])
    return z
end

convert(::Type{BigFloat}, x::Float128) = BigFloat(x)
convert(::Type{Float128}, x::BigFloat) = Float128(x)

# comparison

(==)(x::Float128, y::Float128) =
    ccall((:__eqtf2,quadoplib), Cint, (Cfloat128,Cfloat128), x, y) == 0

(<)(x::Float128, y::Float128) =
    ccall((:__letf2,quadoplib), Cint, (Cfloat128,Cfloat128), x, y) == -1

(<=)(x::Float128, y::Float128) =
    ccall((:__letf2,quadoplib), Cint, (Cfloat128,Cfloat128), x, y) <= 0

# arithmetic

(+)(x::Float128, y::Float128) =
    Float128(ccall((:__addtf3,quadoplib), Cfloat128, (Cfloat128,Cfloat128), x, y))
(-)(x::Float128, y::Float128) =
    Float128(ccall((:__subtf3,quadoplib), Cfloat128, (Cfloat128,Cfloat128), x, y))
(*)(x::Float128, y::Float128) =
    Float128(ccall((:__multf3,quadoplib), Cfloat128, (Cfloat128,Cfloat128), x, y))
(/)(x::Float128, y::Float128) =
    Float128(ccall((:__divtf3,quadoplib), Cfloat128, (Cfloat128,Cfloat128), x, y))
(-)(x::Float128) =
    Float128(ccall((:__negtf2,quadoplib), Cfloat128, (Cfloat128,), x))
(^)(x::Float128, y::Float128) =
    Float128(ccall((:powq,libquadmath), Cfloat128, (Cfloat128,Cfloat128), x, y))

# math

## one argument
for f in (:acos, :acosh, :asin, :asinh, :atan, :atanh, :cosh, :cos,
          :erf, :erfc, :exp, :expm1, :log, :log2, :log10, :log1p,
          :sin, :sinh, :sqrt, :tan, :tanh,
          :ceil, :floor, :trunc, :lgamma, )
    @eval function $f(x::Float128)
        Float128(ccall(($(string(f,:q)), libquadmath), Cfloat128, (Cfloat128, ), x))
    end
end
for (f,fc) in (:abs => :fabs,
               :round => :rint,
               :gamma => :tgamma,
               :besselj0 => :j0,
               :besselj1 => :j1,
               :bessely0 => :y0,
               :bessely1 => :y1,)
    @eval function $f(x::Float128)
        Float128(ccall(($(string(fc,:q)), libquadmath), Cfloat128, (Cfloat128, ), x))
    end
end

## two argument
for f in (:atan2, :copysign, :hypot, )
    @eval function $f(x::Float128, y::Float128)
       Float128(ccall(($(string(f,:q)), libquadmath), Cfloat128, (Cfloat128, Cfloat128), x, y))
    end
end

## misc
fma(x::Float128, y::Float128, z::Float128) =
    Float128(ccall((:fmaq,libquadmath), Cfloat128, (Cfloat128, Cfloat128, Cfloat128), x, y, z))

isnan(x::Float128) =
    0 != ccall((:isnanq,libquadmath), Cint, (Cfloat128, ), x)
isinf(x::Float128) =
    0 != ccall((:isinfq,libquadmath), Cint, (Cfloat128, ), x)
isfinite(x::Float128) =
    0 != ccall((:finiteq,libquadmath), Cint, (Cfloat128, ), x)

besselj(n::Cint, x::Float128) =
       Float128(ccall((:jnq, libquadmath), Cfloat128, (Cint, Cfloat128), n, x))
bessely(n::Cint, x::Float128) =
       Float128(ccall((:ynq, libquadmath), Cfloat128, (Cint, Cfloat128), n, x))

eps(::Type{Float128}) = reinterpret(Float128, 0x3f8f0000000000000000000000000000)
realmin(::Type{Float128}) = reinterpret(Float128, 0x00010000000000000000000000000000)
realmax(::Type{Float128}) = reinterpret(Float128, 0x7ffeffffffffffffffffffffffffffff)
convert(::Type{Float128}, ::Irrational{:π}) =  reinterpret(Float128, 0x4000921fb54442d18469898cc51701b8)
convert(::Type{Float128}, ::Irrational{:e}) =  reinterpret(Float128, 0x40005bf0a8b1457695355fb8ac404e7a)

ldexp(x::Float128, n::Cint) =
    Float128(ccall((:ldexpq, libquadmath), Cfloat128, (Cfloat128, Cint), x, n))
ldexp(x::Float128, n::Integer) =
    ldexp(x, clamp(n, typemin(Cint), typemax(Cint)) % Cint)

function frexp(x::Float128)
    r = Ref{Cint}()
    Float128(ccall((:frexpq, libquadmath), Cfloat128, (Cfloat128, Ptr{Cint}), x, r))
    return x, Int(r[])
end

promote_rule(::Type{Float128}, ::Type{Float16}) = Float128
promote_rule(::Type{Float128}, ::Type{Float32}) = Float128
promote_rule(::Type{Float128}, ::Type{Float64}) = Float128
promote_rule(::Type{Float128}, ::Type{<:Integer}) = Float128

#widen(::Type{Float64}) = Float128
widen(::Type{Float128}) = BigFloat

# TODO: need to do this better
function parse(::Type{Float128}, s::AbstractString)
    Float128(ccall((:strtoflt128, libquadmath), Cfloat128, (Cstring, Ptr{Ptr{Cchar}}), s, C_NULL))
end

function string(x::Float128)
    lng = 64
    buf = Array{UInt8}(undef, lng + 1)
    lng = ccall((:quadmath_snprintf,libquadmath), Cint, (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Cfloat128...), buf, lng + 1, "%.35Qe", x)
    return unsafe_string(pointer(buf), lng)
end

print(io::IO, b::Float128) = print(io, string(b))
show(io::IO, b::Float128) = print(io, string(b))
showcompact(io::IO, b::Float128) = print(io, string(b))

end # modeule Quadmath


text:
using DrWatson
@quickactivate "EPExperiments"

using EPExperiments
using Random
using Parameters
using Distributions
using AnnealedIS
using StatsPlots: savefig
using StatsFuns: logsumexp, logistic
using StatsBase
using Logging
using LoggingExtras
using AdvancedMH
using AdvancedHMC
using Dates
using JLD2
using FileIO

const LOG_FILENAME = "out.log"

function estimate_sir(
    experiment_config::SIRExperimentConfig
)
    @unpack num_samples, num_betas, anis_params, seed, cost_fn_name = experiment_config
    Random.seed!(seed)

    cost_fn = COST_FUNCTIONS[cost_fn_name]

    input_data = load(datadir("sims", experiment_config.data_filename))

    @unpack tmax, total_population, i_0_true, β_true, γ = input_data["params"]
    X, Y = input_data["X"], input_data["Y"]

    result_folder = mkpath(joinpath(
        datadir("inferences"),
        savename(experiment_config)
    ))
    logger = TeeLogger(
        ConsoleLogger(),
        FileLogger(joinpath(result_folder, LOG_FILENAME))
    )

    # Prior predictive simulation
    ode_prior = sample(bayes_sir(Y, total_population, γ, cost_fn).gamma2, Prior(), 1000)
    obstimes = 1.0:1.0:tmax
    savefig(
        plot_predictive(obstimes, X, Y, ode_prior),
        joinpath(result_folder, "prior_pred.png")
    )

    # We first create a geometric spacing between 1 and 1001 because directly
    # doing it between 0 and 1 gives numerical problems.
    betas = (geomspace(1, 1001, num_betas) .- 1) ./ 1000

    anis_alg = get_anis_alg(anis_params, betas)
    tabi = TABI(
        TuringAlgorithm(anis_alg, num_samples),
        TuringAlgorithm(anis_alg, 0),
        TuringAlgorithm(anis_alg, num_samples)
    )

    @time begin
    R0_estimate, ode_tabi = estimate_expectation(
        bayes_sir(Y, total_population, γ, cost_fn), tabi; store_intermediate_samples=true
    )
    end

    # Save experimental data
    @tagsave(joinpath(result_folder, "results.jld2"), Dict(
        "ode_tabi" => ode_tabi,
        "R0_estimate" => R0_estimate,
        "experiment_config" => experiment_config
    ))

    process_results(
        result_folder,
        logger,
        ode_tabi=ode_tabi,
        R0_estimate=R0_estimate,
        tabi=tabi,
        X=X,
        Y=Y,
        obstimes=obstimes,
        i_0_true=i_0_true,
        β_true=β_true,
        cost_fn=((i, b) -> cost_fn(b, γ))
    )
end

estimate_sir(SIRExperimentConfig(
    experiment_name="trial_run",
    num_samples=1000,
    anis_params=AnISHMCConfig{Turing.ForwardDiffAD{2}}(
        proposal=AdvancedHMC.StaticTrajectory(AdvancedHMC.Leapfrog(0.05), 10)
    ),
    cost_fn_name=:cost_fn_sir,
    seed=1298
))


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SourcePath(domain, sources)

A path over a spatial `domain` that starts at given source
locations `sources` and progresses outwards.
"""
struct SourcePath{D<:AbstractDomain} <: AbstractPath{D}
  domain::D
  sources::Vector{Int}

  # state fields
  path::Vector{Int}

  function SourcePath{D}(domain, sources) where {D<:AbstractDomain}
    @assert all(1 .≤ sources .≤ npoints(domain)) "sources must be valid locations"
    @assert length(unique(sources)) == length(sources) "non-unique sources"
    @assert length(sources) ≤ npoints(domain) "more sources than points in domain"

    # coordinate matrix for source points
    S = coordinates(domain, sources)

    # fit search tree
    kdtree = KDTree(S)

    # other locations that are not sources
    others = setdiff(1:npoints(domain), sources)

    # process points in chunks
    chunksize = 10^3
    chunks = Iterators.partition(others, chunksize)

    # pre-allocate memory for coordinates
    X = Matrix{coordtype(domain)}(undef, ndims(domain), chunksize)

    # compute distances to sources
    dists = []
    for chunk in chunks
      coordinates!(X, domain, chunk)
      _, ds = knn(kdtree, view(X,:,1:length(chunk)), length(sources), true)
      append!(dists, ds)
    end

    path = vcat(sources, view(others, sortperm(dists)))

    new(domain, sources, path)
  end
end

SourcePath(domain, sources) = SourcePath{typeof(domain)}(domain, sources)

Base.iterate(p::SourcePath, state=1) = state > npoints(p.domain) ? nothing : (p.path[state], state + 1)

# ------------
# IO methods
# ------------
function Base.show(io::IO, path::SourcePath)
  print(io, "SourcePath")
end

function Base.show(io::IO, ::MIME"text/plain", path::SourcePath)
  println(io, path)
end


text:
#const RETAINED_COLS_BB_SEO = []
const BB_DATE_FORMAT = dateformat"mm/dd/yyyy"
const BB_EVENT_NAMES = ["seo-nasdaq", "seo-otc", "seo-ny"]
const BB_EVENT_RETAINED = [:date, :adate, :offertype, :filedate, :permno, :gvkey, :offersizem,
  :primarysharesm, :securitytype, :primaryexchange, :issuerticker]
const SEO_DAYS = 7
const SEO_RET = :ret
const SEO_MIN_ALLOWED = SEO_DAYS*2

function cleanbbnames!(bb::DataFrame; )
  bbnames::Vector{String} = (string).(names(bb))

  bbnames .= (lowercase).(bbnames)
  bbnames .= (s->replace(s, r"[^a-z0-9]"=>"")).(bbnames)
  names!(bb, (Symbol).(bbnames))

end

function filterbbeventdates!(bb::DataFrame, bbdateformat = BB_DATE_FORMAT,
    yearrange::UnitRange{Int} = YEAR_RANGE[])
  #initial filters to ensure complete data
  filter!(r::DataFrameRow->
    (!ismissing(r.offersizem)) && (r.offersizem > 0.) &&
    (!ismissing(r.offertype)) && ("IPO" ≠ r.offertype[1:3]) &&
    (!occursin("Secondary", r.offertype)) &&
    (!ismissing(r.primaryexchange)) && ("OTC" ≠ r.primaryexchange[1:3]) &&
    (length(r.cusip) == 9) &&
    (!ismissing(r.issuerticker)) &&
    (!ismissing(r.primarysharesm)) && (r.primarysharesm > 0.) &&
    (r.offerstage == "Trading" )
    , bb)

  #begin processing the dates
  for f ∈ [:launchdate, :filingdate, :filingtermdate,
      :pricingdate, :tradedateus, :effectivedate, :announceddate]

    bb[!, f] = (s::MString->ismissing(s) ? missing : Date(s, bbdateformat)).(bb[!,f])
  end

  rename!(bb, :announceddate=>:adate)
  bb[!, :filedate] = (r::DataFrameRow->coalesce(r.filingtermdate, r.filingdate, r.launchdate)).(eachrow(bb))
  bb[!, :date] = (r::DataFrameRow->coalesce(r.pricingdate, r.tradedateus, r.effectivedate)).(eachrow(bb))

  #issue and file date must be at least 10 days appart
  filter!(r::DataFrameRow->
    (!ismissing(r.filedate)) &&  (!ismissing(r.date)) && ((r.date - r.filedate).value ≥ 10) &&
    (!ismissing(r.adate)) && (r.date > r.adate)
    , bb)

  filter!(r::DataFrameRow->year(r.adate) ∈ yearrange, bb) #restrict to desired years
  #WARNING: STUB -Need to acquire data and put in the filter for 10 days between filing and offer date
  return bb
end

function getidentifiersbb(bb::DataFrame)
  local ccm::DataFrame = loadccm()
  ccm = preprocessccm!(ccm)

  bb.eid = collect(1:size(bb,1)) #assign an id to each event
  ccm.cusip = (Symbol).(ccm.cusip)
  bb.cusip = (Symbol).(bb.cusip)

  bb = join(bb, ccm, on=:cusip=>:cusip)
  filter!(r::DataFrameRow -> r.date ∈ r.linkeffdate:Day(1):r.linkenddate, bb)
  @assert allunique(bb.eid) #no dups

  rename!(bb, :lpermno=>:permno) #sanitize permno

  return bb
end

function readbbseo(fnames::Vector{String} = BB_EVENT_NAMES;
  path::String = EVENT_PATH,
  bbretained::Vector{Symbol} = BB_EVENT_RETAINED)

  local bb::DataFrame = DataFrame()

  #read in the bloomberg files
  for (i, f) ∈ enumerate(fnames)
    bbin::DataFrame = CSV.read("$path\\$f.csv") |> DataFrame
    cleanbbnames!(bbin)
    bb = vcat(bb, bbin)
  end

  #println("got here")

  filterbbeventdates!(bb)

  bb = getidentifiersbb(bb)
  select!(bb, bbretained)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
Module for performing a Doppler-constrained PCA analysis

For algorithm information, see:
  - Jones, Stenning, Ford et al. https://arxiv.org/abs/1711.01318
  - See also: Gilbertson, Ford, Jones & Stenning 2020 https://arxiv.org/abs/2009.01085

Author: Eric Ford
Date:   September 2020
"""
module DCPCA

using Statistics, MultivariateStats

#import .RvSpectML : speed_of_light_mps
const speed_of_light_mps = 299792458.0 # TODO: Update value

export doppler_constrained_pca, compute_spectra_perp_doppler_shift

function compute_spectra_perp_doppler_shift(spectra::AA, deriv::V1, rvs::V2) where {
            T1<:Real, AA<:AbstractArray{T1,2}, T2<:Real, V1<:AbstractVector{T2}, T3<:Real, V2<:AbstractVector{T3} }
   @assert size(spectra,1) == length(deriv)
   @assert size(spectra,2) == length(rvs)
   fm_perp = spectra .- rvs' .* deriv./speed_of_light_mps
end

function doppler_constrained_pca(flux::AbstractArray{T1,2}, deriv::AbstractVector{T2}, rvs::AbstractVector{T3} )  where { T1<:Real, T2<:Real, T3<:Real  }#, times::AbstractVector, rvs::AbstractVector)
  rv_shift = rvs .- mean(rvs)
  fm_perp = compute_spectra_perp_doppler_shift(flux,deriv, rv_shift )
  idx_good_obs = 1:length(rvs)
  M = fit(PCA, fm_perp[:,idx_good_obs]; maxoutdim=10)
  pca_out = MultivariateStats.transform(M,fm_perp[:,idx_good_obs])
  return pca_out, M
end

end # DCPCA


text:
const gases = Vector{Symbol}()

const DefaultPerfectGas = PerfectGas()
const Air = PerfectGas()
export DefaultPerfectGas, Air
push!(gases,:Air)

# To add a gas, give it a name, the specific heat ratio, and the molar mass
@gas He  20u"°C"   1.94e-5   5//3    4.002602u"g/mol"
@gas O2  20u"°C"   2.04e-5   1.395   31.999u"g/mol"
@gas CO2 20u"°C"   1.47e-5   1.289   44.01u"g/mol"
@gas H2  20u"°C"   8.84e-6   1.405   2.016u"g/mol"
@gas Methane 20u"°C" 1.10e-5 1.299   16.043u"g/mol"
@gas N2  20u"°C"   1.76e-5   1.40    28.0134u"g/mol"
@gas Ar   300u"K"  2.27e-5   5//3    39.948u"g/mol"

max_stars_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_issues_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_forks_repo_name: UCLAMAEThreads/ThermofluidQuantities


text:
export MultirateRungeKutta

LSRK2N = LowStorageRungeKutta2N

"""
    MultirateRungeKutta(slow_solver, fast_solver; dt, t0 = 0)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q} = f_fast(Q, t) + f_slow(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a multirate Runge-Kutta scheme using two different RK
solvers. This is based on

Currently only the low storage RK methods can be used as slow solvers

### References
 - [Schlegel2012](@cite)
"""
mutable struct MultirateRungeKutta{SS, FS, RT} <: AbstractODESolver
    "slow solver"
    slow_solver::SS
    "fast solver"
    fast_solver::FS
    "time step"
    dt::RT
    "time"
    t::RT
    "elapsed time steps"
    steps::Int

    function MultirateRungeKutta(
        slow_solver::LSRK2N,
        fast_solver,
        Q = nothing;
        dt = getdt(slow_solver),
        t0 = slow_solver.t,
    ) where {AT <: AbstractArray}
        SS = typeof(slow_solver)
        FS = typeof(fast_solver)
        RT = real(eltype(slow_solver.dQ))
        new{SS, FS, RT}(slow_solver, fast_solver, RT(dt), RT(t0), 0)
    end
end

function MultirateRungeKutta(
    solvers::Tuple,
    Q = nothing;
    dt = getdt(solvers[1]),
    t0 = solvers[1].t,
) where {AT <: AbstractArray}
    if length(solvers) < 2
        error("Must specify atleast two solvers")
    elseif length(solvers) == 2
        fast_solver = solvers[2]
    else
        fast_solver = MultirateRungeKutta(solvers[2:end], Q; dt = dt, t0 = t0)
    end

    slow_solver = solvers[1]

    MultirateRungeKutta(slow_solver, fast_solver, Q; dt = dt, t0 = t0)
end

function MultirateRungeKutta(
    mrk,
    op::TimeScaledRHS{2, RT} where {RT},
    Q = nothing;
    dt = 0,
    t0 = 0,
    steps = 0,
) where {AT <: AbstractArray}

    slow_solver = mrk(op.rhs![1], Q, dt = dt, t0 = t0)
    fast_solver = mrk(op.rhs![2], Q, dt = dt, t0 = t0)

    MultirateRungeKutta(slow_solver, fast_solver, Q; dt = dt, t0 = t0)
end

function dostep!(
    Q,
    mrrk::MultirateRungeKutta{SS},
    param,
    time::Real,
    nsteps::Int,
    iStage::Int,
    slow_δ = nothing,
    slow_rv_dQ = nothing,
    slow_scaling = nothing,
) where {SS <: LSRK2N}
    for i in 1:nsteps
        dostep!(Q, mrrk, param, time, slow_δ, slow_rv_dQ, slow_scaling)
        time += mrrk.fast_solver.dt
    end
end

function dostep!(
    Q,
    mrrk::MultirateRungeKutta{SS},
    param,
    time,
    in_slow_δ = nothing,
    in_slow_rv_dQ = nothing,
    in_slow_scaling = nothing,
) where {SS <: LSRK2N}
    dt = mrrk.dt

    slow = mrrk.slow_solver
    fast = mrrk.fast_solver

    slow_rv_dQ = realview(slow.dQ)

    groupsize = 256

    fast_dt_in = getdt(fast)

    for slow_s in 1:length(slow.RKA)
        # Currnent slow state time
        slow_stage_time = time + slow.RKC[slow_s] * dt

        # Evaluate the slow mode
        slow.rhs!(slow.dQ, Q, param, slow_stage_time, increment = true)

        if in_slow_δ !== nothing
            slow_scaling = nothing
            if slow_s == length(slow.RKA)
                slow_scaling = in_slow_scaling
            end
            # update solution and scale RHS
            event = Event(array_device(Q))
            event = update!(array_device(Q), groupsize)(
                slow_rv_dQ,
                in_slow_rv_dQ,
                in_slow_δ,
                slow_scaling;
                ndrange = length(realview(Q)),
                dependencies = (event,),
            )
            wait(array_device(Q), event)
        end

        # Fractional time for slow stage
        if slow_s == length(slow.RKA)
            γ = 1 - slow.RKC[slow_s]
        else
            γ = slow.RKC[slow_s + 1] - slow.RKC[slow_s]
        end

        # RKB for the slow with fractional time factor remove (since full
        # integration of fast will result in scaling by γ)
        slow_δ = slow.RKB[slow_s] / (γ)

        # RKB for the slow with fractional time factor remove (since full
        # integration of fast will result in scaling by γ)
        nsubsteps = fast_dt_in > 0 ? ceil(Int, γ * dt / fast_dt_in) : 1
        fast_dt = γ * dt / nsubsteps

        updatedt!(fast, fast_dt)

        for substep in 1:nsubsteps
            slow_rka = nothing
            if substep == nsubsteps
                slow_rka = slow.RKA[slow_s % length(slow.RKA) + 1]
            end
            fast_time = slow_stage_time + (substep - 1) * fast_dt
            dostep!(Q, fast, param, fast_time, slow_δ, slow_rv_dQ, slow_rka)
        end
    end
    updatedt!(fast, fast_dt_in)
end

@kernel function update!(fast_dQ, slow_dQ, δ, slow_rka = nothing)
    i = @index(Global, Linear)
    @inbounds begin
        fast_dQ[i] += δ * slow_dQ[i]
        if slow_rka !== nothing
            slow_dQ[i] *= slow_rka
        end
    end
end


text:
# Various simple graph functions associated with graph coloring.

export bipartition, two_color, greedy_color, random_greedy_color

# Create a two-coloring of a graph or die trying. Returns a map from
# the vertex set to the set {1,2}, or error if no such mapping exists.

"""
`two_color(G)` creates a two-coloring of the graph or throws an error
if tghe graph is not bipartite. The output is a `Dict` mapping the
vertex set to the values 1 and 2.
"""
function two_color{T}(G::SimpleGraph{T})
    f = Dict{T,Int}()
    for A in components(G)
        a = first(A)
        f[a] = 1
        Q = Deque{T}()
        push!(Q,a)
        while length(Q)>0
            v = pop!(Q)
            Nv = G[v]
            for w in Nv
                if haskey(f,w)
                    if f[w]==f[v]
                        error("Graph is not bipartite")
                    end
                else
                    f[w] = 3-f[v]
                    push!(Q,w)
                end
            end
        end
    end
    return f
end

# Create a bipartition of a graph or die trying. Returns a set {X,Y}
# that is a bipartition of the vertex set of G.

using SimplePartitions
"""
`bipartition(G)` creates a bipartition of the graph (or returns an
error if the graph is not bipartite. Output is a `Partition`.
"""
function bipartition{T}(G::SimpleGraph{T})
  f = two_color(G)
  return Partition(f)
end

# Color a graph by the greedy algorithm in the sequence specified by
# seq. The array seq must be a permutation of G.V. We don't check
# that's true!

"""
`greedy_color(G,seq)` creates a greedy proper coloring of the
graph. The argument `seq` should be a 1-dimensional array containing
every vertex exactly once (the function does not check this for
you). The function follows that order in creating the coloring which
is returned to you as a `Dict` mapping vertices to positive integers
(representing the colors).

If `seq` is omitted, a random permutation of the vertices is used.
"""
function greedy_color{T}(G::SimpleGraph{T}, seq::Array{T,1})
    f = Dict{T,Int}()  # this is the mapping from V to colors
    maxf::Int = 0      # largest color used

    for v in seq
        colors_used = falses(maxf)  # array if color is used by N[v]
        for w in G[v]
            if haskey(f,w)  # w already colored
                colors_used[f[w]]=true  # mark that color is used
            end
        end
        # give first unused color to v
        for k in 1:maxf
            if colors_used[k] == false
                f[v] = k
                break
            end
        end
        # but if that fails, extend the number of colors available
        if !haskey(f,v)
            maxf += 1
            f[v] = maxf
        end
    end
    return f
end

# This function returns a list of the vertices of G in descending
# order by degree. The order of vertices of the same degree is
# indeterminate. NOTE: This is not exported from this module. Should
# it be?
function deg_sorted_vlist(G::SimpleGraph)
    bye = x -> -x[1]
    list = [ (deg(G,v) , v) for v in G.V ]
    sort!(list, by=bye)
    outlist = [ item[2] for item in list ]
    return outlist
end

# Apply greedy_color to the graph visiting the vertices in decreasing
# order of degree.
function greedy_color{T}(G::SimpleGraph{T})
    seq = deg_sorted_vlist(G)
    return greedy_color(G,seq)
end

# Generate multiple random orders of the vertex set and apply
# greedy_color; return one that uses the fewest colors. This do as
# well as or better than greedy_color on some decreasing order of
# degree.

"""
`random_greedy_color(G, reps)` repeatedly invokes `greedy_color` using
random permutations of the vertex set. After `reps` iterations, the
best coloring found is returned.
"""
function random_greedy_color{T}(G::SimpleGraph{T}, reps::Int=1)
    n = NV(G)
    bestf = greedy_color(G)  # degree order default start
    best  = maximum(values(bestf))
    seq = vlist(G)
    println("Initial coloring uses ", best, " colors")

    for k in 1:reps
        shuffle!(seq)
        f = greedy_color(G,seq)
        mx = maximum(values(f))
        if mx < best
            bestf = f
            best = mx
            println("Reduced to ", best, " colors")
        end
    end
    return bestf
end


text:
using GEE
using Test

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


text:
## Julia program for Calculation of Gamma function
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-21
##
## refer to "Coefficients for the Lanczos Approximation to the Gamma Function"
##   https://mrob.com/pub/ries/lanczos-gamma.html

const LG_g = 5.0;
const LG_N = 6;
const lct = [
    1.000000000190015,
    76.18009172947146,
    -86.50532032941677,
    24.01409824083091,
    -1.231739572450155,
    0.1208650973866179e-2,
    -0.5395239384953e-5];
const ln_sqrt_2_pi = 0.91893853320467274178;
const g_pi = 3.14159265358979323846;

function lanczos_ln_gamma(z)
    if z < 0.5
        #Use Euler's reflection formula:
        #Gamma(z) = Pi / [Sin[Pi*z] * Gamma[1-z]];
        return log(g_pi / sin(g_pi * z)) - lanczos_ln_gamma(1.0 - z);
    end
    z = z - 1.0;
    base = z + LG_g + 0.5;  # Base of the Lanczos exponential
    sum = 0;
    ## We start with the terms that have the smallest coefficients and largest denominator.
    for i = 2:LG_N+1
      sum += lct[i] / (z + i - 1);
    end
    sum += lct[1];
    ## This printf is just for debugging
#    @printf("ls2p %7g  l(b^e) %7g   -b %7g  l(s) %7g\n", ln_sqrt_2_pi,
#              log(base)*(z+0.5), -base, log(sum));
    ## Gamma[z] = Sqrt(2*Pi) * sum * base^[z + 0.5] / E^base
    return ((ln_sqrt_2_pi + log(sum)) - base) + log(base)*(z+0.5);
end

function lanczos_gamma(z)
    return(exp(lanczos_ln_gamma(z)));
end

## example
lanczos_gamma(1)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
export solve

"solve(d::Device, pol::Polarization)"
function solve(d::Device, pol::Polarization=TM)
    (ϵ₀, μ₀, c₀) = normalize_parameters(d);

    Nω = length(d.ω);
    pol == TM && ( fields = Array{FieldTM}(undef, Nω) );
    pol == TE && ( fields = Array{FieldTE}(undef, Nω) );

    for i in eachindex(d.ω)
        @info "Frequency: $i/$Nω"
        ω = d.ω[i];

        length(d.modes) > 0 && ( d.src=zeros(Complex, size(d.grid)) ) # Only reset if we are using modes
        for mode in d.modes
            # TODO: handle the polarization here
            setup_mode!(d, TM, ω, mode.neff, mode.pt, mode.dir, mode.width);
        end

        Tϵ = sparse(Diagonal(ϵ₀*d.ϵᵣ[:]));
        Tϵxi = sparse(Diagonal(1 ./ grid_average(ϵ₀*d.ϵᵣ, x̂)[:]));
        Tϵyi = sparse(Diagonal(1 ./ grid_average(ϵ₀*d.ϵᵣ, ŷ)[:]));

        (Sxf, Sxb, Syf, Syb) = S_create(d.grid, ω);

        # Construct derivates
        δxb = Sxb*δ(x̂, Backward, d.grid);
        δxf = Sxf*δ(x̂, Forward,  d.grid);
        δyb = Syb*δ(ŷ, Backward, d.grid);
        δyf = Syf*δ(ŷ, Forward,  d.grid);

        if pol == TM
            # Construct system matrix
            A = δxf*μ₀^-1*δxb + δyf*μ₀^-1*δyb + ω^2*Tϵ;
            b = 1im*ω*d.src[:];

            ez = dolinearsolve(A, b, CSym)

            hx = -1/1im/ω/μ₀*δyb*ez;
            hy = 1/1im/ω/μ₀*δxb*ez;
            fields[i] = FieldTM(d.grid, ω, ez, hx, hy);
        elseif pol == TE
            # Construct system matrix
            A = δxf*Tϵxi*δxb + δyf*Tϵyi*δyb + ω^2*μ₀*speye(length(d.grid));
            b = 1im*ω*d.src[:];

            hz = dolinearsolve(A, b, CSym)

            ex = 1/1im/ω*Tϵyi*δyb*hz;
            ey = 1/1im/ω*Tϵxi*(-δxb*hz);

            fields[i] = FieldTE(d.grid, ω, hz, ex, ey);
        end
    end

    Nω == 1 && return fields[1]
    return fields
end


text:
# Test Parameter Learning routines

@testset "SPN parameter learning" begin

    import SPNetworks: NLL
    import SPNetworks.ParameterLearning: SEM, initialize, converged, update, backpropagate

    @testset "Simple DAG SPN" begin

        # # First create SPN to generate data from
        # SPN = SumProductNetwork(
        #     [
        #         SumNode([2,3,4],[0.2,0.5,0.3]),        # 1
        #         ProductNode([5,7]),                    # 2
        #         ProductNode([5,8]),                    # 3
        #         ProductNode([6,8]),                    # 4
        #         CategoricalDistribution(1, [0.6,0.4]), # 5
        #         CategoricalDistribution(1, [0.1,0.9]), # 6
        #         CategoricalDistribution(2, [0.3,0.7]), # 7
        #         CategoricalDistribution(2, [0.8,0.2])  # 8
        #     ]
        # )
        # # @info " Generating data set..."
        # N = 1000
        # data = rand(SPN,N)

        # Create SPN with indicator distributions at leaves
        SPN = SumProductNetwork(
            [
                SumNode([2,3,4],[0.2,0.5,0.3]),        # 1
                ProductNode([5,7]),                    # 2
                ProductNode([5,8]),                    # 3
                ProductNode([6,8]),                    # 4
                SumNode([9,10],[0.6,0.4]),             # 5
                SumNode([9,10],[0.1,0.9]),             # 6
                SumNode([11,12],[0.3,0.7]),            # 7
                SumNode([11,12],[0.8,0.2]),            # 8
                IndicatorFunction(1, 1.),              # 9
                IndicatorFunction(1, 2.),              # 10
                IndicatorFunction(2, 1.),              # 11
                IndicatorFunction(2, 2.)               # 12
            ]
        )
        # Generate dataset
        N = 3000
        data = rand(SPN,N)

        learner = SEM(SPN)
        #initialize(learner)
        # @info "Running Expectation Maximization until convergence..."
        while !converged(learner) && learner.steps < 10
            update(learner, data)
            # println("It: $(learner.steps) \t NLL: $(learner.score)")
            # println("It: $(learner.steps) \t NLL: $(learner.score) \t MAE: $error")
        end
        # # Reasonable values
        # @test learner.score > -14000 && learner.score < -12000

        # println()
        Z = 0.0
        @testset "Evaluation at $a,$b" for a in 1:2, b in 1:2
            est = SPN([a,b])
            # empirical distribution
            emp = sum(1.0 for i=1:N if data[i,1] == a && data[i,2] == b)/N
            # println("S($a,$b) = $ref ≈ ", est)
            @test est ≈ emp atol=0.05
            # @test ref ≈ est atol=0.01
            Z += est
        end
        @test Z ≈ 1.0

    end
    @testset "Breast-Cancer SPN" begin
        # load larger SPN
        SPN = SumProductNetwork(normpath("$(@__DIR__)/../assets/breast-cancer.spn"); offset = 1)
        #@show summary(SPN)
        # Generate dataset
        N = 3000
        data = rand(SPN, N)
        nll = NLL(SPN, data)
        learner = SEM(SPN)
        #initialize(learner)
        while !converged(learner) && learner.steps < 2
            update(learner, data)
            #println("It: $(learner.steps) \t NLL: $(learner.score)")
        end
        # # Training set NLL must be smaller than sampling distribution's training set NLL
        @test learner.score < nll
    end
    # @testset "NLTCS SPN" begin
    #     using DataFrames, CSV
    #     # load larger SPN
    #     SPN = SumProductNetwork(normpath("$(@__DIR__)/../assets/nltcs.spn"); offset=1)
    #     @show summary(SPN)
    #     # Load datasets
    #     tdata = convert(Matrix,
    #              DataFrame(CSV.File(normpath("$(@__DIR__)/../assets/nltcs.train.csv"),
    #                 header=collect(1:16) # columns names
    #     ))) .+ 1;
    #     @show summary(tdata)
    #     vdata = convert(Matrix,
    #                 DataFrame(CSV.File(normpath("$(@__DIR__)/../assets/nltcs.valid.csv"),
    #                   header=collect(1:16) # columns names
    #     ))) .+ 1;
    #     @show summary(vdata)
    #     # initialize EM learner
    #     learner = EMParamLearner()
    #     initialize(SPN)
    #     # Running Expectation Maximization
    #     while !converged(learner) && learner.steps < 10
    #         if learner.steps % 2 == 0
    #             tnll = NLL(SPN, vdata)
    #             update(learner, SPN, tdata)
    #             println("It: $(learner.steps) \t NLL: $(learner.score) \t held-out NLL: $tnll")
    #         else
    #             update(learner, SPN, tdata)
    #             println("It: $(learner.steps) \t NLL: $(learner.score)")
    #         end
    #     end

    # end
end # END of test set


text:
using StatsPlots, Distributions, CSV, DataFrames, Measures; pyplot()

realData = CSV.read("../data/temperatures.csv")

N = 10^5
include("../data/mvParams.jl")
biNorm = MvNormal(meanVect,covMat)
syntheticData = DataFrame(rand(MvNormal(meanVect,covMat),N)')
rename!(syntheticData,[:x1=>:Brisbane, :x2 => :GoldCoast])

default(c=cgrad([:blue, :red]),
    xlabel="Brisbane Temperature",
    ylabel="Gold Coast Temperature")

p1 = marginalhist(realData.Brisbane, realData.GoldCoast, bins=10:45)
p2 = marginalhist(syntheticData.Brisbane, syntheticData.GoldCoast, bins=10:.5:45)

plot(p1,p2, size = (1000,500), margin = 10mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# A functor for testing
struct Multiplier{T}
    x::T
end
(m::Multiplier)(y) = m.x * y
function ChainRulesCore.rrule(m::Multiplier, y)
    Multiplier_pullback(z̄) = Tangent{typeof(m)}(; x=y * z̄), m.x * z̄
    return m(y), Multiplier_pullback
end

# NoRules - has no rules defined
struct NoRules; end

"A function that outputs a vector from a scalar for testing"
make_two_vec(x) = [x, x]
function ChainRulesCore.rrule(::typeof(make_two_vec), x)
    make_two_vec_pullback(ȳ) = (NoTangent(), sum(ȳ))
    return make_two_vec(x), make_two_vec_pullback
end

@testset "test_helpers.jl" begin
    @testset "Multiplier functor" begin
        test_rrule(Multiplier(4.0), 3.0)
    end
    @testset "make_two_vec" begin
        test_rrule(make_two_vec, 1.5)
    end
end


text:
###############################################
## load the required function: FiltSparseMatrix
###############################################
#include("FiltSparseMatrix.jl")
#include("utils-de-error.jl");
include("FiltSparseMatrix.jl")
include("utils-de-error.jl");

using JLD
using DelimitedFiles ## for using "writedlm"

###########################################
## For Real-Data
###########################################
"""
"real_data_fvec_bettis_PIs" computes
the fvector, betti curves, and persistence intervals of a truncated filtration

Inputs:
(1) data_path: data source path
(2) d_up: dimension upper limit for persistent homology computation
(3) cutoff_size: row removal cutoff threshold
(4) cutoff_ind_ratio: filtration index cutoff

Output:
Nothing; but saves the results (and time peeks) in folders
"""
function real_data_fvec_bettis_PIs(data_path::String, d_up::Int,
        cutoff_size::Int, cutoff_ind_ratio::Float64)
    ## load the source data as a matrix M
    Data = load(data_path)
    M = Data["response"]

    ## the directory name of data source
    dir_data = join(split(data_path, "/")[1:end-1], "/")
    ## the directory name to save results; make directory if necessary
    dir_results = dir_data*"_results_real_0$(Int(cutoff_ind_ratio*100))"
    if !(isdir(dir_results))
        mkdir(dir_results)
    end

    ## the directory name to save progress peek; make directory if necessary
    dir_prog_peek = dir_data*"_prog_peek"
    if !(isdir(dir_prog_peek))
        mkdir(dir_prog_peek)
    end

    base_name = split(data_path, "/")[end][1:end-4] ## get the base name from data_path
    ctrl_name = "dup_$(d_up)_cutoff_$(cutoff_size)_0$(Int(cutoff_ind_ratio*100))"
    Progress_peek_name = dir_prog_peek*"/prog_peek_$(base_name)_$(ctrl_name)"

    ## apply the doubling trick
    M = Array{Float64,2}([M;-M])

    Time_record = [] ## for recording progress peek
    #########################
    ## compute the filtration
    #########################
    println("Working on filtration computation")
    push!(Time_record, "Working on filtration computation")
    writedlm(Progress_peek_name*".txt", Time_record)
    t_filt = @elapsed FS, ind_rem, true_birth_dic, denom = FiltSparseMatrix2(M, cutoff_size);
    push!(Time_record, t_filt)
    writedlm(Progress_peek_name*".txt", Time_record)

    #############################
    ## compute normalized indices
    #############################
    ind_norm = [true_birth_dic[k]/denom for k = 1:FS.birth[end]]
    ## compute the cutoff (un-normalized) index
    cutoff_ind = findlast(ind_norm.<=cutoff_ind_ratio)
    println("LOOK HERE (FOR DEBUG)")
    println(ind_norm[1:cutoff_ind])
    ## compute the index-truncated filtration
    birth_cut = FS.birth[FS.birth.<=cutoff_ind]
    faces_cut = FS.faces[FS.birth.<=cutoff_ind]
    FS_cut = FiltrationOfSimplicialComplexes(faces_cut, birth_cut);

    ####################################
    ## compute the skeleton and fvectors
    ####################################
    println("Working on Skeleton/f_vector computation")
    push!(Time_record, "Working on Skeleton/f_vector computation")
    writedlm(Progress_peek_name*".txt", Time_record)
    t_sk_fvec = @elapsed FS_Skeleton, f_vector = Skeleton_and_fvector(FS_cut, d_up)
    push!(Time_record, t_sk_fvec)
    writedlm(Progress_peek_name*".txt", Time_record)
    ## set the file name to save fvectors
    fvec_name = dir_results* "/$(base_name)_f_vectors_$(ctrl_name)"
    save(fvec_name*".jld",
        "f_vector", f_vector, "true_birth_dic", true_birth_dic, "denominator", denom)

    ####################################
    ## compute the persistence intervals
    ####################################
    println("Working on PI computation")
    push!(Time_record, "Working on PI computation")
    writedlm(Progress_peek_name*".txt", Time_record)
    t_PI = @elapsed PI = PersistenceIntervals(FS_Skeleton, Inf, "$(base_name)_$(ctrl_name)")
    push!(Time_record, t_PI)
    writedlm(Progress_peek_name*".txt", Time_record)
    ## save computed persistence intervals
    println("Saving computed PI")
    push!(Time_record, "Saving computed PI")
    ## set the file name to save PIs
    PIs_name = dir_results*"/$(base_name)_PIs_$(ctrl_name)"
    t_save_PI = @elapsed save(PIs_name*".jld",
        "PI", PI, "true_birth_dic", true_birth_dic, "denominator", denom)
    push!(Time_record, t_save_PI)
    writedlm(Progress_peek_name*".txt", Time_record)

    ###########################
    ## compute the betti curves
    ###########################
    println("Working on betti curves computation")
    push!(Time_record, "Working on betti curves computation")
    writedlm(Progress_peek_name*".txt", Time_record)
    PI_aug = [[zeros(2,2)]; PI]
    t_bettis = @elapsed bettis = Intervals2Bettis(PI_aug, FS_cut.birth[end])
    push!(Time_record, t_bettis)
    writedlm(Progress_peek_name*".txt", Time_record)
    ## set the file name to save betti curves
    bettis_name = dir_results*"/$(base_name)_bettis_$(ctrl_name)"
    save(bettis_name*".jld",
        "bettis", bettis, "true_birth_dic", true_birth_dic, "denominator", denom)
end

# ## Usage demonstration
# data_path = pwd()*"/flyData/fly_020_data.jld" ## data source path

# d_up = 5 ## dimension upper limit for persistent homology computation
# cutoff_size = 60 ## row removal cutoff threshold
# cutoff_ind_ratio = 0.65 ## filtration index cutoff

# real_data_fvec_bettis_PIs(data_path, d_up, cutoff_size, cutoff_ind_ratio)

###########################################
## For Controlled-Data
###########################################

include("unified_dir_pt_sample.jl") ## loading functions for generating random directions and points
"""
"ctrl_data_fvec_bettis_PIs" computes
the fvector, betti curves, and persistence intervals of a truncated filtration

Inputs:
(1) data_path: data source path
(2) d_up: dimension upper limit for persistent homology computation
(3) cutoff_size: row removal cutoff threshold
(4) cutoff_ind_ratio: filtration index cutoff
(5) dir_type: a string, the direction type (in ["pos", "unif"])
(6) pt_type: a string, the point type (in ["uniform", "sphere"])
(7) seed: an integer, for random number generation
(8) d: an integer, the dimension d where the controlled points and directions are in R^d

Output:
Nothing; but saves the results (and time peeks) in folders
"""
function ctrl_data_fvec_bettis_PIs(data_path::String, d_up::Int,
        cutoff_ind_ratio::Float64, cutoff_size::Int,
        dir_type::String, pt_type::String, seed::Int, d::Int,
        theta=pi/4::Float64, r=0.2::Float64)
    ## load the source data as a matrix M
    Data = load(data_path)
    M = Data["response"]

    ## the directory name of data source
    dir_data = join(split(data_path, "/")[1:end-1], "/")
    ## the directory name to save results; make directory if necessary
    dir_results = dir_data*"_results_ctrl_0$(Int(cutoff_ind_ratio*100))"
    if !(isdir(dir_results))
        mkdir(dir_results)
    end

    ## the directory name to save progress peek; make directory if necessary
    dir_prog_peek = dir_data*"_prog_peek"
    if !(isdir(dir_prog_peek))
        mkdir(dir_prog_peek)
    end

    base_name = split(data_path, "/")[end][1:end-4] ## get the base name from data_path
    ctrl_name = "dup_$(d_up)_dir_$(dir_type)_pt_$(pt_type)_d_$(d)_seed_$(seed)_cutoff_$(cutoff_size)_0$(Int(cutoff_ind_ratio*100))"
    if dir_type=="angle"
        name_angle = Int(round(theta*100))
        ctrl_name = "dup_$(d_up)_dir_$(dir_type)_$(name_angle)_pt_$(pt_type)_d_$(d)_seed_$(seed)_cutoff_$(cutoff_size)_0$(Int(cutoff_ind_ratio*100))"
    end
    Progress_peek_name = dir_prog_peek*"/prog_peek_$(base_name)_$(ctrl_name)"

    m, n = size(M)

    ## compute the number of nonzero entries in each row
    Rho = [length(findall(M[i,:].!=0)) for i = 1:m]
    directions = zeros(m,d)
    if dir_type=="angle"
        directions = dir_sample(dir_type::String, d::Int, m::Int, seed::Int, theta::Float64)
    else
        directions = dir_sample(dir_type::String, d::Int, m::Int, seed::Int)
    end
    pts = pt_sample(pt_type::String, d::Int, n::Int, seed::Int, r)

    ## generate the control matrix
    A_ctrl = directions*pts
    println(size(A_ctrl))
    for i = 1:m
        ## compute the column indices in row i with the top Rho_rem[i] values
        ind_keep_i = Array{Int,1}(sortslices([A_ctrl[i,:] collect(1:n)], dims = 1)[end-Rho[i]+1:end, 2])
        ## make the to-keep values all positive
        if ind_keep_i!=[]
            A_ctrl[i, ind_keep_i].-=(minimum(A_ctrl[i,ind_keep_i])-1)
        end
        ## make the to-remove values all zero
        ind_rm_i = setdiff(collect(1:n), ind_keep_i)
        #println("remove_$(i): ", ind_rm_i)
        A_ctrl[i, ind_rm_i] = zeros(n-Rho[i])
    end
    A_ctrl

    ## apply the doubling trick
    M_ctrl = [A_ctrl; -A_ctrl];

    Time_record = [] ## for progress peeking
    ## compute the filtration
    println("Working on filtration computation")
    push!(Time_record, "Working on filtration computation")
    writedlm(Progress_peek_name*".txt", Time_record)
    t_filt = @elapsed FS, ind_rem, true_birth_dic, denom = FiltSparseMatrix2(M_ctrl, cutoff_size);
    push!(Time_record, t_filt)
    writedlm(Progress_peek_name*".txt", Time_record)

    ## compute normalized indices
    ind_norm = [true_birth_dic[k]/denom for k = 1:FS.birth[end]]
    ## compute the cutoff (un-normalized) index
    cutoff_ind = findlast(ind_norm.<=cutoff_ind_ratio)

    birth_cut = FS.birth[FS.birth.<=cutoff_ind]
    faces_cut = FS.faces[FS.birth.<=cutoff_ind]
    FS_cut = FiltrationOfSimplicialComplexes(faces_cut, birth_cut);

    ## compute the skeleton and fvectors
    println("Working on Skeleton/f_vector computation")
    push!(Time_record, "Working on Skeleton/f_vector computation")
    writedlm(Progress_peek_name*".txt", Time_record)
    t_sk_fvec = @elapsed FS_Skeleton, f_vector = Skeleton_and_fvector(FS_cut, d_up)
    push!(Time_record, t_sk_fvec)
    writedlm(Progress_peek_name*".txt", Time_record)
    ## set the file name to save fvectors
    fvec_name = dir_results*"/$(base_name)_f_vectors_$(ctrl_name)"
    save(fvec_name*".jld",
        "f_vector", f_vector, "true_birth_dic", true_birth_dic, "denominator", denom)

    ## compute the persistence intervals
    println("Working on PI computation")
    push!(Time_record, "Working on PI computation")
    writedlm(Progress_peek_name*".txt", Time_record)
    t_PI = @elapsed PI = PersistenceIntervals(FS_Skeleton, Inf, "$(base_name)_$(ctrl_name)")
    push!(Time_record, t_PI)
    writedlm(Progress_peek_name*".txt", Time_record)
    ## save computed persistence intervals
    println("Saving computed PI")
    push!(Time_record, "Saving computed PI")
    ## set the file name to save PIs
    PIs_name = dir_results*"/$(base_name)_PIs_$(ctrl_name)"
    t_save_PI = @elapsed save(PIs_name*".jld",
        "PI", PI, "true_birth_dic", true_birth_dic, "denominator", denom)
    push!(Time_record, t_save_PI)
    writedlm(Progress_peek_name*".txt", Time_record)

    ## compute the betti curves
    println("Working on betti curves computation")
    push!(Time_record, "Working on betti curves computation")
    writedlm(Progress_peek_name*".txt", Time_record)
    PI_aug = [[zeros(2,2)]; PI]
    t_bettis = @elapsed bettis = Intervals2Bettis(PI_aug, FS_cut.birth[end])
    push!(Time_record, t_bettis)
    writedlm(Progress_peek_name*".txt", Time_record)
    ## set the file name to save betti curves
    bettis_name = dir_results*"/$(base_name)_bettis_$(ctrl_name)"
    save(bettis_name*".jld",
        "bettis", bettis, "true_birth_dic", true_birth_dic, "denominator", denom)
end

# ## Usage sample code for controlled data

# include("unified_real_ctrl_fvecs_bettis_PIs.jl");
# ## parameters
# data_path = pwd()*"/flyData/fly_020_data.jld" ## data source path

# d_up = 5 ## dimension upper limit for persistent homology computation
# cutoff_ind_ratio = 0.3 ## filtration index cutoff
# cutoff_size = 15 ## row removal cutoff threshold
# dir_type = "pos" ## direction control, choices in ["pos", "unif"]
# pt_type = "uniform" ## point control, choices in ["uniform", "sphere", "hyperbolic"]
# seed = 1
# d = 2

# ctrl_data_fvec_bettis_PIs(data_path::String, d_up::Int,
#         cutoff_ind_ratio::Float64, cutoff_size::Int,
#         dir_type::String, pt_type::String, seed::Int, d::Int)


text:
# License is MIT: https://github.com/JuliaString/LaTeX_Entities/LICENSE.md

__precompile__()

"""
# Public API (nothing is exported)

* lookupname(str)
* matchchar(char)
* matches(str)
* longestmatches(str)
* completions(str)
"""
module LaTeX_Entities

using StrTables

VER = UInt32(1)

struct LaTeX_Table{T} <: AbstractEntityTable
    ver::UInt32
    tim::String
    inf::String
    base32::UInt32
    base2c::UInt32
    nam::StrTable{T}
    ind::Vector{UInt16}
    val16::Vector{UInt16}
    ind16::Vector{UInt16}
    val32::Vector{UInt16}
    ind32::Vector{UInt16}
    val2c::Vector{UInt32}
    ind2c::Vector{UInt16}
end

function __init__()
    global default =
        LaTeX_Table(StrTables.load(joinpath(@__DIR__, "../data", "latex.dat"))...)
    nothing
end

end # module LaTeX_Entities


text:
import Base.sum

"""
    sum(N::AbstractEcologicalNetwork)

This function will return the sum of all interactions in the network. For
quantitative networks, this is the sum of interaction strengths. For binary
networks, this is the number of interactions. For probabilistic networks, this
is the expected number of realized interactions.
"""
function sum(N::AbstractEcologicalNetwork)
   return sum(N.A)
end

"""
    links(N::BinaryNetwork)

Number of non-zero interactions in a deterministic network.
"""
function links(N::BinaryNetwork)
    return sum(N)
end

"""
    links(N::QuantitativeNetwork)

Number of non-zero interactions in a quantitative network (use `sum` to get the
sum of interaction strengths).
"""
function links(N::QuantitativeNetwork)
    return sum(N.A .!= zero(eltype(N.A)))
end

"""
    links(N::ProbabilisticNetwork)

Expected number of interactions in a probabilistic network. To get the number of
interactions that have a non-zero probability, use *e.g.* `links(N>0.0)`.
"""
function links(N::ProbabilisticNetwork)
    return sum(N)
end

"""
**Variance in the expected number of links**

    links_var(N::ProbabilisticNetwork)

Expected variance of the number of links for a probabilistic network.
"""
function links_var(N::ProbabilisticNetwork)
   return sum(N.A .* (1 .- N.A))
end

"""
    connectance(N::AbstractEcologicalNetwork)

Number of links divided by the number of possible interactions. In unipartite
networks, this is ``L/S^2``. In bipartite networks, this is ``L/(T × B)``. It is
worth noting that while the maximal connectance is always 1 (i.e. the graph is
complete), the minimum value (assuming that the network is not degenerate) is
*not* 0. Instead, the minimum number of interactions in a unipartite network is
`S-1`, and in a bipartite network it is `max(T,B)`.

Connectance can therefore be transformed between 0 and 1, using the following
approach: let `m` be the minimum number of interactions, and Co be the measured
connectance, then the corrected value is `(Co-m)/(1-m)`. To our best knowledge,
this is not standard practice, and therefore is not suggested as a function in
the package.
"""
function connectance(N::AbstractEcologicalNetwork)
    return links(N) / (richness(N; dims=1)*richness(N; dims=2))
end

"""
    linkage_density(N::AbstractEcologicalNetwork)

Number of links divided by species richness.
"""
function linkage_density(N::AbstractEcologicalNetwork)
    return links(N) / richness(N)
end

"""
    connectance_var(N::ProbabilisticNetwork)

Expected variance of the connectance for a probabilistic matrix, measured as the
variance of the number of links divided by the squared size of the matrix.
"""
function connectance_var(N::ProbabilisticNetwork)
   return links_var(N) / (prod(size(N))^2)
end


text:
#===============================================================================

    McPeakOutput.jl

    melting curve analysis output struct
    for a single channel/well

    Author: Tom Price
    Date:   June 2019

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

#===============================================================================
    constants >>
===============================================================================#

const EMPTY_data    = zeros(Float_T, 1,3)[1:0,:]
const EMPTY_peaks   = zeros(Float_T, 1,2)[1:0,:]
const EMPTY_sn      = OrderedDict(:summit => EMPTY_peaks, :nadir => EMPTY_peaks)

#===============================================================================
    Field definitions >>
===============================================================================#

const PEAKOUTPUT_FIELD_DEFS = [
    Field(:observed_data,         Array{Float_T,2},               EMPTY_data),
    Field(:peaks_filtered,        Array{Float_T,2},               EMPTY_peaks),
    Field(:smoothed_data,         Array{Float_T,2},               EMPTY_data),
    Field(:negderiv_midrange,     Float_T,                        NaN_T),
    Field(:extremes,              OrderedDict{Symbol,Array{Float_T,2}},  EMPTY_sn),
    Field(:peaks_raw,             Array{Float_T,2},               EMPTY_peaks),
    Field(:peaks_reported,        Bool,                           false)]

#===============================================================================
    struct and constructor generation >>
===============================================================================#

abstract type McPeakOutput end

SCHEMA = PEAKOUTPUT_FIELD_DEFS
@make_struct_from_SCHEMA McPeakLongOutput McPeakOutput
@make_constructor_from_SCHEMA McPeakLongOutput

SCHEMA = PEAKOUTPUT_FIELD_DEFS[1:2]
@make_struct_from_SCHEMA McPeakShortOutput McPeakOutput
@make_constructor_from_SCHEMA McPeakShortOutput

#===============================================================================
    helper function >>
===============================================================================#

peak_output_format(out_format ::OutputFormat) =
    out_format == full_output ?
        McPeakLongOutput :
        McPeakShortOutput

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr