text:
# export  fit_cubic_spline,
#       sample_spline,
#       sample_spline_derivative,
#       sample_spline_derivative2,
#       sample_spline_speed,
#       sample_spline_theta,
#       sample_spline_curvature,
#       sample_spline_derivative_of_curvature,
#       calc_curve_length,
#       arclength,
#       calc_curve_param_given_arclen

function _integrate_simpsons(f::Function, a::Real, b::Real, n::Int)
    # integrate using Composite Simpson's rule
    # reference: https://en.wikipedia.org/wiki/Simpson%27s_rule

    @assert(n > 0) # number of intervals
    @assert(mod(n,2) == 0) # n must be even

    h = (b-a)/n
    retval = f(a) + f(b)
    flip = true
    for i = 1 : n-1
        retval += f(a+i*h) * (flip ? 4 : 2)
        flip = !flip
    end
    return h/3*retval
end

function _fit_open(pts::AbstractVector{Float64} )
    # fits the 1-D spline such that:
    #   spline goes through each point
    #   first and second derivatives match at each inner point
    #   the second derivative at the ends is zero
    # see: http://mathworld.wolfram.com/CubicSpline.html

    # this function returns a 4×(n-1) spline coefficient matrix, where n = |pts|

    n = length(pts)-1
    @assert(n > 0)

    M = spzeros(n+1,n+1)
    for i = 1 : n
        M[i,i] = 4
        M[i,i+1] = 1
        M[i+1,i] = 1
    end
    M[n+1,n+1] = 2
    M[1,1] = 2

    Y = Array{Float64}(n+1)
    for i = 1 : n+1
        ind_hi = min(i+1,n)
        ind_lo = max(1,i-1)
        Y[i] = 3*(pts[ind_hi] - pts[ind_lo])
    end

    D = M\Y

    spline_coeffs = Array{Float64}(4, n) # col is <a,b,c,d>
    spline_coeffs[1,:] = pts[1:n]
    spline_coeffs[2,:] = D[1:n]
    spline_coeffs[3,:] = 3*(pts[2:n+1] - pts[1:n]) -2*D[1:n]-D[2:n+1]
    spline_coeffs[4,:] = 2*(pts[1:n] - pts[2:n+1]) + D[1:n] + D[2:n+1]

    spline_coeffs
end
function _fit_closed(pts::AbstractVector{Float64} )
    # fits the 1-D spline such that:
    #   spline goes through each point
    #   first and second derivatives match at each inner point
    #   first the second derivative at the ends match
    # see: http://mathworld.wolfram.com/CubicSpline.html

    # this function returns a 4×n spline coefficient matrix, where n = |pts|

    n = length(pts)-1
    @assert(n > 0)

    M = spzeros(n+1,n+1)
    for i = 1 : n
        M[i,i] = 4
        M[i,i+1] = 1
        M[i+1,i] = 1
    end
    M[n+1,n+1] = 4
    M[1,n+1] = 1
    M[n+1,1] = 1

    Y = Array{Float64}(n+1)
    Y[1] = 3*(pts[2] - pts[n+1])
    for i = 2 : n
        Y[i] = 3*(pts[i+1] - pts[i-1])
    end
    Y[end] = 3*(pts[1] - pts[n])

    D = M\Y

    spline_coeffs = Array{Float64}(4, n+1) # col is <a,b,c,d>
    spline_coeffs[1,:] = pts
    spline_coeffs[2,:] = D
    spline_coeffs[3,1:n] = 3*(pts[2:n+1] - pts[1:n]) -2*D[1:n]-D[2:n+1]
    spline_coeffs[4,1:n] = 2*(pts[1:n] - pts[2:n+1]) + D[1:n] + D[2:n+1]
    spline_coeffs[3,n+1] = 3*(pts[1] - pts[n+1]) -2*D[n+1]-D[1]
    spline_coeffs[4,n+1] = 2*(pts[n+1] - pts[1]) + D[n+1] + D[1]

    spline_coeffs
end
function _fit_open(pts::Matrix{Float64}) # 2×n {x,y}

    # see http://mathworld.wolfram.com/CubicSpline.html

    d,n = size(pts)
    n -= 1

    Y = Array{Float64}(n+1)

    M = sparse(Int[], Int[], Float64[], n+1,n+1)
    for i in 1 : n
        M[i,i] = 4.0
        M[i,i+1] = 1.0
        M[i+1,i] = 1.0
    end
    M[n+1,n+1] = 2.0
    M[1,1] = 2.0

    retval = Array{Matrix{Float64}}(d)
    for k in 1 : d

        for i in 1 : n+1
            ind_hi = min(i+1,n)
            ind_lo = max(1,i-1)
            Y[i] = 3*(pts[k,ind_hi] - pts[k,ind_lo])
        end

        D = M \ Y

        spline_coeffs = Array{Float64}(4, n) # col is <a,b,c,d> for a + b⋅t + c⋅t² + d⋅t³
        spline_coeffs[1,:] = pts[k,1:n] # x₀
        spline_coeffs[2,:] = D[1:n]     # x'₀
        spline_coeffs[3,:] = 3*(pts[k,2:n+1]' - pts[k,1:n]') -2*D[1:n] - D[2:n+1] # -3x₀ + 3x₁ - 2x'₀ - x'₁
        spline_coeffs[4,:] = 2*(pts[k,1:n]' - pts[k,2:n+1]') +  D[1:n] + D[2:n+1] #  2x₀ - 2x₁ +  x'₀ + x'₁

        retval[k] = spline_coeffs
    end
    retval
end
function _fit_closed(pts::AbstractMatrix{Float64})
    d = size(pts,1)
    retval = Array{Matrix{Float64}}(d)
    for i = 1 : d
        retval[i] = _fit_closed(vec(pts[i,:]))
    end
    retval
end

function fit_cubic_spline(pts::AbstractArray{Float64}; open::Bool=true)
    if open
        return _fit_open(pts)
    else
        return _fit_closed(pts)
    end
end

function sample_spline(spline_coeffs::AbstractVector{Float64}, t::Float64)
    # here t is generally expected to be t ∈ [0,1]
    return spline_coeffs[1] + t*(spline_coeffs[2] + t*(spline_coeffs[3] + t*spline_coeffs[4]))
end
function sample_spline(spline_coeffs::AbstractMatrix{Float64}, t::Float64)
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    @assert(size(spline_coeffs, 1) == 4)
    col_ind = clamp(ceil(Int, t), 1, size(spline_coeffs,2))
    sample_spline(spline_coeffs[:,col_ind], t-col_ind+1)
end
function sample_spline(spline_coeffs::AbstractVector{Float64}, t_arr::AbstractVector{Float64})
    # here t is generally expected to be t ∈ [0,1]

    a = spline_coeffs[1]
    b = spline_coeffs[2]
    c = spline_coeffs[3]
    d = spline_coeffs[4]

    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        retval[i] = a + t*(b + t*(c + t*d))
    end
    retval
end
function sample_spline(spline_coeffs::AbstractMatrix{Float64}, t_arr::AbstractVector{Float64})
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    @assert(size(spline_coeffs, 1) == 4)
    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        col_ind = clamp(ceil(Int, t), 1, size(spline_coeffs,2))
        retval[i] = sample_spline(spline_coeffs[:,col_ind], t-col_ind+1)
    end
    retval
end

function sample_spline_derivative(spline_coeffs::AbstractVector{Float64}, t::Float64)
    # here t is generally expected to be t ∈ [0,1]
    return spline_coeffs[2] + t*(2spline_coeffs[3] + t*3spline_coeffs[4])
end
function sample_spline_derivative(spline_coeffs::AbstractMatrix{Float64}, t::Float64)
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    @assert(size(spline_coeffs, 1) == 4)
    col_ind = clamp(ceil(Int, t), 1, size(spline_coeffs,2))
    sample_spline_derivative(spline_coeffs[:,col_ind], t-col_ind+1)
end
function sample_spline_derivative(spline_coeffs::AbstractVector{Float64}, t_arr::AbstractVector{Float64})
    # here t is generally expected to be t ∈ [0,1]

    b = spline_coeffs[2]
    c = spline_coeffs[3]
    d = spline_coeffs[4]

    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        retval[i] = b + t*(2c + t*3d)
    end
    retval
end
function sample_spline_derivative(spline_coeffs::AbstractMatrix{Float64}, t_arr::AbstractVector{Float64})
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    @assert(size(spline_coeffs, 1) == 4)
    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        col_ind = clamp(ceil(Int, t), 1, size(spline_coeffs,2))
        retval[i] = sample_spline_derivative(spline_coeffs[:,col_ind], t-col_ind+1)
    end
    retval
end

function sample_spline_derivative2(spline_coeffs::AbstractVector{Float64}, t::Float64)
    # here t is generally expected to be t ∈ [0,1]
    return 2spline_coeffs[3] + t*6spline_coeffs[4]
end
function sample_spline_derivative2(spline_coeffs::AbstractMatrix{Float64}, t::Float64)
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    @assert(size(spline_coeffs, 1) == 4)
    col_ind = clamp(ceil(Int, t), 1, size(spline_coeffs,2))
    sample_spline_derivative2(spline_coeffs[:,col_ind], t-col_ind+1)
end
function sample_spline_derivative2(spline_coeffs::AbstractVector{Float64}, t_arr::AbstractVector{Float64})
    # here t is generally expected to be t ∈ [0,1]

    b = spline_coeffs[2]
    c = spline_coeffs[3]
    d = spline_coeffs[4]

    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        retval[i] = 2c + t*6d
    end
    retval
end
function sample_spline_derivative2(spline_coeffs::AbstractMatrix{Float64}, t_arr::AbstractVector{Float64})
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    @assert(size(spline_coeffs, 1) == 4)
    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        col_ind = clamp(ceil(Int, t), 1, size(spline_coeffs,2))
        retval[i] = sample_spline_derivative2(spline_coeffs[:,col_ind], t-col_ind+1)
    end
    retval
end

function sample_spline_speed(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64}, t::Float64)
    dxdt = sample_spline_derivative(spline_coeffs_x, t)
    dydt = sample_spline_derivative(spline_coeffs_y, t)
    hypot(dxdt, dydt)
end
function sample_spline_speed(spline_coeffs_x::AbstractMatrix{Float64}, spline_coeffs_y::AbstractMatrix{Float64}, t::Float64)
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x, 1) == 4)
    @assert(size(spline_coeffs_y, 1) == 4)
    @assert(n == size(spline_coeffs_y, 2))
    col_ind = clamp(ceil(Int, t), 1, n)::Int
    sample_spline_speed(spline_coeffs_x[:,col_ind], spline_coeffs_y[:,col_ind], t-col_ind+1)
end
function sample_spline_speed(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64}, t_arr::AbstractVector{Float64})
    # here t is generally expected to be t ∈ [0,1]

    bx = spline_coeffs_x[2]
    cx = spline_coeffs_x[3]
    dx = spline_coeffs_x[4]

    by = spline_coeffs_y[2]
    cy = spline_coeffs_y[3]
    dy = spline_coeffs_y[4]

    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        dxdt = bx + t*(2cx + t*3dx)
        dydt = by + t*(2cy + t*3dy)
        retval[i] = hypot(dxdt, dydt)
    end
    retval
end
function sample_spline_speed(spline_coeffs_x::AbstractMatrix{Float64}, spline_coeffs_y::AbstractMatrix{Float64}, t_arr::AbstractVector{Float64})
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.

    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x, 1) == 4)
    @assert(size(spline_coeffs_y, 1) == 4)
    @assert(n == size(spline_coeffs_y, 2))
    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        col_ind = clamp(ceil(Int, t), 1, n)
        retval[i] = sample_spline_speed(spline_coeffs_x[:,col_ind], spline_coeffs_y[:,col_ind], t-col_ind+1)
    end
    retval
end

function sample_spline_theta(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64}, t::Float64;
    stepsize=1e-4
    )

    # compute the angle from positive x-axis (counter-clockwise positive) of the curve in the positive t direction at t
    # uses an approximation via small step size instead of derivative due to zero-derivative issues
    # uses the forward derivative approximation unless it would put it out of range
    # result returned is in radians

    t_lo, t_hi = t, t+stepsize
    if t_hi > 1.0
        t_lo, t_hi = t-min(1000stepsize,0.1), t
    end

    x1 = sample_spline(spline_coeffs_x, t_lo)
    x2 = sample_spline(spline_coeffs_x, t_hi)
    y1 = sample_spline(spline_coeffs_y, t_lo)
    y2 = sample_spline(spline_coeffs_y, t_hi)

    # println("(t, lo, hi)  $t   $t_lo   $t_hi, ($(atan2(y2-y1, x2-x1)))")

    atan2(y2-y1, x2-x1)
end
function sample_spline_theta(spline_coeffs_x::AbstractMatrix{Float64}, spline_coeffs_y::AbstractMatrix{Float64}, t::Float64)
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x, 1) == 4)
    @assert(size(spline_coeffs_y, 1) == 4)
    @assert(n == size(spline_coeffs_y, 2))
    col_ind = clamp(ceil(Int, t), 1, n)
    sample_spline_theta(spline_coeffs_x[:,col_ind], spline_coeffs_y[:,col_ind], t-col_ind+1)
end
function sample_spline_theta(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64}, t_arr::AbstractVector{Float64})
    # here t is generally expected to be t ∈ [0,1]

    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        retval[i] = sample_spline_theta(spline_coeffs_x, spline_coeffs_y, t)
    end
    retval
end
function sample_spline_theta(spline_coeffs_x::AbstractMatrix{Float64}, spline_coeffs_y::AbstractMatrix{Float64}, t_arr::AbstractVector{Float64})
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.

    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x, 1) == 4)
    @assert(size(spline_coeffs_y, 1) == 4)
    @assert(n == size(spline_coeffs_y, 2))
    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        col_ind = clamp(ceil(Int, t), 1, n)
        retval[i] = sample_spline_theta(spline_coeffs_x[:,col_ind], spline_coeffs_y[:,col_ind], t-col_ind+1)
    end
    retval
end

function sample_spline_curvature(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64}, t::Float64)
    # computes the signed curvature

    dx  = sample_spline_derivative( spline_coeffs_x, t)
    dy  = sample_spline_derivative( spline_coeffs_y, t)
    ddx = sample_spline_derivative2(spline_coeffs_x, t)
    ddy = sample_spline_derivative2(spline_coeffs_y, t)

    (dx*ddy - dy*ddx)/(dx*dx + dy*dy)^1.5
end
function sample_spline_curvature(spline_coeffs_x::AbstractMatrix{Float64}, spline_coeffs_y::AbstractMatrix{Float64}, t::Float64)
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x, 1) == 4)
    @assert(size(spline_coeffs_y, 1) == 4)
    @assert(n == size(spline_coeffs_y, 2))
    col_ind = clamp(ceil(Int, t), 1, n)
    sample_spline_curvature(spline_coeffs_x[:,col_ind], spline_coeffs_y[:,col_ind], t-col_ind+1)
end
function sample_spline_curvature(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64}, t_arr::AbstractVector{Float64})
    # here t is generally expected to be t ∈ [0,1]

    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        retval[i] = sample_spline_curvature(spline_coeffs_x, spline_coeffs_y, t)
    end
    retval
end
function sample_spline_curvature(spline_coeffs_x::AbstractMatrix{Float64}, spline_coeffs_y::AbstractMatrix{Float64}, t_arr::AbstractVector{Float64})
    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.

    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x, 1) == 4)
    @assert(size(spline_coeffs_y, 1) == 4)
    @assert(n == size(spline_coeffs_y, 2))
    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        col_ind = clamp(ceil(Int, t), 1, n)
        retval[i] = sample_spline_curvature(spline_coeffs_x[:,col_ind], spline_coeffs_y[:,col_ind], t-col_ind+1)
    end
    retval
end

function sample_spline_derivative_of_curvature(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64}, t::Float64;
    stepsize=1e-4
    )

    # computes the derivative of the signed curvature

    t_lo, t_hi = t, t+stepsize
    if t_hi > 1.0
        t_lo, t_hi = t-stepsize, t
    end

    κ_hi = sample_spline_curvature(spline_coeffs_x, spline_coeffs_y, t_hi)
    κ_lo = sample_spline_curvature(spline_coeffs_x, spline_coeffs_y, t_lo)

    (κ_hi - κ_lo) / stepsize
end
function sample_spline_derivative_of_curvature(spline_coeffs_x::AbstractMatrix{Float64}, spline_coeffs_y::AbstractMatrix{Float64}, t::Float64;
    stepsize=1e-4
    )

    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.
    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x, 1) == 4)
    @assert(size(spline_coeffs_y, 1) == 4)
    @assert(n == size(spline_coeffs_y, 2))
    col_ind = clamp(ceil(Int, t), 1, n)
    sample_spline_derivative_of_curvature(spline_coeffs_x[:,col_ind], spline_coeffs_y[:,col_ind], t-col_ind+1, stepsize=stepsize)
end
function sample_spline_derivative_of_curvature(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64}, t_arr::AbstractVector{Float64};
    stepsize=1e-4
    )

    # here t is generally expected to be t ∈ [0,1]

    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        retval[i] = sample_spline_derivative_of_curvature(spline_coeffs_x, spline_coeffs_y, t, stepsize=stepsize)
    end
    retval
end
function sample_spline_derivative_of_curvature(spline_coeffs_x::AbstractMatrix{Float64}, spline_coeffs_y::AbstractMatrix{Float64}, t_arr::AbstractVector{Float64};
    stepsize=1e-4
    )

    # for t ∈ (-∞,1] we use spline_coeffs[:,1]
    # for t ∈ [1,2] we use spline_coeffs[:,2]
    # etc.

    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x, 1) == 4)
    @assert(size(spline_coeffs_y, 1) == 4)
    @assert(n == size(spline_coeffs_y, 2))
    retval = Array{Float64}(length(t_arr))
    for (i,t) in enumerate(t_arr)
        col_ind = clamp(ceil(Int, t), 1, n)
        retval[i] = sample_spline_derivative_of_curvature(spline_coeffs_x[:,col_ind], spline_coeffs_y[:,col_ind], t-col_ind+1, stepsize=stepsize)
    end
    retval
end

function calc_curve_length(spline_coeffs_x::AbstractVector{Float64}, spline_coeffs_y::AbstractVector{Float64};
    n_intervals::Int = 100
    )

    # integrate using Simpson's rule
    # _integrate_simpsons(t->sample_spline_speed(spline_coeffs_x, spline_coeffs_y, t), 0.0, 1.0, n_intervals)

    a = 0.0
    b = 1.0
    n = n_intervals

    h = (b-a)/n
    retval = sample_spline_speed(spline_coeffs_x, spline_coeffs_y, a) + sample_spline_speed(spline_coeffs_x, spline_coeffs_y, b)
    flip = true
    for i = 1 : n-1
        retval += sample_spline_speed(spline_coeffs_x, spline_coeffs_y, a+i*h) * (flip ? 4 : 2)
        flip = !flip
    end
    return h/3*retval
end
function calc_curve_length(
    spline_coeffs_x::AbstractMatrix{Float64},
    spline_coeffs_y::AbstractMatrix{Float64};
    n_intervals_per_segment::Int = 100
    )

    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_y, 2) == n)
    @assert(size(spline_coeffs_x, 1) == size(spline_coeffs_y, 1) == 4)

    len = 0.0
    for i = 1 : n
        len += calc_curve_length(spline_coeffs_x[:,i], spline_coeffs_y[:,i], n_intervals = n_intervals_per_segment)
    end
    len
end

function arclength(
    spline_coeffs_x::AbstractVector{Float64},
    spline_coeffs_y::AbstractVector{Float64},
    t_min::Real = 0.0,
    t_max::Real = 1.0,
    n_intervals::Int = 100
    )

    if isapprox(t_min, t_max)
        return 0.0
    end

    # _integrate_simpsons(t->sample_spline_speed(spline_coeffs_x, spline_coeffs_y, t), t_min, t_max, n_intervals)

    a = t_min
    b = t_max
    n = n_intervals

    h = (b-a)/n
    retval = sample_spline_speed(spline_coeffs_x, spline_coeffs_y, a) + sample_spline_speed(spline_coeffs_x, spline_coeffs_y, b)
    flip = true
    for i = 1 : n-1
        retval += sample_spline_speed(spline_coeffs_x, spline_coeffs_y, a+i*h) * (flip ? 4 : 2)
        flip = !flip
    end
    return h/3*retval
end
function arclength(
    spline_coeffs_x::AbstractMatrix{Float64},
    spline_coeffs_y::AbstractMatrix{Float64},
    t_min::Real = 0.0,
    t_max::Real = size(spline_coeffs_x, 2),
    n_intervals_per_segment::Int = 100
    )

    n = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_y, 2) == n)
    @assert(size(spline_coeffs_x, 1) == size(spline_coeffs_y, 1) == 4)

    if isapprox(t_min, t_max)
        return 0.0
    end

    # println("tmin/tmax: $t_min / $t_max")

    len = 0.0
    for i = floor(Int, t_min) : min(floor(Int, t_max), n-1)
        t_lo, t_hi = float(i), i+1.0

        spline_ind = i+1
        t_in_min = max(t_lo, t_min) - t_lo
        t_in_max = min(t_hi, t_max) - t_lo
        # println("($i) t_lo: $t_lo, t_hi: $t_hi, : $t_in_min → $t_in_max")
        len += arclength(spline_coeffs_x[:,spline_ind], spline_coeffs_y[:,spline_ind], t_in_min, t_in_max, n_intervals_per_segment)
    end
    # println("len: ", len)
    len
end

function calc_curve_param_given_arclen(
    spline_coeffs_x :: AbstractVector{Float64},
    spline_coeffs_y :: AbstractVector{Float64},
    s :: Float64;
    max_iterations :: Int=100,
    curve_length :: Float64 = calc_curve_length(spline_coeffs_x, spline_coeffs_y),
    epsilon::Float64 = 1e-4 # tolerance required before termination
    )

    # finds t such that p(t) is a distance s from start of curve
    # returns t=0 if s ≤ 0.0 and t=1 if s > L
    if s ≤ 0.0
        return 0.0
    elseif s ≥ curve_length
        return 1.0
    end

    t = s/curve_length
    lo, hi = 0.0, 1.0

    # @printf("%10s %10s %10s %10s %10s %10s\n", "iter", "lo", "hi", "t", "s", "F")
    # println("-"^65)

    for iter = 1 : max_iterations
        F = arclength(spline_coeffs_x, spline_coeffs_y, 0.0, t) - s

        # @printf("%10d %10.5f %10.5f %10.5f %10.5f %10.5f\n", iter-1, lo, hi, t, s, F)

        if abs(F) < epsilon
            # |F(t)| is close enough to zero, report it
            return t
        end

        DF = sample_spline_speed(spline_coeffs_x, spline_coeffs_y, t)
        tCandidate = t - F/DF
        if F > 0
            hi = t
            t = tCandidate ≤ lo ? 0.5*(lo+hi) : tCandidate
        else
            lo = t
            t = tCandidate ≥ hi ? 0.5*(lo+hi) : tCandidate
        end
    end

    # @printf("%10d %10.5f %10.5f %10.5f %10.5f %10s\n", max_iterations, lo, hi, t, s, "-")

    t
end
function calc_curve_param_given_arclen(
    spline_coeffs_x :: AbstractMatrix{Float64},
    spline_coeffs_y :: AbstractMatrix{Float64},
    s :: Float64;
    max_iterations :: Int=100,
    n_integration_intervals :: Int=100, # must be multiple of 2
    curve_length :: Float64 = calc_curve_length(spline_coeffs_x, spline_coeffs_y),
    epsilon::Float64 = 1e-4 # tolerance required before termination
    )

    # finds t such that p(t) is a distance s from start of curve
    # returns t=0 if s ≤ 0.0 and t=t_max if s > L

    n_segments = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x,1) == size(spline_coeffs_y,1) == 4)
    @assert(size(spline_coeffs_y,2) == n_segments)

    if s ≤ 0.0
        return 0.0
    elseif s ≥ curve_length
        return float(n_segments)
    end

    t = s/curve_length
    lo, hi = 0.0, float(n_segments)

    # @printf("%10s %10s %10s %10s %10s %10s\n", "iter", "lo", "hi", "t", "s", "F")
    # println("-"^65)

    for iter = 1 : max_iterations
        F = arclength(spline_coeffs_x, spline_coeffs_y, 0.0, t, n_integration_intervals) - s

        # @printf("%10d %10.5f %10.5f %10.5f %10.5f %10.5f\n", iter-1, lo, hi, t, s, F)

        if abs(F) < epsilon
            return t
        end

        DF = sample_spline_speed(spline_coeffs_x, spline_coeffs_y, t)
        tCandidate = t - F/DF
        if F > 0
            hi = t
            t = tCandidate ≤ lo ? 0.5*(lo+hi) : tCandidate
        else
            lo = t
            t = tCandidate ≥ hi ? 0.5*(lo+hi) : tCandidate
        end
    end

    # @printf("%10d %10.5f %10.5f %10.5f %10.5f %10s\n", max_iterations, lo, hi, t, s, "-")

    t
end
function calc_curve_param_given_arclen(
    spline_coeffs_x :: AbstractVector{Float64},
    spline_coeffs_y :: AbstractVector{Float64},
    s_arr :: AbstractVector{Float64}; # assumes s_arr is sorted
    max_iterations :: Int=100,
    curve_length :: Float64 = calc_curve_length(spline_coeffs_x, spline_coeffs_y),
    epsilon::Float64 = 1e-4 # tolerance required before termination
    )

    n = length(s_arr)
    t_arr = Array{Float64}(n)

    s = s_arr[1]
    t = s/curve_length
    if s ≤ 0.0
        t = 0.0
    elseif s ≥ curve_length
        t = 1.0
    end

    lo =  0.0

    for (i,s) in enumerate(s_arr)

        if s ≤ 0.0
            t = 0.0
            t_arr[i], lo = t, t
            continue
        elseif s ≥ curve_length
            t = 1.0
            t_arr[i], lo = t, t
            continue
        end

        hi = 1.0
        for iter = 1 : max_iterations
            F = arclength(spline_coeffs_x, spline_coeffs_y, 0.0, t) - s

            if abs(F) < epsilon
                t_arr[i], lo = t, t
                continue
            end

            DF = sample_spline_speed(spline_coeffs_x, spline_coeffs_y, t)
            tCandidate = t - F/DF
            if F > 0
                hi = t
                t = tCandidate ≤ lo ? 0.5*(lo+hi) : tCandidate
            else
                lo = t
                t = tCandidate ≥ hi ? 0.5*(lo+hi) : tCandidate
            end
        end

        t_arr[i], lo = t, t
    end

    t_arr
end
function calc_curve_param_given_arclen(
    spline_coeffs_x :: AbstractMatrix{Float64},
    spline_coeffs_y :: AbstractMatrix{Float64},
    s_arr :: AbstractVector{Float64}; # assumes s_arr is sorted
    max_iterations :: Int = 50,
    curve_length :: Float64 = calc_curve_length(spline_coeffs_x, spline_coeffs_y),
    epsilon::Float64 = 1e-4, # tolerance required before termination
    n_intervals_in_arclen::Int = 100
    )

    n_segments = size(spline_coeffs_x, 2)
    @assert(size(spline_coeffs_x,1) == size(spline_coeffs_y,1) == 4)
    @assert(size(spline_coeffs_y,2) == n_segments)

    n = length(s_arr)
    t_arr = Array{Float64}(n)

    s = s_arr[1]
    t = s/curve_length
    if s ≤ 0.0
        t = 0.0
    elseif s ≥ curve_length
        return float(n_segments)
    end

    lo =  0.0
    # println("L: ", curve_length)
    # println("s_max: ", s_arr[end])
    for (i,s) in enumerate(s_arr)

        # println("\ns: ", s)

        if s ≤ 0.0
            t = 0.0
            t_arr[i] = lo = t
            continue
        elseif s ≥ curve_length
            t = float(n_segments)
            t_arr[i] = lo = t
            continue
        end

        hi = float(n_segments)
        for iter = 1 : max_iterations
            F = arclength(spline_coeffs_x, spline_coeffs_y, 0.0, t, n_intervals_in_arclen) - s

            if abs(F) < epsilon
                break
            end

            DF = sample_spline_speed(spline_coeffs_x, spline_coeffs_y, t)
            tCandidate = t - F/DF
            if F > 0
                hi = t
                t = tCandidate ≤ lo ? 0.5*(lo+hi) : tCandidate
            else
                lo = t
                t = tCandidate ≥ hi ? 0.5*(lo+hi) : tCandidate
            end
        end

        t_arr[i] = lo = t
    end

    t_arr
end


text:
import Base.Operators.(*)

import Base: scale!, norm, vecdot

import Base: A_mul_B!, At_mul_B,  A_mul_Bt,  Ac_mul_B,  A_mul_Bc,  At_mul_Bt,  Ac_mul_Bc,  At_mul_Bt,
                       At_mul_B!, A_mul_Bt!, Ac_mul_B!, A_mul_Bc!, At_mul_Bt!, Ac_mul_Bc!, At_mul_Bt!

cublas_size(t::Char, M::CuVecOrMat) = (size(M, t=='N' ? 1:2), size(M, t=='N' ? 2:1))

###########
#
# BLAS 1
#
###########

#######
# SCAL
#######
scale!(x::CuArray{T}, k::Number) where {T<:CublasFloat} = CUBLAS.scal!(length(x), k, x, 1)

#######
# DOT
#######
function dot(x::CuVector{T}, rx::Union{UnitRange{TI},Range{TI}}, y::CuVector{T}, ry::Union{UnitRange{TI},Range{TI}}) where {T <: CublasFloat, TI<:Integer}
    if length(rx) != length(ry)
        throw(DimensionMismatch("length of rx, $(length(rx)), does not equal length of ry, $(length(ry))"))
    end
    if minimum(rx) < 1 || maximum(rx) > length(x)
        throw(BoundsError(x, rx))
    end
    if minimum(ry) < 1 || maximum(ry) > length(y)
        throw(BoundsError(y, ry))
    end
    dot(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end

At_mul_B(x::CuVector{T}, y::CuVector{T}) where {T<:CublasReal} = [CUBLAS.dot(x, y)]
At_mul_B(x::CuVector{T}, y::CuVector{T}) where {T<:CublasComplex} = [CUBLAS.dotu(x, y)]
Ac_mul_B(x::CuVector{T}, y::CuVector{T}) where {T<:CublasComplex} = [CUBLAS.dotc(x, y)]

vecdot(x::CuVector{T}, y::CuVector{T}) where {T<:CublasReal} = dot(x, y)
vecdot(x::CuVector{T}, y::CuVector{T}) where {T<:CublasComplex} = dotc(x, y)

#######
# NRM2
#######
norm(x::CuArray) = nrm2(x)

############
#
# BLAS 2
#
############

#########
# GEMV
##########
function gemv_wrapper!(y::CuVector{T}, tA::Char, A::CuMatrix{T}, x::CuVector{T},
                       alpha = one(T), beta = zero(T)) where T<:CublasFloat
    mA, nA = cublas_size(tA, A)
    if nA != length(x)
        throw(DimensionMismatch("second dimension of A, $nA, does not match length of x, $(length(x))"))
    end
    if mA != length(y)
        throw(DimensionMismatch("first dimension of A, $mA, does not match length of y, $(length(y))"))
    end
    if mA == 0
        return y
    end
    if nA == 0
        return scale!(y, 0)
    end
    gemv!(tA, alpha, A, x, beta, y)
end

A_mul_B!(y::CuVector{T}, A::CuMatrix{T}, x::CuVector{T}) where {T<:CublasFloat} = gemv_wrapper!(y, 'N', A,  x)
At_mul_B!(y::CuVector{T}, A::CuMatrix{T}, x::CuVector{T}) where {T<:CublasFloat} = gemv_wrapper!(y, 'T', A, x)
Ac_mul_B!(y::CuVector{T}, A::CuMatrix{T}, x::CuVector{T}) where {T<:CublasFloat} = gemv_wrapper!(y, 'T', A, x)
Ac_mul_B!(y::CuVector{T}, A::CuMatrix{T}, x::CuVector{T}) where {T<:CublasComplex} = gemv_wrapper!(y, 'C', A, x)

function (*)(A::CuMatrix{T}, x::CuVector{T}) where T<:CublasFloat
    A_mul_B!(similar(x, T, size(A,1)), A, x)
end

function At_mul_B(A::CuMatrix{T}, x::CuVector{T}) where T<:CublasFloat
    At_mul_B!(similar(x, T, size(A,2)), A, x)
end

function Ac_mul_B(A::CuMatrix{T}, x::CuVector{T}) where T<:CublasFloat
    Ac_mul_B!(similar(x, T, size(A,2)), A, x)
end

############
#
# BLAS 3
#
############

########
# GEMM
########
function gemm_wrapper!(C::CuVecOrMat{T}, tA::Char, tB::Char,
                   A::CuVecOrMat{T},
                   B::CuVecOrMat{T},
                   alpha = one(T),
                   beta = zero(T)) where T <: CublasFloat
    mA, nA = cublas_size(tA, A)
    mB, nB = cublas_size(tB, B)

    if nA != mB
        throw(DimensionMismatch("A has dimensions ($mA,$nA) but B has dimensions ($mB,$nB)"))
    end

    if C === A || B === C
        throw(ArgumentError("output matrix must not be aliased with input matrix"))
    end

    if mA == 0 || nA == 0 || nB == 0
        if size(C) != (mA, nB)
            throw(DimensionMismatch("C has dimensions $(size(C)), should have ($mA,$nB)"))
        end
        return scale!(C, 0)
    end

    gemm!(tA, tB, alpha, A, B, beta, C)
end

# Mutating
A_mul_B!(C::CuMatrix{T}, A::CuMatrix{T}, B::CuMatrix{T}) where {T <: CublasFloat} = gemm_wrapper!(C, 'N', 'N', A, B)
At_mul_B!(C::CuMatrix, A::CuMatrix, B::CuMatrix) = gemm_wrapper!(C, 'T', 'N', A, B)
A_mul_Bt!(C::CuMatrix, A::CuMatrix, B::CuMatrix) = gemm_wrapper!(C, 'N', 'T', A, B)
At_mul_Bt!(C::CuMatrix, A::CuMatrix, B::CuMatrix) = gemm_wrapper!(C, 'T', 'T', A, B)
Ac_mul_B!(C::CuMatrix{T}, A::CuMatrix{T}, B::CuMatrix{T}) where {T<:CublasReal} = At_mul_B!(C, A, B)
Ac_mul_B!(C::CuMatrix, A::CuMatrix, B::CuMatrix) = gemm_wrapper!(C, 'C', 'N', A, B)

function A_mul_B!(C::CuMatrix{T}, A::CuVecOrMat{T}, B::CuVecOrMat{T}) where T
    gemm_wrapper!(C, 'N', 'N', A, B)
end

# Non mutating

# A_mul_Bx
function (*)(A::CuMatrix{T}, B::CuMatrix{T}) where T <: CublasFloat
    A_mul_B!(similar(B, T,(size(A,1), size(B,2))), A, B)
end

function A_mul_Bt(A::CuMatrix{T}, B::CuMatrix{T}) where T
    A_mul_Bt!(similar(B, T, (size(A,1), size(B,1))), A, B)
end

function A_mul_Bc(A::CuMatrix{T}, B::CuMatrix{T}) where T
    A_mul_Bc!(similar(B, T,(size(A,1),size(B,1))),A, B)
end

# At_mul_Bx
function At_mul_B(A::CuMatrix{T}, B::CuMatrix{T}) where T
    At_mul_B!(similar(B, T, (size(A,2), size(B,2))), A, B)
end

function At_mul_Bt(A::CuMatrix{T}, B::CuMatrix{T}) where T
    At_mul_Bt!(similar(B, T, (size(A,2), size(B,1))), A, B)
end

# Ac_mul_Bx
function Ac_mul_B(A::CuMatrix{T}, B::CuMatrix{T}) where T
    Ac_mul_B!(similar(B, T, (size(A,2), size(B,2))), A, B)
end

function Ac_mul_Bt(A::CuMatrix{T}, B::CuMatrix{S}) where {T,S}
    Ac_mul_Bt(similar(B, T, (size(A,2), size(B,1))), A, B)
end

function Ac_mul_Bc(A::CuMatrix{T}, B::CuMatrix{S}) where {T,S}
    Ac_mul_Bc!(similar(B, T, (size(A,2), size(B,1))), A, B)
end

ext: jl


text:
Base.@kwdef struct Platform
    arch::String
    rootfs_tag::String
    rootfs_treehash::String
    allow_fail::Bool
    commit_status::Bool
end

struct Platforms
    ps::Vector{Platform}
end

Base.length(platforms::Platforms)         = Base.length(platforms.ps)
Base.iterate(platforms::Platforms)        = Base.iterate(platforms.ps)
Base.iterate(platforms::Platforms, state) = Base.iterate(platforms.ps, state)

const platforms = Platforms(
    [
        Platform(;
            arch            = "x86_64",
            rootfs_tag      = "v5.14",
            rootfs_treehash = "72dab0734a234d2277e3a204263ec7d038b4087e",
            allow_fail      = false,
            commit_status   = true,
        ),
        Platform(;
            arch            = "aarch64",
            rootfs_tag      = "v5.14",
            rootfs_treehash = "20db175a1cc643f6d05f172eae921171e95326d3",
            allow_fail      = true,
            commit_status   = false,
        ),
    ]
)

max_stars_repo_name: mozilla/rr

max_issues_repo_name: mozilla/rr

max_forks_repo_name: mozilla/rr


text:
using Test
using TightBindingLattice

@testset "irrep" begin
    unitcell = make_unitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
    addsite!(unitcell, "Ox", FractCoord([0,0], [0.5, 0.0]))
    addsite!(unitcell, "Oy", FractCoord([0,0], [0.0, 0.5]))

    lattice = make_lattice(unitcell, [4 0; 0 4])
    tsym = TranslationSymmetry(lattice)
    psym = project(PointSymmetryDatabase.get(13), [1 0 0; 0 1 0])

    tsymbed = embed(lattice, tsym)
    psymbed = embed(lattice, psym)

    @testset "translation" begin
        @test_throws ArgumentError IrrepComponent(tsym, 99999)
        @test_throws ArgumentError IrrepComponent(tsym, 1, 2)
        @test IrrepComponent(tsym, 1, 1) == IrrepComponent(tsym, 1, 1)
        @test IrrepComponent(tsym, 1, 1) != IrrepComponent(tsym, 2, 1)

        for tsym_irrep in 1:num_irreps(tsym)
            tsic = IrrepComponent(tsym, tsym_irrep)
            tsicbed = IrrepComponent(tsymbed, tsym_irrep)

            @test group_order(tsic) == group_order(tsym)
            psym_little = little_symmetry(tsic, psym)
            psymbed_little = little_symmetry(tsicbed, psymbed)
            @test psym_little.group == symmetry(psymbed_little).group

            k = tsym.coordinates[tsym_irrep]
            @test iscompatible(tsym, tsym_irrep, psym) == (k in [[0,0], [2,2]])
            @test iscompatible(tsym, tsym_irrep, psym_little)

            @test iscompatible(tsic, psym) == (k in [[0,0], [2,2]])
            @test iscompatible(tsic, psym_little)

            @test iscompatible(tsicbed, psymbed) == (k in [[0,0], [2,2]])
            @test iscompatible(tsicbed, psymbed_little)

            @test little_group(tsic, psym) == little_group(tsicbed, psymbed)

            lge1 = little_group_elements(tsic, psym)
            lge2 = little_group_elements(tsicbed, psymbed)
            @test lge1 == lge2
            lg_matrep = psym.matrix_representations[lge1]
            @test !isnothing(group_isomorphism(little_group(tsym, tsym_irrep, psym),
                                  FiniteGroup(group_multiplication_table(lg_matrep))))
            let psic = IrrepComponent(psym, 1, 1)
                if k in [[0,0], [2,2]]
                    SymmorphicIrrepComponent(tsic, psic)
                else
                    @test_throws ArgumentError SymmorphicIrrepComponent(tsic, psic)
                end
            end
        end # for tsym_irrep

        for tsic in get_irrep_components(tsym)
            k = tsym.coordinates[tsic.irrep_index]
            psym_little = little_symmetry(tsic, psym)

            @test iscompatible(tsic, psym) == (k in [[0,0], [2,2]])
            @test iscompatible(tsic, psym_little)
            lg_matrep = psym.matrix_representations[little_group_elements(tsic, psym)]
            @test !isnothing(group_isomorphism(little_group(tsic, psym),
                                FiniteGroup(group_multiplication_table(lg_matrep))))
            let psic = IrrepComponent(psym, 1, 1)
                if k in [[0,0], [2,2]]
                    SymmorphicIrrepComponent(tsic, psic)
                else
                    @test_throws ArgumentError SymmorphicIrrepComponent(tsic, psic)
                end
            end
        end
    end

    @testset "point" begin
        @test_throws ArgumentError IrrepComponent(psym, 99999, 1)
        @test_throws ArgumentError IrrepComponent(psym, 1, 10)
        @test length(collect(get_irrep_components(psym))) == 1 + 1 + 1 + 1 + 2
        @test IrrepComponent(psym, 1, 1) == IrrepComponent(psym, 1, 1)
        @test IrrepComponent(psym, 1, 1) != IrrepComponent(psym, 2, 1)

        for psic in get_irrep_components(psym)
            @test group_order(psic) == group_order(psym)

            permphase = collect(get_irrep_iterator(psic))
            @test length(permphase) == group_order(psym)

            perms = [p for (p, ϕ) in permphase]
            @test ishomomorphic(psym.group, perms)

            mtab = group_multiplication_table(perms)
            @test mtab == psym.group.multiplication_table
        end
    end

    @testset "space" begin
        ssics1 = []
        for tsic in get_irrep_components(tsym)
            psym_little = little_symmetry(tsic, psym)
            for psic in get_irrep_components(psym_little)
                push!(ssics1, (tsic, psic))
            end
        end
        ssics2 = []
        for ssic in get_irrep_components(tsym ⋊ psym)
            push!(ssics2, ssic)
        end
        @test length(ssics1) == length(ssics2)
        for ((tsic, psic), ssic) in zip(ssics1, ssics2)
            @test group_order(ssic) == group_order(tsic) * group_order(psic)

            @test tsic.symmetry.orthocube == ssic.normal.symmetry.orthocube
            @test psic.symmetry.hermann_mauguin == ssic.rest.symmetry.hermann_mauguin
            @test tsic.irrep_index == ssic.normal.irrep_index
            @test psic.irrep_index == ssic.rest.irrep_index
            @test tsic.irrep_component == ssic.normal.irrep_component
            @test psic.irrep_component == ssic.rest.irrep_component

            @test sum(1 for x in collect(get_irrep_iterator(ssic))) == group_order(tsic) * group_order(psic)
        end
    end

end # testset little_symmetry


text:
using BasisFunctions: UnitInterval, PolynomialBasis
using BasisFunctions: hasderivative, hasantiderivative, ordering, support
using BasisFunctions: supremum, infimum
using BasisFunctions.Test: suitable_function

const PBSplineInterval = UnitInterval
const PBSplineIndex = NativeIndex{:pbspline}

# standard grid generation
get_pbspline_nodes(T,n) = collect(LinRange{T}(0.0, 1.0, n))

function check_equispaced(ξ::Vector{T}) where {T}
    Δ = ξ[2] - ξ[1]
    for i in 1:(length(ξ)-1)
        abs(ξ[i+1] - ξ[i] - Δ) > 1/length(ξ) * 1e-2 ? (return false) : nothing
    end
    true
end

# constructors
struct PBSpline{T} <: PolynomialBasis{T}
    p :: Int                    # spline degree
    nodes :: ScatteredGrid{T}   # knot points
    # n :: Int                    # number of splines
    # L :: T                      # domain length
    # h :: T                      # L/n

    function PBSpline(p::Int, nodes::ScatteredGrid{T}) where {T}
        @assert length(nodes) > p
        @assert check_equispaced(nodes.points)

        new{T}(p,nodes)
    end
end

# constructors for unit interval grid
PBSpline(p::Int, n::Int) = PBSpline(p, get_pbspline_nodes(Float64,n))
PBSpline(p::Int, ξ::Vector{T}) where {T} = PBSpline(p, ScatteredGrid(ξ, PBSplineInterval{T}()))

# n splines of degree p=n-1
PBSpline(n::Int) = PBSpline(n-1, get_pbspline_nodes(Float64,n))
PBSpline{T}(n::Int) where {T} = PBSpline(n-1, get_pbspline_nodes(T,n))
PBSpline(ξ::Vector{T}) where {T} = PBSpline(length(ξ)-1, ξ)

# constructor for interval [a,b)
PBSpline(p, ξ, a::Number, b::Number) = rescale(PBSpline(p, ξ), a, b)

# basis properties
nodes(b::PBSpline)  = b.nodes.points
nnodes(b::PBSpline) = length(b.nodes.points)
degree(b::PBSpline) = b.p
Base.size(b::PBSpline) = (nnodes(b),)

# convenience
L_domain(b::PBSpline) = nodes(b)[end] - nodes(b)[1]
invh_elements(b::PBSpline) = nnodes(b)/L_domain(b)

BasisFunctions.native_index(b::PBSpline, idx) = PBSplineIndex(idx)
BasisFunctions.ordering(b::PBSpline) = Base.OneTo(nnodes(b))
BasisFunctions.period(b::PBSpline{T}) where {T} = one(T)
BasisFunctions.support(b::PBSpline{T}) where {T} = PBSplineInterval{T}()

BasisFunctions.interpolation_grid(b::PBSpline{T}) where {T} = PeriodicEquispacedGrid(length(b), PBSplineInterval{T}())
BasisFunctions.hasinterpolationgrid(b::PBSpline) = true

BasisFunctions.similar(::PBSpline, ::Type{T}, n::Int) where {T} = PBSpline{T}(n)
BasisFunctions.similar(::PBSpline, ::Type{T}, ξ::Vector{T}) where {T} = PBSpline(ξ)

# for testing purposes
BasisFunctions.Test.suitable_function(s::PBSpline) = x->exp(cos(2π*x/(supremum(support(s)) - infimum(support(s)))))

# evaluation at point x
function _pbspline(b::PBSpline{T}, i::PBSplineIndex, x::T,
                    n::Int = nnodes(b), invh::T = invh_elements(b)) where {T}
    x -= nodes(b)[1]; x *= invh; x -= (i-1)
    # periodicity
    while x < zero(T)  x += n end
    while x > n        x -= n end
    return _unitpbspline(b.p,x)
end

_pbspline(b::PBSpline, i, x) = _pbspline(b, native_index(b,i), x)

# evaluation of "unit" spline with support (0,n+1)
function _unitpbspline(p::Int, x::T) where {T}
    if p == 0
        return zero(T) <= x < one(T) ? one(T) : zero(T)
    elseif x > p+1
        return zero(T)
    else
        return x/p * _unitpbspline(p-1,x) + (p+1-x)/p * _unitpbspline(p-1,x-1)
    end
end

function BasisFunctions.unsafe_eval_element(b::PBSpline, i, x)
    _pbspline(b, i, x)
end

function _pbspline_derivative(b::PBSpline{T}, i::PBSplineIndex, x::T,
                    n::Int = nnodes(b), invh::T = invh_elements(b)) where {T}
    x -= nodes(b)[1]; x *= invh; x -= (i-1)
    # periodicity
    while x < zero(T)  x += n end
    while x > n        x -= n end
    invh * (_unitpbspline(b.p-1, x) - _unitpbspline(b.p-1, x-1 < 0 ? x-1+n : x-1) )
end

function BasisFunctions.unsafe_eval_element_derivative(b::PBSpline, i, x, nᵢ::Int = nnodes(b))
    _pbspline_derivative(b, i, x)
end


text:
__precompile__()
module Random
using Compat #xor
export ParisiRapuano, randperm!, getRNG,
        randchoice, randchoice!, randint,
        DiscreteDistribution

using Base.Random
import Base: rand, srand

include("parisi_rapuano.jl")
include("discrete_distribution.jl")

# Knuth shuffles (see Wikipedia)
"""
    randperm!([rng,] v::AbstractArray)

In-place random permutation of the elements of `v`.
"""
function randperm!(rng::AbstractRNG, v::AbstractArray)
    n = length(v)
    for i=1:n-1
        r = rand(rng, 0:n-i)
        v[i], v[i+r] = v[i+r], v[i]
    end
end

randperm!(v::AbstractArray) = randperm!(GLOBAL_RNG, v)

"""
randchoice!([rng,] a, k; exclude=nothing)

Sample `k` elements from array `a` without repetition and eventually excluding elements in `exclude`.
Pay attention, it has complecity `O(length(a))`, since it creates a copy of `a`
internally.
Use [`randchoice!`](@ref) if you want a `O(k)` sampler.
"""
randchoice(rng::AbstractRNG, a::AbstractArray{T}, k::Integer; exclude = T[]) where {T} = randchoice!(rng, collect(a), k; exclude = exclude)
randchoice(a::AbstractArray{T}, k::Integer; exclude = T[]) where {T} = randchoice(GLOBAL_RNG, a, k; exclude = exclude)

"""
randchoice!([rng,] a, k; exclude=nothing)

Sample `k` elements from array `a` without repetition and eventually excluding elements in `exclude`.
Pay attention, it changes the order of the elements in `a`.
Use [`randchoice`](@ref) for the non-mutating version.
"""
function randchoice!(rng::AbstractRNG, a::AbstractArray{T}, k::Integer; exclude::Vector{T} = T[]) where T
    length(a) < k + length(exclude) && error("Array too short.")
    res = Vector{eltype(a)}()
    sizehint!(res, k)
    n = length(a)
    i = 1
    while length(res) < k
        r = rand(rng, 1:n-i+1)
        if !(a[r] in exclude)
            push!(res, a[r])
            a[r], a[n-i+1] = a[n-i+1], a[r]
            i += 1
        end
    end
    return res
end

randchoice!(a::AbstractArray{T}, k::Integer; exclude = T[]) where {T} = randchoice!(GLOBAL_RNG, a, k; exclude = exclude)

"""
    randchoice([rng,] itr)

Returns a random element of `itr`. This implementation
is more efficient than the one in Base for `rand(itr)`.
"""
randchoice(rng, itr) = nth(itr, _randchoice(rng, length(itr)))
randchoice(itr) = randchoice(GLOBAL_RNG, itr)

_randchoice(rng, n::Int) = floor(Int, rand(rng) * n) + 1

#
#    nth(xs, n::Integer)
#
#Return the n'th element of xs. Mostly useful for non indexable collections.
function nth(xs, n::Integer)
    n > 0 || throw(BoundsError(xs, n))
    # catch, if possible
    applicable(length, xs) && (n ≤ length(xs) || throw(BoundsError(xs, n)))
    s = start(xs)
    i = 0
    while !done(xs, s)
        (val, s) = next(xs, s)
        i += 1
        i == n && return val
    end
    # catch iterators with no length but actual finite size less then n
    throw(BoundsError(xs, n))
end

nth(xs::AbstractArray, n::Integer) = xs[n]

"""
    getRNG(seed)

Returns `MersenneTwister(seed)` if `seed` is greater than `0`, `GLOBAL_RNG`
otherwise.
"""
getRNG(seed::Integer) = seed > 0 ? MersenneTwister(seed) : GLOBAL_RNG

end #submodule

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
module trlib
    type trlib_data
        H
        invM
        grad          :: Array{Float64, 1}
        itmax         :: Int64
        equality      :: Int64
        itmax_lanczos :: Int64
        tol_rel_i     :: Float64
        tol_abs_i     :: Float64
        tol_rel_b     :: Float64
        tol_abs_b     :: Float64
        zero          :: Float64
        obj_lb        :: Float64
        ctl_invariant :: Int64
        convexify     :: Int64
        earlyterm     :: Int64
        refine        :: Int64
        verbose       :: Int64
        h_pointer     :: Int64
        iwork         :: Array{Int64, 1}
        fwork         :: Array{Float64, 1}
        timing        :: Array{Int64, 1}
        init          :: Int64
        sol           :: Array{Float64, 1}
        g             :: Array{Float64, 1}
        v             :: Array{Float64, 1}
        gm            :: Array{Float64, 1}
        p             :: Array{Float64, 1}
        Hp            :: Array{Float64, 1}
        Q             :: Array{Float64, 2}
        ret           :: Int64
        obj           :: Float64
        lam           :: Float64
        trlib_data(H, invM, grad) = (itmax = min(round(Int, 2e8/size(grad,1)), 2*size(grad, 1));
                               iwork_size = Array{Clong}(1);
                               fwork_size = Array{Clong}(1);
                               h_pointer = Array{Clong}(1);
                               ccall(
                                     (:trlib_krylov_memory_size, "libtrlib"),
                                     Clong,
                                     (Clong, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}),
                                     itmax, iwork_size, fwork_size, h_pointer);
                               fwork = Array{Cdouble}(fwork_size[1]);
                               ccall(
                                     (:trlib_krylov_prepare_memory, "libtrlib"),
                                     Clong,
                                     (Clong, Ptr{Cdouble}),
                                     itmax, fwork);
                               new(H, invM, grad, itmax,
                                   0, 100, -2.0, 0.0, -3.0, 0.0, 2e-16, -1e20, 0, 1, 1, 1, 0,
                                   h_pointer[1], Array{Clong}(iwork_size[1]), fwork,
                                   Array{Clong}(ccall((:trlib_krylov_timing_size, "libtrlib"), Clong, ())),
                                   1, similar(grad), similar(grad), similar(grad),
                                   similar(grad), similar(grad), similar(grad),
                                   Array{Float64}(size(grad,1), itmax+1), 0, 0.0, 0.0 ) )
        trlib_data(H, grad) = (itmax = min(round(Int, 2e8/size(grad,1)), 2*size(grad, 1));
                               iwork_size = Array{Clong}(1);
                               fwork_size = Array{Clong}(1);
                               h_pointer = Array{Clong}(1);
                               ccall(
                                     (:trlib_krylov_memory_size, "libtrlib"),
                                     Clong,
                                     (Clong, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}),
                                     itmax, iwork_size, fwork_size, h_pointer);
                               fwork = Array{Cdouble}(fwork_size[1]);
                               ccall(
                                     (:trlib_krylov_prepare_memory, "libtrlib"),
                                     Clong,
                                     (Clong, Ptr{Cdouble}),
                                     itmax, fwork);
                               new(H, speye(size(grad,1)), grad, itmax,
                                   0, 100, -2.0, 0.0, -3.0, 0.0, 2e-16, -1e20, 0, 1, 1, 1, 0,
                                   h_pointer[1], Array{Clong}(iwork_size[1]), fwork,
                                   Array{Clong}(ccall((:trlib_krylov_timing_size, "libtrlib"), Clong, ())),
                                   1, similar(grad), similar(grad), similar(grad),
                                   similar(grad), similar(grad), similar(grad),
                                   Array{Float64}(size(grad,1), itmax+1), 0, 0.0, 0.0 ) )
    end

    function trlib_solve(TR::trlib_data, radius::Float64)
        g_dot_g = 0.0
        v_dot_g = 0.0
        p_dot_Hp = 0.0

        action = Array{Clong}(1)
        it     = Array{Clong}(1)
        ityp   = Array{Clong}(1)
        flt1   = Array{Cdouble}(1)
        flt2   = Array{Cdouble}(1)
        flt3   = Array{Cdouble}(1)

        while true
            TR.ret = ccall(
                          (:trlib_krylov_min, "libtrlib"),
                          Clong,
                          (Clong, Cdouble, Clong, Clong, Clong, Cdouble, Cdouble,
                           Cdouble, Cdouble, Cdouble, Cdouble, Clong, Clong, Clong,
                           Cdouble, Cdouble, Cdouble, Ptr{Clong}, Ptr{Cdouble},
                           Clong, Clong, Clong, Cstring, Ptr{Clong}, Ptr{Clong},
                           Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Cdouble},
                           Ptr{Cdouble}, Ptr{Cdouble}),
                          TR.init, radius, TR.equality, TR.itmax, TR.itmax_lanczos,
                          TR.tol_rel_i, TR.tol_abs_i, TR.tol_rel_b, TR.tol_abs_b,
                          TR.zero, TR.obj_lb, TR.ctl_invariant, TR.convexify,
                          TR.earlyterm, g_dot_g, v_dot_g, p_dot_Hp, TR.iwork,
                          TR.fwork, TR.refine, 0, 0, "", C_NULL,
                          TR.timing, action, it, ityp, flt1, flt2, flt3)
            TR.init = 0
            if action[1] == 1
                TR.sol[:] = 0.0
                TR.gm[:] = 0.0
                TR.g = TR.grad
                TR.v = TR.invM * TR.g
                g_dot_g = dot(TR.g, TR.g)
                v_dot_g = dot(TR.v, TR.g)
                TR.p = -TR.v
                TR.Hp = TR.H*TR.p
                p_dot_Hp = dot(TR.p, TR.Hp)
                TR.Q[:,1] = TR.v/sqrt(v_dot_g)
            end
            if action[1] == 2
                TR.sol = TR.Q[:,1:1+it[1]] * TR.fwork[1+TR.h_pointer:1+TR.h_pointer+it[1]]
            end
            if action[1] == 3
                if ityp[1] == 1
                    TR.sol += flt1[1] * TR.p
                end
            end
            if action[1] == 4
                if ityp[1] == 1
                    TR.Q[:,1+it] = flt2[1] * TR.v
                    TR.gm = TR.g
                    TR.g += flt1[1]*TR.Hp
                end
                if ityp[1] == 2
                    TR.sol = TR.Hp + flt1[1] * TR.g + flt2[1] * TR.gm
                    TR.gm = flt3[1] * TR.g
                    TR.g = TR.sol
                end
                TR.v = TR.invM * TR.g
                g_dot_g = dot(TR.g, TR.g)
                v_dot_g = dot(TR.v, TR.g)
            end
            if action[1] == 5
                TR.p = flt1[1]*TR.v + flt2[1]*TR.p
                TR.Hp = TR.H*TR.p
                p_dot_Hp = dot(TR.p, TR.Hp)
                if ityp[1] == 2
                    TR.Q[:,1+it] = TR.p
                end
            end
            if action[1] == 8
                g_dot_g = .5 * dot(TR.sol, TR.H*TR.sol) + dot(TR.sol, TR.grad)
            end
            if TR.ret < 10
                break
            end
        end
        TR.obj = TR.fwork[9]
        TR.lam = TR.fwork[8]
        TR.init = 2
    end
end

max_stars_repo_name: JonasZehn/trlib

max_issues_repo_name: JonasZehn/trlib

max_forks_repo_name: JonasZehn/trlib


text:
using JuMP, GLPK

# ---------------------------------------------------------------------------
# modelisation SPP
# T la matrice (nbVar x nbSC) de coefficients de variables dans les contraintes, soit 0 soit 1
# c le vecteur de coefficients de variables dans la fonction objective

function model_spp(C, A)
    m::Model = Model(GLPK.Optimizer)
    set_optimizer_attribute(m, "msg_lev", GLPK.GLP_MSG_ALL)

    # nombre de variable
    nbVar::Int64 = length(C)
    # nombre de contrainte
    nbSC::Int64 = size(A,1)

    # declaration des variables
    @variable(m, x[1:nbVar], Bin)

    # declaration fonction objective
    @objective(m, Max, sum(C[j]x[j] for j = 1:nbVar))

    # declaration contraintes
    @constraint(m, sc[i = 1:nbSC], sum(A[i,j]x[j] for j in 1:nbVar) <= 1)

    return m
end

# ---------------------------------------------------------------------------
# Construction gloutonne d'une solution admissible de SPP

function GreedyConstruction(C, A)
    x = zeros(Int64, length(C))
    S = 1:length(C)
    M = 1:size(A, 1)
    E = Array{Array{Int64,1},1}(undef, length(C))
    Einit = Array{Array{Int64,1},1}(undef, length(C))
    F = Array{Array{Int64,1},1}(undef, length(C))
    Cout = Array{Float64,1}(undef, length(C))
    firstloop = true
    while (length(S) != 0)
        fill!(E, Int64[])
        fill!(F, Int64[])
        fill!(Cout, -Inf)
        if length(S) == 1
            imax = S[1]
        else
            for j in M
                current = intersect(findall(isequal(1), A[j, :]), S)
                for i in current
                    E[i] = union(E[i], setdiff(current, i))
                    F[i] = union(F[i], j)
                end
            end
            for i in S
                if length(E[i]) != 0
                    Cout[i] = C[i] - sum(C[u] for u in E[i])
                end
            end
            if maximum(Cout) == -Inf        #pas de conflit
                imax = S[1]
            else
                imax = findfirst(isequal(maximum(Cout)), Cout)
            end
        end
        if firstloop
            Einit = copy(E)
            firstloop = false
        end
        x[imax] = 1
        S = setdiff(S, union(E[imax], imax))
        M = setdiff(M, F[imax])
    end
    z = sum(C[i]x[i] for i in 1:length(C))
    return x, z, Einit
end

# ---------------------------------------------------------------------------
# Amelioration gloutonne par recherche locale d’une solution de SPP

function findAllCombination(S, E)
    common = Int64[]
    while S != Int64[]
        if (intersect(S, E[S[1]]) != [])
            S1 = setdiff(S,E[S[1]])
            S2 = setdiff(S,S[1])
            res1 = findAllCombination(S1,E)
            for i in 1:length(res1)
                res1[i] = union(res1[i], common)
            end
            res2 = findAllCombination(S2,E)
            for i in 1:length(res2)
                res2[i] = union(res2[i], common)
            end
            return union(res1,res2)
        else
            common = union(common, S[1])
            S = setdiff(S, S[1])
        end
    end
    return [common]
end

function GreedyImprovement(C, A, xInit, zInit, Einit)
    meilleur = true
    xMax = xInit
    zMax = zInit
    while meilleur
        meilleur = false
        Y = findall(isequal(1), xMax)
        zCurrent = zMax
        xCurrent = copy(xMax)
        for i in Y
            Ek = Int64[]
            for k in setdiff(Y,i)
                Ek = union(Ek,Einit[k])
            end
            S = setdiff(Einit[i],Ek)
            Combis = findAllCombination(S,Einit)
            if Combis != [Int64[]]
                for u in 1:length(Combis)
                    Cv = sum(C[v] for v in Combis[u])
                    Ci = C[i]
                    zPrime = zCurrent - Ci + Cv
                    if zPrime > zMax
                        meilleur = true
                        zMax = zPrime
                        xMax = copy(xCurrent)
                        for v in Combis[u]
                            xMax[v] = 1
                        end
                        xMax[i] = 0
                    end
                end
            end
        end
    end
    return xMax, zMax
end

max_stars_repo_name: ducanh-le/metaSPPstu

max_issues_repo_name: ducanh-le/metaSPPstu

max_forks_repo_name: ducanh-le/metaSPPstu


text:
# Methods related to all things GP regression

mutable struct GPInfo
    gp
    γ_bound::Float64
    RKHS_bound::Float64
    bound_type::String
    post_scale_factor::Float64
    Kinv
end

"""
Select datapoints in a neighborhood of a point.
"""
function get_points_in_neighborhood(center, radius, x_train, y_train)
    norm_array = [norm(row) for row in eachrow(center .- x_train)]
    bool_array = norm_array .<= radius
    x_sub = @view x_train[bool_array, :]
    y_sub = @view y_train[bool_array, :]
    return x_sub, y_sub
end

"""
Create GP From Given Data
"""
function train_gps(x_train, y_train; se_params=[0., 0.65], optimize_hyperparameters=false, lnoise=nothing, opt_fraction=1.0)
    gp_set = Dict()
    dim_keys = ["x$i" for i=1:size(x_train,1)]
    for (i, out_dim) in enumerate(dim_keys)
        # Handle data dependency here
        # x_train_sub = x_train[:, findall(.>(0), data_deps[out_dim])[:]]
        x_train_sub = x_train
        gp = train_gp_1dim(x_train, y_train[i,:]; se_params=se_params, optimize_hyperparameters=optimize_hyperparameters, lnoise=lnoise, opt_fraction=opt_fraction)

        gp_set["x$i"] = deepcopy(gp)
    end
    return gp_set
end

function generate_gp_filename(params; filename_appendix=nothing)
    gps_dir = @sprintf("%s/gps", params.experiment_directory)
    !isdir(gps_dir) && mkpath(gps_dir)
    gps_filename = @sprintf("%s-m%d-σ%1.3f-rs%d-gps", params.system_params.mode_tag, params.data_params.data_num, params.data_params.noise_sigma, params.random_seed)
    gps_filename = isnothing(filename_appendix) ? gps_filename : "$gps_filename-$filename_appendix"
    return "$gps_dir/$gps_filename.bin"
end

function train_gp_1dim(x_train, y_train; se_params=[0., 0.65], optimize_hyperparameters=false, lnoise=nothing, opt_fraction=1.0)
    m_prior = MeanZero()
    k_prior = SE(se_params[2], se_params[1])

    if isnothing(lnoise)
        lnoise = log(sqrt(1+2/size(x_train, 2))) # Generalize to handle any bound
    end

    if optimize_hyperparameters
        @assert 0 < opt_fraction <= 1
        num_opt = Int(opt_fraction*length(x_train))
        opt_idx = StatsBase.sample(1:length(y_train), num_opt, replace = false)
        gp_pre = GP(x_train[:, opt_idx], y_train[opt_idx], m_prior, k_prior, lnoise)
        optimize!(gp_pre)
        gp = update_gp_data(gp_pre, x_train, y_train)
    else
        gp = GP(x_train, y_train, m_prior, k_prior, lnoise)
        if optimize_hyperparameters
            optimize!(gp)
        end
    end

    return gp
end
"""
Generate GP regression given the dataset.
"""
function generate_estimates(params, x_train, y_train)
    data_deps = params.system_params.dependency_dims
    x_train_sub = x_train #[:, findall(.>(0), data_deps[out_dim])[:]]
    optimize_hyperparameters = params.data_params.data_frac_optimize > 0. ? true : false
    gp_set = train_gps(x_train, y_train; se_params=[0., 0.65], optimize_hyperparameters=optimize_hyperparameters, lnoise=nothing, opt_fraction=params.data_params.data_frac_optimize)

    gp_info_dict = Dict()
    for dim_key in keys(data_deps)
        gp_info_dict[dim_key] = create_gp_info(params, gp_set, dim_key)
    end

    return gp_set, gp_info_dict
end

"""
Create Local GP From Global GP
"""
function create_local_gp_info_from_global(x_data, y_data, global_gp_info_dict)

    local_gp_info_dict = copy(global_gp_info_dict)

    for (i, gp_key) in enumerate(keys(global_gp_info_dict))
        gp = global_gp_info_dict[gp_key].gp
        gp_local = update_gp_data(gp, x_data, y_data)
        local_gp_info_dict[gp_key].gp = gp_local
        update_gp_info(local_gp_info_dict[gp_key])
    end

    return local_gp_info_dict
end

"""
Create the GP Info structure componenets.
"""
function create_gp_info(params, gp_set, dim_key)
    gp = gp_set[dim_key]

    scale_factor = params.data_params.bound_type == "rkhs-tight" ? params.data_params.noise_sigma/sqrt(1. + 2. / gp.nobs) : 1.
    domain = params.domain
    diam_domain = 0.
    for dim_key in keys(domain)
        diam_domain += (domain[dim_key][1] - domain[dim_key][2])^2
    end
    diam_domain = sqrt(diam_domain)
    σ_inf = sqrt(gp.kernel.σ2*exp(-1/2*(diam_domain)^2/gp.kernel.ℓ2))

    # Calculating the RKHS parameter bounds
    RKHS_bound = abs(domain[dim_key][2] + params.system_params.lipschitz_bound*diam_domain)/ σ_inf
    B = 1 + (1 + 2/(gp.nobs))^(-1)
    γ = 0.5*gp.nobs*log(B)
    K_inv = inv(gp.cK.mat + exp(gp.logNoise.value)^2*I)
    gp_info = GPInfo(gp, γ, RKHS_bound, params.data_params.bound_type, scale_factor, K_inv)

    return gp_info
end

"""
Predict mean function
"""
function predict_N_steps(gp_info_dict, x0, N)
    xp_trajectory = []
    xc = x0
    dim_keys = keys(gp_info_dict)

    for i=1:N
        xp = zeros(length(dim_keys), 1)
        for (j, key) in enumerate(dim_keys)
            μ, _ = predict_f(gp_info_dict[key].gp, xc)
            xp[j] = μ[1]
        end

        push!(xp_trajectory, xp)
        xc = xp
    end

    return xp_trajectory
end

"""
Predict mean and vairance over one step
"""
function predict_one_step_full(gp_info_dict, x0; known_flag=true)
    μ_dict = Dict()
    σ_dict = Dict()
    dim_keys = keys(gp_info_dict)
    for (i, key) in enumerate(dim_keys)
        μ, σ = predict_f(gp_info_dict[key].gp, x0)
        μ_dict[key] = μ[1] + x0[i]
        σ_dict[key] = σ[1]
    end
return μ_dict, σ_dict
end

"""
Update the GP Info after updating the gp
"""
function update_gp_info(gp_info)
    gp = gp_info.gp
    bound_type = gp_info.bound_type

    old_nobs = size(gp_info.Kinv)[1]
    scale_factor = gp_info.bound_type == "rkhs-tight" ? gp_info.post_scale_factor*sqrt(1. + 2. / (old_nobs))/sqrt(1. + 2. / gp.nobs) : 1.

    B = 1 + (1 + 2/(gp.nobs))^(-1)
    γ = 0.5*gp.nobs*log(B)

    K_inv = inv(gp.cK.mat + exp(gp.logNoise.value)^2*I)

    gp_info.post_scale_factor = scale_factor
    gp_info.Kinv = K_inv
    gp_info.γ_bound = γ

    return nothing
end

"""
Save the GP info and objects.
"""
function save_gp_info(params, gp_set, gp_info_dict; filename_appendix=nothing)
    exp_dir = params.experiment_directory

    gp_save_info = Dict()
    gp_save_info[:gp_set] = gp_set
    gp_save_info[:gp_info] = gp_info_dict

    gps_filename = generate_gp_filename(params, filename_appendix=filename_appendix)

    # Save the GPs for further analysis.
    @info "Saving GP regressions to experiment directory..."
    open(gps_filename, "w") do f
        serialize(f, gp_save_info)
    end
end

function load_gp_info(filename)
    gp_save_info = nothing
    open(filename, "r") do f
        gp_save_info = deserialize(f)
    end
    return gp_save_info[:gp_set], gp_save_info[:gp_info]
end

"""
Update the GP with the specified dataset and return a new GP object
"""
function update_gp_data(gp, x_data, y_data)
    new_gp = deepcopy(gp)
    GaussianProcesses.fit!(new_gp, x_data, y_data)
    return new_gp
end

"""
Update the GP by adding a single datapoint
"""
function add_datapoint_to_gp(gp::GPE, datapoint::Tuple)
    newx = [gp.x datapoint[1]]
    newy = [gp.y; datapoint[2]]
    GaussianProcesses.fit!(gp, newx, newy)
end

"""
Load a GP bin file.
"""
function load_gp(filename)
    f = open(filename)
    gp_set = deserialize(f)
    close(f)
    return gp_set
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Print function of the text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Low-level function to print the table using the text backend.
function _pt_text(
    r_io::Ref{Any}, pinfo::PrintInfo;
    border_crayon::Crayon = Crayon(),
    header_crayon::Union{Crayon, Vector{Crayon}} = Crayon(bold = true),
    subheader_crayon::Union{Crayon,Vector{Crayon}} = Crayon(foreground = :dark_gray),
    rownum_header_crayon::Crayon = Crayon(bold = true),
    text_crayon::Crayon = Crayon(),
    omitted_cell_summary_crayon::Crayon = Crayon(foreground = :cyan),
    alignment_anchor_fallback::Symbol = :l,
    alignment_anchor_fallback_override::Dict{Int, Symbol} = Dict{Int, Symbol}(),
    alignment_anchor_regex::Dict{Int, T} where T <:AbstractVector{Regex} = Dict{Int, Vector{Regex}}(),
    autowrap::Bool = false,
    body_hlines::Vector{Int} = Int[],
    body_hlines_format::Union{Nothing, NTuple{4, Char}} = nothing,
    continuation_row_alignment::Symbol = :c,
    crop::Symbol = get(r_io.x, :limit, false) ? :both : :none,
    crop_subheader::Bool = false,
    crop_num_lines_at_beginning::Int = 0,
    columns_width::Union{Int, AbstractVector{Int}} = 0,
    display_size::Tuple{Int, Int} = displaysize(r_io.x),
    equal_columns_width::Bool = false,
    ellipsis_line_skip::Integer = 0,
    highlighters::Union{Highlighter, Tuple} = (),
    hlines::Union{Nothing, Symbol, AbstractVector} = nothing,
    linebreaks::Bool = false,
    maximum_columns_width::Union{Int, AbstractVector{Int}} = 0,
    minimum_columns_width::Union{Int, AbstractVector{Int}} = 0,
    newline_at_end::Bool = true,
    overwrite::Bool = false,
    noheader::Bool = false,
    nosubheader::Bool = false,
    row_name_crayon::Crayon = Crayon(bold = true),
    row_name_header_crayon::Crayon = Crayon(bold = true),
    row_number_alignment::Symbol = :r,
    show_omitted_cell_summary::Bool = true,
    sortkeys::Bool = false,
    tf::TextFormat = tf_unicode,
    title_autowrap::Bool = false,
    title_crayon::Crayon = Crayon(bold = true),
    title_same_width_as_table::Bool = false,
    vcrop_mode::Symbol = :bottom,
    vlines::Union{Nothing, Symbol, AbstractVector} = nothing
)
    # `r_io` must always be a reference to `IO`. Here, we unpack it. This is
    # done to improve inference and reduce compilation time. Ideally, we need to
    # add the `@nospecialize` annotation to `io`. However, it returns the
    # message that this annotation is not supported with more than 32 arguments.
    io = r_io.x

    # Unpack fields of `pinfo`.
    data                    = pinfo.data
    header                  = pinfo.header
    id_cols                 = pinfo.id_cols
    id_rows                 = pinfo.id_rows
    num_rows                = pinfo.num_rows
    num_cols                = pinfo.num_cols
    num_printed_cols        = pinfo.num_printed_cols
    num_printed_rows        = pinfo.num_printed_rows
    header_num_rows         = pinfo.header_num_rows
    header_num_cols         = pinfo.header_num_cols
    show_row_number         = pinfo.show_row_number
    row_number_column_title = pinfo.row_number_column_title
    show_row_names          = pinfo.show_row_names
    row_names               = pinfo.row_names
    row_name_alignment      = pinfo.row_name_alignment
    row_name_column_title   = pinfo.row_name_column_title
    alignment               = pinfo.alignment
    cell_alignment          = pinfo.cell_alignment
    formatters              = pinfo.formatters
    compact_printing        = pinfo.compact_printing
    title                   = pinfo.title
    title_alignment         = pinfo.title_alignment
    header_alignment        = pinfo.header_alignment
    header_cell_alignment   = pinfo.header_cell_alignment
    cell_first_line_only    = pinfo.cell_first_line_only
    renderer                = pinfo.renderer
    limit_printing          = pinfo.limit_printing

    # Input variables verification and initial setup
    # ==========================================================================

    # Let's create a `IOBuffer` to write everything and then transfer to `io`.
    io_has_color = get(io, :color, false)::Bool
    buf_io       = IOBuffer()
    buf          = IOContext(buf_io, :color => io_has_color)
    display      = Display(has_color = io_has_color)

    # Check which type of cropping the user wants.
    if crop == :both
        display.size = display_size
    elseif crop == :vertical
        display.size = (display_size[1], -1)
    elseif crop == :horizontal
        display.size = (-1, display_size[2])
    else
        # If the table will not be cropped, then we should never show an omitted
        # cell summary.
        show_omitted_cell_summary = false
    end

    # Make sure that `vcrop_mode` is valid.
    vcrop_mode ∉ [:bottom, :middle] && (vcrop_mode = :bottom)

    crop_num_lines_at_beginning < 0 && (crop_num_lines_at_beginning = 0)

    if !noheader && num_cols != header_num_cols
        error("The header length must be equal to the number of columns.")
    end

    # Additional processing necessary if the user wants to print the header.
    if !noheader
        # If the user do not want to print the sub-header but wants to print the
        # header, then just force the number of rows in header to be 1.
        if nosubheader
            # Now, `header` will be a view of the first line of the matrix that
            # has the header.
            header_num_rows = 1
        end

        # Transform some keywords that are single elements to vectors.
        if header_crayon isa Crayon
            header_crayon = fill(header_crayon, num_cols)
        else
            if length(header_crayon) != num_cols
                error("The length of `header_crayon` must be the same as the number of columns.")
            end
        end

        if subheader_crayon isa Crayon
            subheader_crayon = fill(subheader_crayon, num_cols)
        else
            if length(subheader_crayon) != num_cols
                error("The length of `subheader_crayon` must be the same as the number of columns.")
            end
        end
    end

    # Make sure that `highlighters` is always a Ref{Any}(Tuple).
    if !(highlighters isa Tuple)
        highlighters = Ref{Any}((highlighters,))
    else
        highlighters = Ref{Any}(highlighters)
    end

    # Make sure that `maximum_columns_width` is always a vector.
    if maximum_columns_width isa Integer
        maximum_columns_width = fill(maximum_columns_width, num_cols)
    end

    # Make sure that `minimum_columns_width` is always a vector.
    if minimum_columns_width isa Integer
        minimum_columns_width = fill(minimum_columns_width, num_cols)
    end

    # Check which columns must have fixed sizes.
    columns_width isa Integer && (columns_width = fill(columns_width, num_cols))

    if length(columns_width) != num_cols
        error("The length of `columns_width` must be the same as the number of columns.")
    end

    # The number of lines that must be skipped from printing ellipsis must be
    # greater of equal 0.
    (ellipsis_line_skip < 0) && (ellipsis_line_skip = 0)

    # Increase the number of printed columns if the user wants to show
    # additional columns.
    Δc = show_row_names + show_row_number
    num_printed_cols += Δc

    # If the user wants to horizontally crop the printing, then it is not
    # necessary to process all the lines. We will to process, at most, the
    # number of lines in the display.
    if display.size[1] > 0
        num_printed_rows = min(num_printed_rows, display.size[1])
    end

    # If the user wants to vertically crop the printing, then it is not
    # necessary to process all the columns. However, we cannot know at this
    # stage the size of each column. Thus, this initial algorithm uses the fact
    # that each column printed will have at least 4 spaces.
    if display.size[2] > 0
        num_printed_cols = min(num_printed_cols, ceil(Int, display.size[2] / 4))
    end

    # We need to process at least the row number and row name columns to avoid
    # errors.
    num_printed_cols < Δc && (num_printed_cols = Δc)

    # Number of rows and columns after filtering.
    num_filtered_rows = length(id_rows)
    num_filtered_cols = length(id_cols)

    # Create the string matrices that will be printed
    # ==========================================================================

    # Get the string which is printed when `print` is called in each element of
    # the matrix. Notice that we must create only the matrix with the printed
    # rows and columns.
    header_str = Matrix{String}(undef, header_num_rows, num_printed_cols)
    data_str   = Matrix{Vector{String}}(undef, num_printed_rows, num_printed_cols)

    # Auxiliary variables to adjust `alignment` based on the new columns added
    # to the table.
    alignment_add = Symbol[]

    # Vector that must contain the width of each column.
    cols_width = ones(Int, num_printed_cols)

    # Vector that must contain the number of lines in each printed row.
    num_lines_in_row = ones(Int, num_printed_rows)

    # NOTE: Algorithm to compute the table size and number of lines in the rows.
    # Previously, the algorithm to compute the table size and the number of
    # lines in the rows was a function called after the conversion of the input
    # matrix to the string matrix. Although this approach was cleaner, we had
    # problems when computing how many columns we can fit on the display to
    # avoid unnecessary processing. Hence, the functions that fill the data are
    # also responsible to compute the size of each column.

    # Vectors with the order that the rows must be filled
    # --------------------------------------------------------------------------

    # The vector `jvec` contains the indices in the printed matrix whereas the
    # vector `jrvec` contains the indices in the data matrix.
    #
    # NOTE: The vector `jvec` is not sorted when crop mode is `:middle`.
    # This format is required to make sure we always fill one row at the top and
    # one row in the bottom. Otherwise, we can have access to an undefined
    # reference.
    jvec, jrvec = _compute_row_fill_vectors(
        id_rows,
        num_printed_rows,
        vcrop_mode
    )

    # Cell alignment override
    # --------------------------------------------------------------------------

    # Compute the cells in which the alignment is overridden.
    cell_alignment_override = _compute_cell_alignment_override(
        data,
        id_cols,
        id_rows,
        Δc,
        num_printed_cols,
        num_printed_rows,
        # Configurations.
        cell_alignment
    )

    # Row numbers
    # --------------------------------------------------------------------------

    if show_row_number
        _fill_row_number_column!(
            header_str,
            data_str,
            cols_width,
            jvec,
            jrvec,
            noheader,
            row_number_column_title
        )

        # Add information about the row number column to the variables.
        push!(alignment_add, row_number_alignment)
    end

    # Row names
    # --------------------------------------------------------------------------

    if show_row_names
        _fill_row_name_column!(
            header_str,
            data_str,
            cols_width,
            row_names,
            jvec,
            jrvec,
            Δc,
            compact_printing,
            renderer,
            row_name_column_title
        )

        # Add information about the row name column to the variables.
        push!(alignment_add, row_name_alignment)
    end

    !isempty(alignment_add) && (alignment = vcat(alignment_add, alignment))

    # Table data
    # --------------------------------------------------------------------------

    num_printed_cols, num_printed_rows = _fill_matrix_data!(
        header_str,
        data_str,
        cols_width,
        num_lines_in_row,
        id_cols,
        jvec,
        jrvec,
        Δc,
        data,
        header,
        formatters,
        display,
        # Configuration options.
        autowrap,
        cell_first_line_only,
        columns_width,
        compact_printing,
        crop_subheader,
        limit_printing,
        linebreaks,
        maximum_columns_width,
        minimum_columns_width,
        noheader,
        renderer,
        vcrop_mode
    )

    # Column alignment regex
    # --------------------------------------------------------------------------

    _apply_alignment_anchor_regex!(
        data_str,
        cols_width,
        id_cols,
        id_rows,
        Δc,
        # Configurations.
        alignment_anchor_fallback,
        alignment_anchor_fallback_override,
        alignment_anchor_regex,
        cell_alignment_override,
        columns_width,
        maximum_columns_width,
        minimum_columns_width
    )

    # If the user wants all the columns with the same size, then select the
    # larger.
    equal_columns_width && (cols_width = fill(maximum(cols_width), num_printed_cols))

    # Compute where the horizontal and vertical lines must be drawn
    # --------------------------------------------------------------------------

    # Create the format of the horizontal lines.
    if body_hlines_format === nothing
        body_hlines_format = (
            tf.left_intersection,
            tf.middle_intersection,
            tf.right_intersection,
            tf.row
        )
    end

    if hlines === nothing
        hlines = _process_hlines(tf.hlines, body_hlines, num_filtered_rows, noheader)
    else
        hlines = _process_hlines(hlines, body_hlines, num_filtered_rows, noheader)
    end

    # Check if the last horizontal line must be drawn. This is required when
    # computing the moment that the display will be cropped.
    draw_last_hline = (num_filtered_rows + !noheader) ∈ hlines

    # Process `vlines`.
    if vlines === nothing
        vlines = _process_vlines(tf.vlines, num_printed_cols)
    else
        vlines = _process_vlines(vlines, num_printed_cols)
    end

    # Compute the table width
    # --------------------------------------------------------------------------

    table_width = sum(cols_width) + 2length(cols_width) + length(vlines)

    if (display.size[2] > 0) && (table_width > display.size[2])
        table_width = display.size[2]
    end

    # Process the title
    # --------------------------------------------------------------------------

    title_tokens = _tokenize_title(
        title,
        display.size[2],
        table_width,
        # Configurations
        title_alignment,
        title_autowrap,
        title_same_width_as_table
    )

    # Sum the number of lines in the title to the number of lines that must be
    # available at the beginning of the table.
    crop_num_lines_at_beginning += length(title_tokens)

    # Compute the table data printing recipe
    # --------------------------------------------------------------------------

    # Number of additional lines that must be consider to crop the display
    # vertically.
    Δdisplay_lines = 1 + newline_at_end + draw_last_hline + crop_num_lines_at_beginning

    row_printing_recipe, col_printing_recipe, num_omitted_rows, num_omitted_cols =
        _create_printing_recipe(
            display,
            header_num_rows,
            num_filtered_rows,
            num_filtered_cols,
            num_printed_rows,
            num_printed_cols,
            num_lines_in_row,
            cols_width,
            id_rows,
            hlines,
            vlines,
            Δdisplay_lines,
            Δc,
            # Configurations
            crop,
            noheader,
            show_omitted_cell_summary,
            vcrop_mode
        )

    #                           Print the table
    # ==========================================================================

    # Title
    # ==========================================================================

    _print_title(buf, title_tokens, display.has_color, title_crayon)

    # If there is no column or row to be printed, then just exit.
    if (num_printed_cols == 0) || (num_printed_rows == 0)
        @goto print_to_output
    end

    # Top table line
    # ==========================================================================

    if 0 ∈ hlines
        _draw_line!(
            display,
            buf,
            tf.up_left_corner,
            tf.up_intersection,
            tf.up_right_corner,
            tf.row,
            border_crayon,
            cols_width,
            vlines
        )
    end

    # Header
    # ==========================================================================

    # Header and sub-header texts
    # --------------------------------------------------------------------------

    if !noheader
        _print_table_header!(
            buf,
            display,
            header,
            header_str,
            id_cols,
            id_rows,
            num_printed_cols,
            Δc,
            cols_width,
            vlines,
            # Configurations.
            alignment,
            header_alignment,
            header_cell_alignment,
            show_row_names,
            show_row_number,
            tf,
            # Crayons.
            border_crayon,
            header_crayon,
            subheader_crayon,
            rownum_header_crayon,
            row_name_header_crayon
        )

        # Bottom header line
        #-----------------------------------------------------------------------

        if 1 ∈ hlines
            _draw_line!(
                display,
                buf,
                tf.left_intersection,
                tf.middle_intersection,
                tf.right_intersection,
                tf.row,
                border_crayon,
                cols_width,
                vlines
            )
        end
    end

    # Data
    # ==========================================================================

    @inbounds _print_table_data(
        buf,
        display,
        data,
        data_str,
        id_cols,
        id_rows,
        Δc,
        cols_width,
        vlines,
        row_printing_recipe,
        col_printing_recipe,
        # Configurations.
        alignment,
        body_hlines_format,
        cell_alignment_override,
        continuation_row_alignment,
        ellipsis_line_skip,
        highlighters,
        noheader,
        show_row_names,
        show_row_number,
        tf,
        # Crayons.
        border_crayon,
        row_name_crayon,
        text_crayon
    )

    # Bottom table line
    # ==========================================================================

    if draw_last_hline
        _draw_line!(
            display,
            buf,
            tf.bottom_left_corner,
            tf.bottom_intersection,
            tf.bottom_right_corner,
            tf.row,
            border_crayon,
            cols_width,
            vlines
        )
    end

    # Summary of the omitted cells
    # ==========================================================================

    _print_omitted_cell_summary(
        buf,
        display.has_color,
        num_omitted_cols,
        num_omitted_rows,
        omitted_cell_summary_crayon,
        show_omitted_cell_summary,
        table_width
    )

    @label print_to_output

    # Print the buffer
    # ==========================================================================

    _flush_buffer!(io, buf_io, overwrite, newline_at_end, display.row)

    return nothing
end


text:
type LinearInterpolation <: Interpolation
  x_vals::Vector{Float64}
  y_vals::Vector{Float64}
  s::Vector{Float64}
end

# Linear initialize
function initialize!(interp::LinearInterpolation, x_vals::Vector{Float64}, y_vals::Vector{Float64})
  interp.x_vals = x_vals
  interp.y_vals = y_vals
  interp.s = zeros(length(y_vals))

  return interp
end

# Linear Interpolation update
function update!(interp::LinearInterpolation, idx::Int)
  for i = 2:idx
    @inbounds dx = interp.x_vals[i] - interp.x_vals[i - 1]
    @inbounds interp.s[i - 1] = (interp.y_vals[i] - interp.y_vals[i - 1]) / dx
  end

  return interp
end

update!(interp::LinearInterpolation) = update!(interp, length(interp.y_vals))

function value(interp::LinearInterpolation, val::Float64)
  i = locate(interp, val)
  # println("I is ", i)
  # println("Val is ", val)
  # println("Y vals ", interp.y_vals)
  # println("X vals ", interp.x_vals)
  # println("S vals ", interp.s)
  return interp.y_vals[i] + (val - interp.x_vals[i]) * interp.s[i]
end

function derivative(interp::LinearInterpolation, val::Float64)
  i = locate(interp, val)
  return interp.s[i]
end


text:
"""
    LinearView(offset, stride)

Subtype of `IndexTransform` that provides linear indexing for `Base.FastSubArray` and
`FastContiguousSubArray`.
"""
struct LinearView{O<:CanonicalInt,S<:CanonicalInt} <: IndexTransform{1,1}
    offset::O
    stride::S
end

const Offset{O} = LinearView{O,StaticInt{1}}
Offset(o::CanonicalInt) = LinearView(o, static(1))

(t::LinearView)(x::CanonicalInt) = getfield(t, :offset) + getfield(t, :stride) * x
@inline function (t::LinearView)(x::AbstractRange{Int})
    t(static_first(i)):t(static_step(i)):t(static_last(i))
end
@inline function (t::LinearView)(x::AbstractArray{Int})
    TransformedIndex{Int,ndims(i),typeof(i),typeof(t)}(i, t)
end
@inline function (t1::LinearView)(t2::LinearView)
    LinearView(
        getfield(t2, :offset) + getfield(t1, :offset) * getfield(t2, :stride),
        getfield(t1, :stride) * getfield(t2, :stride)
    )
end

function Base.show(io::IO, ::MIME"text/plain", @nospecialize(x::LinearView))
    if static(1) === x.stride
        print(io, "Offset($(x.offset))")
    else
        print(io, "LinearView(offset=$(x.offset), stride=$(x.stride))")
    end
end


text:
using DelayEmbeddings, ChaosTools, OrdinaryDiffEq
using Test

println("\nTesting R-param. estimation...")

"""
    saturation_point(x, y; threshold = 0.01, dxi::Int = 1, tol = 0.2)
Decompose the curve `y(x)` into linear regions using `linear_regions(x, y; dxi, tol)`
and then attempt to find a saturation point where the the first slope
of the linear regions become `< threshold`.

Return the `x` value of the saturation point.
"""
function saturation_point(Ds, E1s; threshold = 0.01, kwargs...)
    lrs, slops = ChaosTools.linear_regions(Ds, E1s; kwargs...)
    i = findfirst(x -> x < threshold, slops)
    return i == 0 ? Ds[end] : Ds[lrs[i]]
end

test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax

@testset "Estimate Dimension" begin
    s_sin = sin.(0:0.1:1000)
    τ = 15
    Ds = 1:5
    E1s = DelayEmbeddings.estimate_dimension(s_sin, τ, Ds)
    E2s = DelayEmbeddings.stochastic_indicator(s_sin, τ, Ds)
    @test saturation_point(Ds,E1s; threshold=0.01) == 1

    # @test minimum(E2s) < 0.1 # THIS TEST FAILS

    ds = Systems.roessler();τ=15; dt=0.1
    data = trajectory(ds,1000;dt=dt,diffeq...)
    s_roessler = data[:,1]
    Ds = 1:5
    E1s = DelayEmbeddings.estimate_dimension(s_roessler, τ, Ds)
    @test saturation_point(Ds,E1s; threshold=0.1) ∈ [2, 3]
    E2s = DelayEmbeddings.stochastic_indicator(s_roessler, τ, Ds)
    @test minimum(E2s) < 0.3

    E1s = DelayEmbeddings.estimate_dimension(s_roessler, τ, Ds;metric = Chebyshev())
    @test saturation_point(Ds,E1s; threshold=0.1) ∈ [2, 3]

    ds = Systems.lorenz();τ=5; dt=0.01
    data = trajectory(ds,500;dt=dt,diffeq...)
    s_lorenz = data[:,1]
    Ds = 1:5
    E1s = DelayEmbeddings.estimate_dimension(s_lorenz, τ, Ds)
    E2s = DelayEmbeddings.stochastic_indicator(s_lorenz, τ, Ds)
    @test saturation_point(Ds,E1s; threshold=0.1) ∈ [2, 3]

    # @test minimum(E2s) < 0.1 # THIS TEST FAILS

    #Test against random signal
    E2s = DelayEmbeddings.stochastic_indicator(rand(100000), 1, 1:5)
    @test minimum(E2s) > 0.9

    # Test `fnn` method

    τ = 15
    Ds = 1:5
    number_fnn = DelayEmbeddings.estimate_dimension(s_sin, τ, Ds, "fnn")
    @test findfirst(number_fnn .≈ 0.0) == 1

    τ = 15
    Ds = 1:5
    number_fnn = DelayEmbeddings.estimate_dimension(s_roessler, τ, Ds, "fnn"; rtol=15)
    @test findfirst(number_fnn .≈ 0.0) ∈ [2, 3]

    τ = 1
    Ds = 1:5
    number_fnn = DelayEmbeddings.estimate_dimension(s_lorenz, τ, Ds, "fnn"; atol=1, rtol=3.0)
    @test findfirst(number_fnn .≈ 0.0) ∈ [2, 3]

    # Test `f1nn` method

    τ = 15
    Ds = 1:5
    ffnn_ratio = DelayEmbeddings.estimate_dimension(s_sin, τ, Ds, "f1nn")
    @test saturation_point(Ds,(1 .- ffnn_ratio); threshold=0.1) == 1

    τ = 15
    Ds = 1:5
    ffnn_ratio = DelayEmbeddings.estimate_dimension(s_roessler, τ, Ds, "f1nn")
    @test saturation_point(Ds,(1 .- ffnn_ratio); threshold=0.1) ∈ [2, 3]
    ffnn_ratio = DelayEmbeddings.estimate_dimension(s_roessler, τ, Ds, "f1nn";metric=Chebyshev())
    @test saturation_point(Ds,(1 .- ffnn_ratio); threshold=0.1) ∈ [2, 3]

    τ = 15
    Ds = 1:5
    ffnn_ratio = DelayEmbeddings.estimate_dimension(s_lorenz, τ, Ds, "f1nn")
    @test saturation_point(Ds,(1 .- ffnn_ratio); threshold=0.1) ∈ [2, 3]

end


text:
using BenchmarkTools: @btime

function memory(start, target)
  numbers = zeros(Int32, target)
  for (i,n) in enumerate(start[1:end-1])
    numbers[n+1]=i
  end
  turn = length(start)
  number = start[end]
  while turn < target
    seen = numbers[number+1]
    numbers[number+1] = turn
    number = seen == 0 ? 0 : turn - seen
    turn += 1
  end
  number
end

input = [0,5,4,1,10,14,7]
part1 = memory(input, 2020)
println("$part1")

part2 = @btime memory(input, 30000000)
println("$part2")


text:
module ParamEstim

const MODEL_PATH = "models/Nakakuki_Cell_2010_ODE"

export
    C,
    V,
    observables,
    observables_index,
    Sim,
    Exp,
    Viz,
    param_values,
    initial_values,
    objective,
    get_search_index,
    get_search_region,
    decode_gene2val,
    encode_val2gene,
    encode_bestIndivVal2randGene,
    update_param,
    simulate_all,
    ga_v1,
    ga_v1_continue,
    ga_v2,
    ga_v2_continue,
    MODEL_PATH

using Printf
using LinearAlgebra
using Random
using StatsBase
using Statistics
using DelimitedFiles
using PyPlot

import Seaborn

include(MODEL_PATH * "/name2idx/parameters.jl")
include(MODEL_PATH * "/name2idx/species.jl")
include(MODEL_PATH * "/set_model.jl")
include(MODEL_PATH * "/observable.jl")
include(MODEL_PATH * "/experimental_data.jl")
include(MODEL_PATH * "/simulation.jl")
include(MODEL_PATH * "/viz.jl")
include(MODEL_PATH * "/fitness.jl")
include(MODEL_PATH * "/set_search_param.jl")

include("ga/initial_population.jl")
include("ga/undxmgg.jl")
include("ga/converging.jl")
include("ga/local_search.jl")
include("ga/v1.jl")
include("ga/v2.jl")

include("dynamics/signaling_systems.jl")
include("dynamics/temporal_dynamics.jl")

end  # module

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim

max_forks_repo_name: himoto/paramEstim