text:
# function combinevectors(x::Vector{Vector{T}})::Vector{T} where T

#     if isempty(x)
#         return Vector{T}(undef, 0)
#     end

#     N = sum(length(x[i]) for i = 1:length(x))

#     y = Vector{T}(undef,N)

#     st_ind = 0
#     fin_ind = 0
#     for i = 1:length(x)
#         st_ind = fin_ind + 1
#         fin_ind = st_ind + length(x[i]) - 1

#         y[st_ind:fin_ind] = x[i]
#     end

#     return y
# end

function array2matrix(X::Vector{Vector{T}})::Matrix{T} where T

    N = length(X)
    D = length(X[1])

    out = Matrix{T}(undef,D,N)
    for n = 1:N
        out[:,n] = X[n]
    end

    return out
end

###########  density fit-related

# Boyd's convex optimization convention: log(-2) := Inf.
function evalcvxlog(x::T)::T where T <: Real
    #println("x = ", x)
    if x < zero(T)
        return Inf
    end

    return log(x)
end

function RKHSfitdensitycostfunc(α::Vector{T}, Kp::Matrix{T}, y::Vector{T}, μ::T)::T where T
    N = length(y)
    @assert length(α) == size(Kp,1) == size(Kp,2)

    r = Kp*α-y
    term1 = dot(r,r)
    term2 = μ*dot(α, Kp, α)
    obj = term1 + term2

    return obj
end
# # test code.
# θ = θ_canonical
# y = vec(𝑌)
# μ = σ²
#
# K = constructkernelmatrix(𝑋, θ)
# #α_solution = fitpdfviaSDP(y, K, μ, max_iters)
#
# # verify objective.
# obj = mycostfunc(α_SDP, K, y, μ)
# # this objective scaoe is the same as the cone solver's. 1.69e-01 answer.

########### end of density fit-related

function findboundingbox(X::Vector{Vector{T}}) where T
    D = length(X[1])
    N = length(X)

    limit_a = Vector{T}(undef, D)
    limit_b = Vector{T}(undef, D)

    for d = 1:D
        Xd = collect( X[n][d] for n = 1:N )
        sorted_Xd = sort(Xd)

        limit_a[d] = sorted_Xd[1]
        limit_b[d] = sorted_Xd[end]
    end

    return limit_a, limit_b
end

function bandpasstotalvariations(LP_percentages,
                        LP_attenuations,
                        HP_percentages,
                        HP_attenuations,
                        y::Vector{T}) where T
    #
    N = length(LP_percentages)
    @assert N == length(LP_attenuations) ==
                length(HP_percentages) ==
                length(HP_attenuations)

    #
    tv_array = Vector{T}(undef, N)

    for i = 1:N
        y_bp = SignalTools.applybandpass( LP_percentages[i],
                                LP_attenuations[i],
                                HP_percentages[i],
                                HP_attenuations[i],
                                y )
        #
        tv_array[i] = sum( abs.(y_bp) )
    end

    return tv_array
end

function bandpassmaxmagnitudes(LP_percentages,
                        LP_attenuations,
                        HP_percentages,
                        HP_attenuations,
                        y::Vector{T}) where T
    #
    N = length(LP_percentages)
    @assert N == length(LP_attenuations) ==
                length(HP_percentages) ==
                length(HP_attenuations)

    #
    max_magnitude = Vector{T}(undef, N)

    for i = 1:N
        y_bp = SignalTools.applybandpass( LP_percentages[i],
                                LP_attenuations[i],
                                HP_percentages[i],
                                HP_attenuations[i],
                                y )
        #
        max_magnitude[i] = maximum( abs.(fft(y_bp)) )
    end

    return max_magnitude
end

# 1. FGPGM for parameter identification in systems of nonlinear ODEs,
#       supplementary.
# Let V be the derivative variables of X.
# D_mat is equation 34 of [1]
function evalderivative(x::Vector{T},
                        μ_GP::Vector{T},
                        μ_GP_dot::Vector{T},
                        θ::KT) where {T <: Real, KT}

    #
    K_XX, K_XV, K_VV = constructkernelematrices(X, θ)

    #K_VX = copy(K_XV')
    μ_dot, S_dot = Utilities.getMVNmarginalparams(x, μ_GP_dot, μ_GP,
                                    K_VV, K_XV, K_XX)
    #
    D_mat = K_XV*inv(K_XX)

    return μ_dot, S_dot, D_mat
end

# b splines

# constant (flat) boundary condition.
# see https://github.com/JuliaMath/Interpolations.jl/blob/master/src/b-splines/cubic.jl
function cubicitp1D(c, x::T, a::T, b::T) where T <: Real
    if one(T) <= x <= length(c)

        δx = abs(x - floor(x))
        i = round(Int, x-δx)

        # out = 1/6*c[i-1]*(1+i-δx)^3 + c[i]*(2/3 - (δx-i)^2 + 1/2*(δx-i)^3) +
        #          c[i+1]*(2/3 - (1+i-δx)^2 + 1/2*(1+i-δx)^3) + 1/6*c[i+2]*(δx-i)^3

        out = 1/6*c[i-1]*(1+i-x)^3 + c[i]*(2/3 - (x-i)^2 + 1/2*(x-i)^3) +
                 c[i+1]*(2/3 - (1+i-x)^2 + 1/2*(1+i-x)^3) + 1/6*c[i+2]*(x-i)^3

        return out
    end

    if x < one(T)
        return a
    end

    return b
end

# # ignore boundary conditions for now.
function cubicitp2D(c, x::Vector{T}, N::Vector{Int}, a, b) where T <: Real

    @assert one(T) <= x[1] <= convert(T,N[1])
    @assert one(T) <= x[2] <= convert(T,N[2])

    δx1 = abs(x[1] - floor(x[1]))
    i = round(Int, x[1]-δx1)

    δx2 = abs(x[2] - floor(x[2]))
    j = round(Int, x[2]-δx2)

    running_sum = zero(T)
    for 𝑖 in Iterators.product(i-1:i+2, j-1:j+2)
        𝑖_vec = collect(𝑖)
        running_sum += c[𝑖_vec...]*evalBspline3(x-𝑖_vec)
    end

    return running_sum
end

function evalBspline3(x::Vector{T}) where T <: Real
    return prod( evalBspline3(x[d]) for d = 1:length(x) )
end

function evalBspline3(x::T) where T <: Real
    abs_x = abs(x)

    if zero(T) <= abs_x < one(T)
        return 2/3 - abs_x^2 + 0.5*abs_x^3
    end

    if one(T) <= abs_x < 2.0
        return 1/6 * (2-abs_x)^3
    end

    return zero(T)
end

function cubicitp2Dexplicit(c, x::Vector{T}, N::Vector{Int}) where T <: Real

    @assert one(T) <= x[1] <= convert(T,N[1])
    @assert one(T) <= x[2] <= convert(T,N[2])

    # determine anchor, 𝑖 = (i1,i2)
    δx1 = abs(x[1] - floor(x[1]))
    i1 = round(Int, x[1]-δx1)

    δx2 = abs(x[2] - floor(x[2]))
    i2 = round(Int, x[2]-δx2)

    x1_lag = x[1]-i1
    x2_lag = x[2]-i2

    x1_adv = 1 - x1_lag
    x2_adv = 1 - x2_lag

    out  = c[i1-1, i2-1]*evalBspline3far(x1_lag)*evalBspline3far(x2_lag)
    out += c[i1, i2-1]*evalBspline3near(x1_lag)*evalBspline3far(x2_lag)
    out += c[i1-1, i2]*evalBspline3far(x1_lag)*evalBspline3near(x2_lag)
    out += c[i1, i2]*evalBspline3near(x1_lag)*evalBspline3near(x2_lag)

    out += c[i1-1, i2+1]*evalBspline3far(x1_lag)*evalBspline3near(x2_adv)
    out += c[i1, i2+1]*evalBspline3near(x1_lag)*evalBspline3near(x2_adv)
    out += c[i1-1, i2+2]*evalBspline3far(x1_lag)*evalBspline3far(x2_adv)
    out += c[i1, i2+2]*evalBspline3near(x1_lag)*evalBspline3far(x2_adv)

    out += c[i1+1, i2-1]*evalBspline3near(x1_adv)*evalBspline3far(x2_lag)
    out += c[i1+2, i2-1]*evalBspline3far(x1_adv)*evalBspline3far(x2_lag)
    out += c[i1+1, i2]*evalBspline3near(x1_adv)*evalBspline3near(x2_lag)
    out += c[i1+2, i2]*evalBspline3far(x1_adv)*evalBspline3near(x2_lag)

    out += c[i1+1, i2+1]*evalBspline3near(x1_adv)*evalBspline3near(x2_adv)
    out += c[i1+2, i2+1]*evalBspline3far(x1_adv)*evalBspline3near(x2_adv)
    out += c[i1+1, i2+2]*evalBspline3near(x1_adv)*evalBspline3far(x2_adv)
    out += c[i1+2, i2+2]*evalBspline3far(x1_adv)*evalBspline3far(x2_adv)

    return out
end

function evalBspline3far(x::T) where T <: Real

    return 1/6 * (1-x)^3 # note this isn't (2-abs(x)), since we're using xk_adv.
    # taken from https://github.com/JuliaMath/Interpolations.jl/blob/master/src/b-splines/cubic.jl
end

function evalBspline3near(x::T) where T <: Real
    abs_x = abs(x)

    return 2/3 - x^2 + 0.5*x^3
end

# with respect to the second coordinate.
function cubicitp2Dexplicit2(c, x::Vector{T}, N::Vector{Int}) where T <: Real

    @assert one(T) <= x[1] <= convert(T,N[1])
    @assert one(T) <= x[2] <= convert(T,N[2])

    # determine anchor, 𝑖 = (i1,i2)
    δx1 = abs(x[1] - floor(x[1]))
    i1 = round(Int, x[1]-δx1)

    δx2 = abs(x[2] - floor(x[2]))
    i2 = round(Int, x[2]-δx2)

    x1_lag = x[1]-i1
    x2_lag = x[2]-i2

    x1_adv = 1 - x1_lag
    x2_adv = 1 - x2_lag

    a3 = zero(T)
    a2 = zero(T)
    a1 = zero(T)
    a0 = zero(T)

    k = i2

    ### x2_lag.

    # c[i1-1, i2-1]*evalBspline3far(x1_lag)*evalBspline3far(x2_lag)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3farlag(a3, a2, a1, a0,
                        c[i1-1, i2-1]*evalBspline3far(x1_lag), k)

    # c[i1, i2-1]*evalBspline3near(x1_lag)*evalBspline3far(x2_lag)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3farlag(a3, a2, a1, a0,
                        c[i1, i2-1]*evalBspline3near(x1_lag), k)

    # c[i1-1, i2]*evalBspline3far(x1_lag)*evalBspline3near(x2_lag)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3nearlag(a3, a2, a1, a0,
                        c[i1-1, i2]*evalBspline3far(x1_lag), k)

    # c[i1, i2]*evalBspline3near(x1_lag)*evalBspline3near(x2_lag)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3nearlag(a3, a2, a1, a0,
                        c[i1, i2]*evalBspline3near(x1_lag), k)

    ##
    # c[i1-1, i2+1]*evalBspline3far(x1_lag)*evalBspline3near(x2_adv)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3nearadv(a3, a2, a1, a0,
                        c[i1-1, i2+1]*evalBspline3far(x1_lag), k)

    # c[i1, i2+1]*evalBspline3near(x1_lag)*evalBspline3near(x2_adv)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3nearadv(a3, a2, a1, a0,
                        c[i1, i2+1]*evalBspline3near(x1_lag), k)

    # c[i1-1, i2+2]*evalBspline3far(x1_lag)*evalBspline3far(x2_adv)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3faradv(a3, a2, a1, a0,
                        c[i1-1, i2+2]*evalBspline3far(x1_lag), k)

    # c[i1, i2+2]*evalBspline3near(x1_lag)*evalBspline3far(x2_adv)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3faradv(a3, a2, a1, a0,
                        c[i1, i2+2]*evalBspline3near(x1_lag), k)

    ##
    # c[i1+1, i2-1]*evalBspline3near(x1_adv)*evalBspline3far(x2_lag)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3farlag(a3, a2, a1, a0,
                        c[i1+1, i2-1]*evalBspline3near(x1_adv), k)

    # c[i1+2, i2-1]*evalBspline3far(x1_adv)*evalBspline3far(x2_lag)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3farlag(a3, a2, a1, a0,
                        c[i1+2, i2-1]*evalBspline3far(x1_adv), k)

    # c[i1+1, i2]*evalBspline3near(x1_adv)*evalBspline3near(x2_lag)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3nearlag(a3, a2, a1, a0,
                        c[i1+1, i2]*evalBspline3near(x1_adv), k)

    # c[i1+2, i2]*evalBspline3far(x1_adv)*evalBspline3near(x2_lag)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3nearlag(a3, a2, a1, a0,
                        c[i1+2, i2]*evalBspline3far(x1_adv), k)

    ##
    # c[i1+1, i2+1]*evalBspline3near(x1_adv)*evalBspline3near(x2_adv)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3nearadv(a3, a2, a1, a0,
                        c[i1+1, i2+1]*evalBspline3near(x1_adv), k)

    # c[i1+2, i2+1]*evalBspline3far(x1_adv)*evalBspline3near(x2_adv)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3nearadv(a3, a2, a1, a0,
                        c[i1+2, i2+1]*evalBspline3far(x1_adv), k)

    # c[i1+1, i2+2]*evalBspline3near(x1_adv)*evalBspline3far(x2_adv)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3faradv(a3, a2, a1, a0,
                        c[i1+1, i2+2]*evalBspline3near(x1_adv), k)

    # c[i1+2, i2+2]*evalBspline3far(x1_adv)*evalBspline3far(x2_adv)
    a3, a2, a1, a0 = getmonomialcoefficientsBspline3faradv(a3, a2, a1, a0,
                        c[i1+2, i2+2]*evalBspline3far(x1_adv), k)

    return a3, a2, a1, a0
end

# monomial coefficients for evalBspline3far(adv)
## for β_far(x_k,adv): expand ((x-k)^3)/6
# x^3/6 -k*x^2/2 +k^2*x/2 -k^3/6
function getmonomialcoefficientsBspline3faradv(a3, a2, a1, a0,
                                                B::T,
                                                k::Int)::Tuple{T,T,T,T} where T <: Real
    b3 = 1/6
    b2 = -k/2
    b1 = k^2/2
    b0 = -k^3/6

    return B*b3+a3, B*b2+a2, B*b1+a1, B*b0+a0
end

# monomial coefficients for evalBspline3far(lag)
## for β_far(x_k,lag): expand ((1-(x-k))^3)/6
# k^3/6 - k^2*x/2 + k^2/2 + k*x^2/2 - k*x + k/2 - x^3/6 + x^2/2 - x/2 + 1/6
function getmonomialcoefficientsBspline3farlag(a3, a2, a1, a0,
                                                B::T,
                                                k::Int)::Tuple{T,T,T,T} where T <: Real
    b3 = -1/6
    b2 = (1 + k)/2
    b1 = -(1 + 2*k + k^2)/2
    b0 = (1 + 3*k + 3*k^2 + k^3)/6

    return B*b3+a3, B*b2+a2, B*b1+a1, B*b0+a0
end

# monomial coefficients for evalBspline3near(adv)
## for β_near(x_k,adv): expand 2/3 -(1-(x-k))^2+((1-(x-k))^3)/2
# 2/3 + k^3/2 - 3*k^2*x/2 + k^2/2 + 3*k*x^2/2 -k*x -k/2 -x^3/2 +x^2/2 +x/2 -0.5
function getmonomialcoefficientsBspline3nearadv(a3, a2, a1, a0,
                                                B::T,
                                                k::Int)::Tuple{T,T,T,T} where T <: Real
    b3 = -0.5
    b2 = (1 + 3*k )/2
    b1 = (1 - 2*k - 3*k^2)/2
    b0 = (-1 -k + k^2 + k^3)/2 + 2/3

    return B*b3+a3, B*b2+a2, B*b1+a1, B*b0+a0
end

# monomial coefficients for evalBspline3near(lag)
## for β_near(x_k,lag): expand 2/3 -(x-k)^2 +((x-k)^3)/2
# -k^3/2 + 3*k^2*x/2 -k^2 -3*k*x^2/2 +2*k*x +x^3/2 -x^2 +2/3
function getmonomialcoefficientsBspline3nearlag(a3, a2, a1, a0,
                                                B::T,
                                                k::Int)::Tuple{T,T,T,T} where T <: Real
    b3 = 0.5
    b2 = -1 - 1.5*k
    b1 = 2*k + 1.5*k^2
    b0 = 2/3 - k^2 - k^3/2

    return B*b3+a3, B*b2+a2, B*b1+a1, B*b0+a0
end

### drawing RVs.
function drawcategorical(w::Vector{T})::Int where T <: Real
    x = rand()

    tmp = 0.0
    for k = 1:length(w)-1
        if x < w[k]+tmp
            return k
        end
        tmp += w[k]
    end

    return length(w)
end

### iterated Brownian Bridge

function Bernoullipolynomial0(x::T)::T where T <: Real
    return one(T)
end

function Bernoullipolynomial1(x::T)::T where T <: Real
    return x-0.5
end

function Bernoullipolynomial2(x::T)::T where T <: Real
    return x^2 -x +1/6
end

function Bernoullipolynomial3(x::T)::T where T <: Real
    return x^3 -1.5*x^2 + 0.5*x
end

function Bernoullipolynomial4(x::T)::T where T <: Real
    return x^4 -2*x^3 +x^2 -1/30
end

function Bernoullipolynomial5(x::T)::T where T <: Real
    return x^5 -2.5*x^4 +5/3*x^3 -1/6*x
end

function Bernoullipolynomial6(x::T)::T where T <: Real
    return x^6 -3*x^5 +5/2*x^4 -0.5*x^2 +1/42
end

### visualization helpers.

# options for warp map samples.
# this is the case when the target density is unknown, but realizations are available.
function getwarpmapsamplecustom( y::Array{T,D},
                            ω_set,
                            pass_band_factor) where {T,D}
    #
    N_bands = length(ω_set)

    Y = y

    #### Split-band analysis.
    ϕY, ψY = SignalTools.runsplitbandanalysis(Y, ω_set, SignalTools.getGaussianfilters)
    ηY = SignalTools.runbandpassanalysis(Y, ω_set, pass_band_factor, SignalTools.getGaussianfilters)

    # #### Riesz transform on the different filtered signals.
    # H, ordering = gethigherorderRTfilters(Y,order)
    #
    # 𝓡ϕY = collect( RieszAnalysisLimited(ϕY[s],H) for s = 1:N_bands)
    # 𝓡ψY = collect( RieszAnalysisLimited(ψY[s],H) for s = 1:N_bands)
    # 𝓡ηY = collect( RieszAnalysisLimited(ηY[s],H) for s = 1:N_bands)

    ϕ_set = ηY

    ϕ = reduce(+,ϕ_set)./N_bands

    return ϕY, ψY, ηY
end

function getwarpmaplinear(ϕ::Array{T,D}) where {T,D}

    itp_ϕ = Interpolations.interpolate(ϕ,
                Interpolations.BSpline(Interpolations.Linear()))

    etp_ϕ = Interpolations.extrapolate(itp_ϕ, Interpolations.Line())

    return etp_ϕ
end

# 1-D CDF: integrate from -Inf to b.
# f is the probability density function.
# check the Jacobian: derivative of b-x/(1-x)-0.5.
# check the b bound becomes -1: evaluate b-x/(1-x)-0.5, x = -1
# for use with HCubature, so vector-valued input.
function get1DCDFintegrand(b::T, f::Function) where T <: Real

    return tt->f(b - tt[1]/(1-tt[1]) -0.5)*(1/(1-tt[1])^2)
end

function getKRCDFintegrand(b::T, f::Function) where T <: Real

    return tt->f(b - tt[end]/(1-tt[end]) -0.5)*(1/(1-tt[end])^2) # Sign flip because we have the bounds reversed.
end

# compute the normalizing constant.
# for use with HCubature, so vector-valued input.
function get1DZintegrand(f::Function) where T <: Real

    return tt->f(tt[1]/(1-tt[1]^2))*((1+tt[1]^2)/(1-tt[1]^2)^2)
end

# computes 1D normalizing constant in the last input component of f.
function getKRZintegrand(f::Function) where T <: Real

    return tt->f(tt[end]/(1+tt[end])^2)*(1+tt[end]^2)/(1-tt[end]^2)^2
end

# computes 1D normalizing constant in the last input component of f.
function getKRZintegrandlogspace(f::Function) where T <: Real

    return tt->exp( log(f(tt[end]/(1+tt[end])^2)) + log(1+tt[end]^2) - 2*log(1-tt[end]^2) )
end

# returns positions, weights. Not working.
function gettanhsinhquad(N::Int, m::Int)::Tuple{Vector{Float64},Vector{Float64}}
    half_π = π/2

    j_range = -N:N
    h = 1/2^m

    g = tt->tanh(half_π*sinh(tt))
    dg = tt->half_π *cosh(tt)*sech(half_π*sinh(tt))^2 # derivative of tanh(pi/2*sinh(t))

    x = collect( g(h*j) for j in j_range )
    w = collect( dg(h*j) for j in j_range )

    return x, w
end

# evaluates a tanh-sinh quadature integral.
function computetsq(x, w, f, m)
    N_nodes = length(w)
    h = 1/2^m

    return h*sum( w[i]*f(x[i]) for i = 1:N_nodes )
    #return sum( w[i]*f(x[i]) for i = 1:N_nodes )
end

# coordinates.
# x ∈ prod( [a[d], b[d]] for d = 1:D ).
# same sampling interval for each dimension.
# M[d] is the number of samples in dimension d.
# out ∈ prod( 1:M[d] for d = 1:D )
function convert2itpindex(  x::Vector{T},
                            a::Vector{T},
                            b::Vector{T},
                            M::Vector{Int})::Vector{T} where T <: Real

    D = length(x)
    @assert D == length(a) == length(b)

    out = Vector{T}(undef,D)
    for d = 1:D
        offset = b[d]-a[d]
        u = (x[d]-a[d])/offset # ∈ [0,1]
        i = u*(M[d]-1) # ∈ [0,M-1]
        out[d] = i + 1
    end

    return out
end
# x =[-0.01; 2]
# a = [-0.1; -1 ]
# b = [0.0; 3]
# M = [10; 5]
# p = convert2itpindex([-0.01; 3], a, b, M) # should be [9;5]

# marginalize over ℝ for the dimensions specified in u_dims.
# Uses the change-of-variable here:
# https://github.com/stevengj/cubature/blob/master/README.md#infinite-intervals
# function marginalize(f::Function, u_dims::Vector{Float64}, ϵ = 1e-4)
#
# (1+tt[d]^2)/(1-tt[d]^2)^2)
#
#     d = u[1]
#     integrand_func = tt->f(tt[d]/(1-tt[d]^2)
#     for d in u[1:end]
#         integrand_func = tt->integrand_func(tt[d]/(1-tt[d]^2)
#     end
#
#     x_min = -one(Float64)
#     x_max = one(Float64)
#     (val,err) = HCubature.hcubature(integrand_func, [x_min], [x_max];
#                     norm = LinearAlgebra.norm, rtol = sqrt(eps(Float64)), atol = 0,
#                     maxevals = typemax(Int), initdiv=1)
#
#     return val, err
# end

function iterativesin()
    N = 3
    f_array = Vector{Function}(undef,N)

    f_array[1] = xx->sin(xx)
    for i = 2:N
        f_array[i] = xx->f_array[i-1](sin(xx))
    end

    return f_array[end]
end

# # f is the joint pdf.
# # marginalizes over variables N:D, where D is the total number of variables.
# # assumes the N variables have support ℝ.
# function getintegrandformarginalization(v::Vector{T},
#                                         f::Function,
#                                         N::Int,
#                                         D::Int) where T <: Real
#     @assert 0 < N < D # at least one variable remains.
#     @assert length(v) == D-N
#
#     g = tt->tt/(1-tt^2)
#     u = xx->collect( g(xx[d]) for d = N:D )
#
#     J = tt->( (1+tt^2)/(1-tt^2)^2 )
#
#     integrand_func = xx->f( [v; u(xx)] )*prod( J(xx[d]) for d = N:D )
#
#     return integrand_func
# end
#
# # f is joint pdf.
# # based on numerical integration.
# function evalmarginalpdf(f::Function, v::Vector{Float64}, N::Int, D::Int, max_integral_evals::Int = 1000)
#
#     integrand_func = getintegrandformarginalization(v, f, N, D)
#
#     x_min = collect( -1.0 for n = N:D )
#     x_max = collect( 1.0 for n = N:D )
#     (val,err) = HCubature.hcubature(integrand_func, x_min, x_max;
#                     norm = LinearAlgebra.norm, rtol = sqrt(eps(Float64)), atol = 0,
#                     maxevals = max_integral_evals, initdiv = 1)
#     #
#     return val
# end

# f is the joint pdf.
# marginalizes over variables 1:N. D is the total number of variables.
# assumes the N variables have support ℝ.
function getintegrandformarginalization(v::Vector{T},
                                        f::Function,
                                        N::Int,
                                        D::Int) where T <: Real
    @assert N < D # at least one variable remains.
    @assert length(v) == D-N

    g = tt->tt/(1-tt^2)
    u = xx->collect( g(xx[d]) for d = 1:N )

    J = tt->( (1+tt^2)/(1-tt^2)^2 )

    integrand_func = xx->f( [u(xx); v] )*prod( J(xx[d]) for d = 1:N )
    #integrand_func = xx->f( [u(xx); u(v)] )*prod( J(xx[d]) for d = 1:length(xx) )*prod( J(v[d]) for d = 1:length(v) )

    return integrand_func
end
# f is joint pdf.
# based on numerical integration.
function evalmarginalpdf(f::Function, v::Vector{Float64}, N::Int, D::Int, max_integral_evals::Int = 1000)

    integrand_func = getintegrandformarginalization(v, f, N, D)

    x_min = collect( -1.0 for n = 1:N )
    x_max = collect( 1.0 for n = 1:N )
    (val,err) = HCubature.hcubature(integrand_func, x_min, x_max;
                    norm = LinearAlgebra.norm, rtol = sqrt(eps(Float64)), atol = 0,
                    maxevals = max_integral_evals, initdiv = 1)
    #
    return val
end

# f is the joint pdf.
# marginalizes over the last N variables, where D is the total number of variables.
# assumes the N variables have support ℝ.
function getintegrandformarginalizationrev(v::Vector{T},
                                        f::Function,
                                        N::Int,
                                        D::Int) where T <: Real
    @assert 0 < N < D
    @assert length(v) == D- length((D-N+1):D)

    g = tt->tt/(1-tt^2)
    u = xx->collect( g(xx[d]) for d = 1:length(xx) )

    J = tt->( (1+tt^2)/(1-tt^2)^2 )

    integrand_func = xx->f( [v; u(xx)] )*prod( J(xx[d]) for d = 1:length(xx) )

    return integrand_func
end
# f is joint pdf.
# based on numerical integration.
function evalmarginalpdfrev(f::Function, v::Vector{Float64}, N::Int, D::Int, max_integral_evals::Int = 1000)

    integrand_func = getintegrandformarginalizationrev(v, f, N, D)

    x_min = collect( -1.0 for n = (D-N+1):D )
    x_max = collect( 1.0 for n = (D-N+1):D )
    (val,err) = HCubature.hcubature(integrand_func, x_min, x_max;
                    norm = LinearAlgebra.norm, rtol = sqrt(eps(Float64)), atol = 0,
                    maxevals = max_integral_evals, initdiv = 1)
    #
    return val
end

# f is the joint pdf.
# marginalizes over the first N variables.
# assumes the N variables have support ℝ.
function getintegrandforZ(  f::Function,
                            D::Int) where T <: Real

    g = tt->tt/(1-tt^2)
    u = xx->collect( g(xx[d]) for d = 1:D )

    J = tt->( (1+tt^2)/(1-tt^2)^2 )

    integrand_func = xx->f( u(xx) )*prod( J(xx[d]) for d = 1:D )

    return integrand_func
end

# f is joint pdf.
# based on numerical integration.
function evalZ9(f::Function, D::Int, max_integral_evals::Int = 10000)

    integrand_func = getintegrandforZ(f, D)

    x_min = collect( -1.0 for d = 1:D )
    x_max = collect( 1.0 for d = 1:D )
    (val,err) = HCubature.hcubature(integrand_func, x_min, x_max;
                    norm = LinearAlgebra.norm, rtol = sqrt(eps(Float64)), atol = 0,
                    maxevals = max_integral_evals, initdiv = 1)
    #
    return val, err
end

# Uses the change-of-variable here:
# https://github.com/stevengj/cubature/blob/master/README.md#infinite-intervals
function evalindefiniteintegral(f::Function, max_integral_evals = typemax(Int))

    integrand_func = tt->(f(tt[1]/(1-tt[1]^2))*(1+tt[1]^2)/(1-tt[1]^2)^2)
    x_min = -one(Float64)
    x_max = one(Float64)
    (val,err) = HCubature.hcubature(integrand_func, [x_min], [x_max];
                    norm = LinearAlgebra.norm, rtol = sqrt(eps(Float64)), atol = 0,
                    maxevals = max_integral_evals, initdiv=1)

    return val, err
end

# for p(A | B) as a function of f(A,B).
function setupCDFkernel(hyperparam_A::T,
                        hyperparam_canonical_B::T,
                        ϕ_map::Function)::CPDF1DKernelType{T,Spline34KernelType{T}} where T <: Real

    θ_A = RationalQuadraticKernelType(hyperparam_A)

    θ_canonical = Spline34KernelType(hyperparam_canonical_B) # for adaptive kernel.
    θ_B = AdaptiveKernelType(θ_canonical, xx->ϕ_map(xx))

    return CPDF1DKernelType(θ_A, θ_B)
end

max_stars_repo_name: RoyCCWang/RKHSRegularization

max_issues_repo_name: RoyCCWang/RKHSRegularization

max_forks_repo_name: RoyCCWang/RKHSRegularization


text:
export SplitExplicitSolver

using ..BalanceLaws:
    initialize_states!,
    tendency_from_slow_to_fast!,
    cummulate_fast_solution!,
    reconcile_from_fast_to_slow!

LSRK2N = LowStorageRungeKutta2N

@doc """
    SplitExplicitSolver(slow_solver, fast_solver; dt, t0 = 0, coupled = true)

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

```math
  \\dot{Q_{fast}} = f_{fast}(Q_{fast}, Q_{slow}, t)
  \\dot{Q_{slow}} = f_{slow}(Q_{slow}, Q_{fast}, t)
```

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

This method performs an operator splitting to timestep the vertical average
of the model at a faster rate than the full model. This results in a first-
order time stepper.

""" SplitExplicitSolver
mutable struct SplitExplicitSolver{SS, FS, RT, MSA} <: AbstractODESolver
    "slow solver"
    slow_solver::SS
    "fast solver"
    fast_solver::FS
    "time step"
    dt::RT
    "time"
    t::RT
    "elapsed time steps"
    steps::Int
    "storage for transfer tendency"
    dQ2fast::MSA

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

        dQ2fast = similar(slow_solver.dQ)
        dQ2fast .= -0
        MSA = typeof(dQ2fast)

        return new{SS, FS, RT, MSA}(
            slow_solver,
            fast_solver,
            RT(dt),
            RT(t0),
            0,
            dQ2fast,
        )
    end
end

function dostep!(
    Qslow,
    split::SplitExplicitSolver{SS},
    param,
    time::Real,
) where {SS <: LSRK2N}
    slow = split.slow_solver
    fast = split.fast_solver

    Qfast = slow.rhs!.modeldata.Q_2D

    dQslow = slow.dQ
    dQ2fast = split.dQ2fast

    slow_bl = slow.rhs!.balance_law
    fast_bl = fast.rhs!.balance_law

    groupsize = 256

    slow_dt = getdt(slow)
    fast_dt_in = getdt(fast)

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

        # Initialize fast model and tendency adjustment
        # before evalution of slow mode
        initialize_states!(slow_bl, fast_bl, slow.rhs!, fast.rhs!, Qslow, Qfast)

        # Evaluate the slow mode
        # --> save tendency for the fast
        slow.rhs!(dQ2fast, Qslow, param, slow_stage_time, increment = false)

        # vertically integrate slow tendency to advance fast equation
        # and use vertical mean for slow model (negative source)
        # ---> work with dQ2fast as input
        tendency_from_slow_to_fast!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
            dQ2fast,
        )

        # Compute (and RK update) slow tendency
        slow.rhs!(dQslow, Qslow, param, slow_stage_time, increment = true)

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

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

        updatedt!(fast, fast_dt)

        for substep in 1:nsubsteps
            fast_time = slow_stage_time + (substep - 1) * fast_dt
            dostep!(Qfast, fast, param, fast_time)
            cummulate_fast_solution!(
                slow_bl,
                fast_bl,
                fast.rhs!,
                Qfast,
                fast_time,
                fast_dt,
                substep,
            )
        end

        # Update (RK-stage) slow state
        event = Event(array_device(Qslow))
        event = update!(array_device(Qslow), groupsize)(
            realview(dQslow),
            realview(Qslow),
            slow.RKA[slow_s % length(slow.RKA) + 1],
            slow.RKB[slow_s],
            slow_dt,
            nothing,
            nothing,
            nothing;
            ndrange = length(realview(Qslow)),
            dependencies = (event,),
        )
        wait(array_device(Qslow), event)

        # reconcile slow equation using fast equation
        reconcile_from_fast_to_slow!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
        )
    end
    updatedt!(fast, fast_dt_in)

    return nothing
end

ext: jl


text:
using Test
using BigArrays
using BigArrays.GSDicts
using OffsetArrays

@testset "\ntest 3D image reading and saving" begin
    d = GSDict( "gs://seunglab/jpwu/test/image" )
    ba = BigArray(d; mode=:sequential)

    @show ndims(ba)
    @show size(ba)
    @show eltype(ba)
    a = rand(UInt8, 256,256,16)

    # ba = BigArray(d, UInt8, (128,128,8))
    @time ba[257:512, 257:512, 17:32] = a
    # BigArrays.mysetindex!(ba, a, (201:400, 201:400, 161:116))
    @time b = ba[257:512, 257:512, 17:32]
    @test all(a.==b |> parent)
end

@testset "\ntest UInt32 segmentation" begin
    d = GSDict( "gs://seunglab/jpwu/test/segmentation/" )
    # ba = BigArray( d, configDict )
    ba = BigArray(d)
    a = rand(UInt32, 256,256,16)
    @time ba[257:512, 257:512, 17:32] = a
    @time b = ba[257:512, 257:512, 17:32]
    @test all(a.==b |> parent)
end

@testset "\ntest UInt64 segmenation with uint64" begin
    d = GSDict( "gs://seunglab/jpwu/test/segmentation-uint64/" )
    ba = BigArray(d)
    a = rand(UInt64, 256,256,16)
    @time ba[257:512, 257:512, 17:32] = a
    @time b = ba[257:512, 257:512, 17:32]
    @test all(a.==b |> parent)
end

@testset "test affinity map" begin
    d = GSDict( "gs://seunglab/jpwu/test/affinitymap/" )
    ba = BigArray(d; mode=:sequential)

    a = rand(Float32, 256,256,16,3)
    @time ba[257:512, 257:512, 17:32, 1:3] = a
    @time b = ba[257:512, 257:512, 17:32, 1:3] |> parent

    @show size(a)
    @test all(a.==b)
end

@testset "test semantic map" begin
    d = GSDict( "gs://seunglab/jpwu/test/semanticmap/" )
    a = rand(Float32, 256,256,16,4)
    ba = BigArray(d)

    @time ba[257:512, 257:512, 17:32, 1:4] = a
    @time b = ba[257:512, 257:512, 17:32, 1:4] |> parent

    @show size(a)
    @show size(b)
    @test all(a.==b)
end


text:
using Distributed
# @everywhere using DelimitedFiles
@everywhere using DifferentialEquations
@everywhere using ProgressMeter
@everywhere using MATLAB
@everywhere using Printf

include("BASIN_PSS_function_parallel.jl")

# @everywhere using PyCall
# pygui(:qt)
# @everywhere using PyPlot
# pygui(true)

@everywhere max_hit=5
@everywhere Energy=.25
@everywhere Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
@everywhere ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
@everywhere ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
@everywhere function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=0
return
end
@everywhere condition(u,t,integrator)= u[5]>max_hit
@everywhere function condition2(u,t,integrator) # Event when event_f(u,t) == 0
u[1]
end
@everywhere affect!(integrator) = terminate!(integrator)
@everywhere function affect2!(integrator)
integrator.u[5]=integrator.u[5]+1
end
@everywhere cb2 = ContinuousCallback(condition2,affect2!,rootfind=false)
@everywhere cb1 = DiscreteCallback(condition,affect!)

@everywhere cb=CallbackSet(cb2,cb1)
@everywhere H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))
@everywhere ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
@everywhere ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
@everywhere function Eq_of_M(du,u,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=0
return
end

Q_bound=zeros(0)
P_bound=zeros(0)

Q_0=zeros(0)
P_0=zeros(0)

Q_1=zeros(0)
P_1=zeros(0)

Q_2=zeros(0)
P_2=zeros(0)

Q_3=zeros(0)
P_3=zeros(0)

Q_4=zeros(0)
P_4=zeros(0)

Q_5=zeros(0)
P_5=zeros(0)

# ArrQ[k]=.1*rand()
# ArrP[j]=.01*rand()
@everywhere Q_start=-1.5
@everywhere Q_end=1.5
@everywhere P_start=-2.5
@everywhere P_end=2.5
@everywhere n_iter_P=2000
@everywhere n_iter_Q=2001
@everywhere N=n_iter_P*n_iter_Q

@everywhere delta_P=(P_end-P_start)/n_iter_P
@everywhere delta_Q=(Q_end-Q_start)/n_iter_Q

@everywhere ArrP=P_start: delta_P:P_end
@everywhere ArrQ=Q_start: delta_Q: Q_end

Energy=.25
location="/mnt/bdd38f66-9ece-451a-b915-952523c139d2/Escape/"
h=replace(@sprintf("%.13f",Energy),"."=>"_")
file_name=location*"Escape_"*h*".fig"

T=@showprogress pmap(BASIN_PSS_function_parallel,1:N)
for i=1:N
    Q=ArrQ[Int(ceil(i/n_iter_Q))]
    P=ArrP[mod(i,n_iter_P)+1]
    if floor(T[i])==0
        push!(Q_0,Q)
        push!(P_0,P)
    elseif floor(T[i])==1
        push!(Q_1,Q)
        push!(P_1,P)
    elseif floor(T[i])==2
        push!(Q_2,Q)
        push!(P_2,P)
    elseif floor(T[i])==3
        push!(Q_3,Q)
        push!(P_3,P)
    elseif floor(T[i])==4
        push!(Q_4,Q)
        push!(P_4,P)
    elseif floor(T[i])==5
        push!(Q_5,Q)
        push!(P_5,P)
    elseif T[i]==-1

        push!(Q_bound,Q)
        push!(P_bound,P)
    end
end

mat"figure(); hold on;"
mat"axis([ -$width,$width,-$height,$height ])"
# mat"plot($P_bound,$Q_bound ,'k.','MarkerSize',10)"
mat"plot($P_0,$Q_0 ,'b.','MarkerSize',$size)"
mat"plot($P_1,$Q_1 ,'r.','MarkerSize',$size)"
mat"plot($P_2,$Q_2 ,'g.','MarkerSize',$size)"
mat"plot($P_4,$Q_4 ,'c.','MarkerSize',$size)"
mat"plot($P_5,$Q_5 ,'y.','MarkerSize',$size)"
mat"savefig($file_name)"
mat"close"


text:
using RecipesBase

@recipe function f(filt::PhotometricFilter)
    label --> filt.name
    yguide --> "throughput"
    xguide --> "wavelength"

    wave(filt), throughput(filt)
end

@recipe function f(filts::Vector{<:PhotometricFilter})
    yguide --> "throughput"
    xguide --> "wavelength"

    for filt in filts
        @series begin
            label --> filt.name
            wave(filt), throughput(filt)
        end
    end
end


text:
#####
##### Sequential Access Kernels
#####

function vector_sum(
        A::AbstractArray{T},
        ::Val{N},
        aligned = Val(false),
        nontemporal = Val(false)
    ) where {T, N}

    unroll = 4

    # Make sure we can successfully chunk up this array into SIMD sizes
    @assert iszero(mod(length(A), unroll * N))

    # If we've passed alignment flags, make sure the base pointer of this array is in fact
    # aligned correctly.
    if aligned == Val{true}()
        @assert iszero(mod(Int(pointer(A)), sizeof(T) * N))
    end

    s = Vec{N,T}(zero(T))
    @inbounds for i in 1:(unroll*N):length(A)
        _v1 = vload(Vec{N,T}, pointer(A, i),         aligned, nontemporal)
        _v2 = vload(Vec{N,T}, pointer(A, i + N),     aligned, nontemporal)
        _v3 = vload(Vec{N,T}, pointer(A, i + 2 * N), aligned, nontemporal)
        _v4 = vload(Vec{N,T}, pointer(A, i + 3 * N), aligned, nontemporal)

        s += _v1 + _v2 + _v3 + _v4
    end
    return sum(s)
end

function vector_write(
        A::AbstractArray{T},
        ::Val{N},
        aligned = Val(false),
        nontemporal = Val(false),
    ) where {T, N}

    unroll = 4

    # Make sure we can successfully chunk up this array into SIMD sizes
    @assert iszero(mod(length(A), unroll * N))

    # If we've passed alignment flags, make sure the base pointer of this array is in fact
    # aligned correctly.
    if aligned == Val{true}()
        @assert iszero(mod(Int(pointer(A)), sizeof(T) * N))
    end

    # We need to do some schenanigans to get LLVM to emit the correct code.
    # Just doing something like
    #    vstore(s, pointer(A, i),         aligned, nontemporal)
    #    vstore(s, pointer(A, i + N),     aligned, nontemporal)
    #    vstore(s, pointer(A, i + 2 * N), aligned, nontemporal)
    #    vstore(s, pointer(A, i + 3 * N), aligned, nontemporal)
    #
    # Results in spurious `mov` instructions between the vector stores for pointer
    # conversions, even though these are really not needed.
    #
    # Instead, we perform the pointer arithmetic manually.
    s = Vec{N,T}(zero(T))
    base = pointer(A)
    @inbounds for i in 0:(unroll*N):(length(A) - 1)
        vstore(s, base + sizeof(T) * i,           aligned, nontemporal)
        vstore(s, base + sizeof(T) * (i + N),     aligned, nontemporal)
        vstore(s, base + sizeof(T) * (i + (2*N)), aligned, nontemporal)
        vstore(s, base + sizeof(T) * (i + (3*N)), aligned, nontemporal)
    end
    return nothing
end

function vector_increment(
        A::AbstractArray{T},
        ::Val{N},
        aligned = Val(false),
        nontemporal = Val(false),
    ) where {T, N}

    unroll = 4

    # Make sure we can successfully chunk up this array into SIMD sizes
    @assert iszero(mod(length(A), unroll * N))

    # If we've passed alignment flags, make sure the base pointer of this array is in fact
    # aligned correctly.
    if aligned == Val{true}()
        @assert iszero(mod(Int(pointer(A)), sizeof(T) * N))
    end

    s = Vec{N,T}(one(T))
    base = pointer(A)
    @inbounds for i in 0:(unroll*N):(length(A) - 1)
        _v1 = vload(Vec{N,T}, base + sizeof(T) * i,           aligned, nontemporal)
        _v2 = vload(Vec{N,T}, base + sizeof(T) * (i + N),     aligned, nontemporal)
        _v3 = vload(Vec{N,T}, base + sizeof(T) * (i + (2*N)), aligned, nontemporal)
        _v4 = vload(Vec{N,T}, base + sizeof(T) * (i + (3*N)), aligned, nontemporal)

        _u1 = _v1 + s
        _u2 = _v2 + s
        _u3 = _v3 + s
        _u4 = _v4 + s

        vstore(_u1, base + sizeof(T) * i,           aligned, nontemporal)
        vstore(_u2, base + sizeof(T) * (i + N),     aligned, nontemporal)
        vstore(_u3, base + sizeof(T) * (i + (2*N)), aligned, nontemporal)
        vstore(_u4, base + sizeof(T) * (i + (3*N)), aligned, nontemporal)
    end
    return nothing
end

#####
##### Random Access Kernels
#####

function hop_sum(A::AbstractArray{T}, ::Val{N}, lfsr, nontemporal = Val{false}()) where {T,N}
    return _hop_sum(reinterpret(Vec{N,T}, A), lfsr, nontemporal)
end

@inline function _hop_sum(A::AbstractArray{Vec{N,T}}, lfsr, nontemporal = Val{false}()) where {N,T}
    s = zero(eltype(A))
    base = Base.unsafe_convert(Ptr{T}, pointer(A))
    @inbounds for i in lfsr
        ptr = base + sizeof(eltype(A)) * (i-1)
        v = vload(Vec{N,T}, ptr, Val{true}(), nontemporal)
        s += v
    end
    return s
end

function hop_write(A::AbstractArray{T}, ::Val{N}, lfsr, nontemporal = Val{false}()) where {T,N}
    return _hop_write(reinterpret(Vec{N,T}, A), lfsr, nontemporal)
end

@inline function _hop_write(A::AbstractArray{Vec{N,T}}, lfsr, nontemporal = Val{false}()) where {N,T}
    s = zero(eltype(A))
    base = Base.unsafe_convert(Ptr{T}, pointer(A))
    @inbounds for i in lfsr
        ptr = base + sizeof(eltype(A)) * (i-1)
        vstore(s, ptr, Val{true}(), nontemporal)
    end
    return nothing
end

function hop_increment(A::AbstractArray{T}, ::Val{N}, lfsr, nontemporal = Val{false}()) where {T,N}
    return _hop_increment(reinterpret(Vec{N,T}, A), lfsr, nontemporal)
end

@inline function _hop_increment(A::AbstractArray{Vec{N,T}}, lfsr, nontemporal = Val{false}()) where {N,T}
    base = Base.unsafe_convert(Ptr{T}, pointer(A))
    @inbounds for i in lfsr
        ptr = base + sizeof(eltype(A)) * (i-1)
        v = vload(Vec{N,T}, ptr, Val{true}(), nontemporal)
        vstore(v + one(eltype(A)), ptr, Val{true}(), nontemporal)
    end
end

#####
##### Streamed copy from A to B
#####

function stream_copy!(A::AbstractArray{T}, B::AbstractArray{T}, valn::Val{N}) where {T, N}
    unroll = 4

    # Alignment and size checking
    @assert iszero(mod(length(A), unroll * N))
    @assert iszero(mod(Int(pointer(A)), sizeof(T) * N))
    @assert iszero(mod(Int(pointer(B)), sizeof(T) * N))

    # Forward to the one without bounds checking so we can checkout out the generated
    # code more easily.
    return unsafe_stream_copy!(A, B, valn)
end

function unsafe_stream_copy!(A::AbstractArray{T}, B::AbstractArray{T}, ::Val{N}) where {N,T}
    # These are streaming stores after all.
    unroll = 4
    aligned = Val{true}()
    nontemporal = Val{true}()

    # Again, this pointer arithmetic thing seems to be necessary to get the best native code.
    pa = pointer(A)
    pb = pointer(B)
    @inbounds for i in 0:(unroll*N):(length(A) - 1)
        _v1 = vload(Vec{N,T}, pb + sizeof(T) * i,           aligned, nontemporal)
        _v2 = vload(Vec{N,T}, pb + sizeof(T) * (i + N),     aligned, nontemporal)
        _v3 = vload(Vec{N,T}, pb + sizeof(T) * (i + (2*N)), aligned, nontemporal)
        _v4 = vload(Vec{N,T}, pb + sizeof(T) * (i + (3*N)), aligned, nontemporal)

        vstore(_v1, pa + sizeof(T) * i,           aligned, nontemporal)
        vstore(_v2, pa + sizeof(T) * (i + N),     aligned, nontemporal)
        vstore(_v3, pa + sizeof(T) * (i + (2*N)), aligned, nontemporal)
        vstore(_v4, pa + sizeof(T) * (i + (3*N)), aligned, nontemporal)
    end
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
"""
# Summary ScaleChannels
    mutable struct ScaleChannels <: Scaler
# Fields
    scale::VarOrNil
    views::Tuple
Applies vector multiplication over a N-dimensional input.
This scalar called `scale` is a learnable vector parameter.
"""
mutable struct ScaleChannels <: Scaler
    scale::VarOrNil
    views::Tuple
    function ScaleChannels(scalar::AbstractFloat;
                           ndims::Int,
                           keptdim::Int,
                           keptsize::Int,
                           type::Type=Array{Float32})
        @assert ndims > 0 "ndims > 0, but got ndims=$ndims"
        @assert ndims >= keptdim >= 1 "keptdim in [1, $ndims], but got keptdim=$keptdim"

        shape = ntuple(i -> i==keptdim ? keptsize : 1, ndims);
        views = ntuple(i -> i>=keptdim ? i+1 : i, ndims-1);

        scale = Variable{type}(Zeros(type, shape) .+ eltype(type)(scalar), true, true, true);
        new(scale, views)
    end
    function ScaleChannels(views::Tuple)
        new(nothing, views)
    end
end

function clone(this::ScaleChannels; type::Type=Array{Float32})
    cloned = ScaleChannels(this.views)
    cloned.scale = clone(this.scale, type=type)
    return cloned
end

function Base.show(io::IO, m::ScaleChannels)
    print(io, "ScaleChannels(size(scale)=$(size(m.scale)); type=$(typeof(m.scale.value)))")
end

function paramsof(m::ScaleChannels)
    params = Vector{Variable}(undef,1)
    params[1] = m.scale
    return params
end

function xparamsof(m::ScaleChannels)
    xparams = Vector{XVariable}(undef,1)
    xparams[1] = ('w', m.scale)
    return xparams
end

function nparamsof(m::ScaleChannels)
    return 1
end

function bytesof(m::ScaleChannels, unit::String="MB")
    return blocksize(sizeof(m.scale), uppercase(unit))
end

function forward(m::ScaleChannels, x::Variable{T}) where T
    k = m.scale
    y = Variable{T}(ᵛ(x) .* ᵛ(k), x.backprop)

    if x.backprop
        function ScaleChannelsBackward()
            if need2computeδ!(x) δ(x) .+=     δ(y) .* ᵛ(k)                end
            if need2computeδ!(k) δ(k) .+= sum(δ(y) .* ᵛ(x), dims=m.views) end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, ScaleChannelsBackward)
    end
    return y
end

function predict(m::ScaleChannels, x::AbstractArray)
    k = m.scale.value
    return x .* k
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using Test, StatsModels
using DataFrames
using StatsBase

using Unfold
include("test_utilities.jl")

data, evts = loadtestdata("test_case_3a") #
data_r = reshape(data, (1, :))

data_e, times = Unfold.epoch(data = data_r, tbl = evts, τ = (0.0, 1), sfreq = 20) # cut the data into epochs
basisfunction = firbasis(τ = (0.0, 1), sfreq = 20, name = "basisA")

# f  = @formula 0~ 1 * beta_0 + contA*beta_1
# fit => beta_0 = -1.5, beta_1 = 2
# predict => (contA=0.5) y_hat = 1 * -1.5 + 0.5 * 2 = -0.5

f = @formula 0 ~ 1 # 1
m_mul = fit(UnfoldModel, f, evts, data_e, times)
m_tul = fit(UnfoldModel, f, evts, data_r, basisfunction)
m_mul_results = coeftable(m_mul)
m_tul_results = coeftable(m_tul)

conditionA = [0, 1.0]
continuousA = [-1.0, 0, 1.0]

tmp = reshape(
    [[x, y] for x in conditionA, y in continuousA],
    length(conditionA) * length(continuousA),
)
evts_grid = DataFrame(collect(hcat(tmp...)'), ["conditionA", "continuousA"])

yhat_tul = predict(m_tul, evts_grid)
yhat_mul = predict(m_mul, evts_grid, times)

@test unique(yhat_mul.times) == times

@test size(m_tul_results)[1] * size(evts_grid)[1] == size(yhat_tul)[1]
@test size(m_mul_results)[1] * size(evts_grid)[1] == size(yhat_mul)[1]

@test yhat_mul[end-3, :yhat] ≈ mean(data[data.!=0])
@test yhat_tul[end-3, :yhat] ≈ mean(data[data.!=0])

@test yhat_mul.times[1] == 0.0

f = @formula 0 ~ 1 + conditionA + continuousA# 1
m_tul = fit(UnfoldModel, f, evts, data_r, basisfunction)
m_mul = fit(UnfoldModel, f, evts, data_e, times)

m_mul_results = coeftable(m_mul)
m_tul_results = coeftable(m_tul)

yhat_tul = predict(m_tul, evts_grid)
yhat_mul = predict(m_mul, evts_grid, times)
#@test predict(m_mul,evts).yhat ≈ yhat.yhat

@test isapprox(sum(coef(m_tul)[[5, 25, 45]]), yhat_tul.yhat[end-3], atol = 0.001)

@test isapprox(
    yhat_tul[yhat_tul.times.==0.5, :yhat],
    yhat_mul[yhat_mul.times.==0.5, :yhat],
    atol = 0.001,
)
@test isapprox(yhat_tul[yhat_tul.times.==0.5, :yhat], [-2, 1, 2, 5, 6, 9.0], atol = 0.0001)

if 1 == 0
    #for visualization
    using AlgebraOfGraphics
    yhat_mul.conditionA = categorical(yhat_mul.conditionA)
    m = mapping(:times, :yhat, color = :continuousA, linestyle = :conditionA)
    df = yhat_mul
    AlgebraOfGraphics.data(df) * visual(Lines) * m |> draw
end


text:
#####
##### Diffusivities
#####

@inline κᶠᶜᶜ(i, j, k, grid, clock, κ::Number) = κ
@inline κᶜᶠᶜ(i, j, k, grid, clock, κ::Number) = κ
@inline κᶜᶜᶠ(i, j, k, grid, clock, κ::Number) = κ

@inline κᶠᶜᶜ(i, j, k, grid, clock, κ::F) where F<:Function = κ(xnode(Face(),   i, grid), ynode(Center(), j, grid), znode(Center(), k, grid), clock.time)
@inline κᶜᶠᶜ(i, j, k, grid, clock, κ::F) where F<:Function = κ(xnode(Center(), i, grid), ynode(Face(),   j, grid), znode(Center(), k, grid), clock.time)
@inline κᶜᶜᶠ(i, j, k, grid, clock, κ::F) where F<:Function = κ(xnode(Center(), i, grid), ynode(Center(), j, grid), znode(Face(),   k, grid), clock.time)

# Assumes that `κ` is located at cell centers
@inline κᶠᶜᶜ(i, j, k, grid, clock, κ::AbstractArray) = ℑxᶠᵃᵃ(i, j, k, grid, κ)
@inline κᶜᶠᶜ(i, j, k, grid, clock, κ::AbstractArray) = ℑyᵃᶠᵃ(i, j, k, grid, κ)
@inline κᶜᶜᶠ(i, j, k, grid, clock, κ::AbstractArray) = ℑzᵃᵃᶠ(i, j, k, grid, κ)

#####
##### Convenience diffusive flux function
#####

@inline diffusive_flux_x(i, j, k, grid, clock, κ, c) = - κᶠᶜᶜ(i, j, k, grid, clock, κ) * ∂xᶠᶜᵃ(i, j, k, grid, c)
@inline diffusive_flux_y(i, j, k, grid, clock, κ, c) = - κᶜᶠᶜ(i, j, k, grid, clock, κ) * ∂yᶜᶠᵃ(i, j, k, grid, c)
@inline diffusive_flux_z(i, j, k, grid, clock, κ, c) = - κᶜᶜᶠ(i, j, k, grid, clock, κ) * ∂zᵃᵃᶠ(i, j, k, grid, c)

#####
##### Diffusive flux divergence
#####

"""
    ∇_dot_qᶜ(i, j, k, grid, clock, closure::AbstractTurbulenceClosure, c, ::Val{tracer_index}, args...)

Calculates the divergence of the diffusive flux `qᶜ` for a tracer `c` via

    1/V * [δxᶜᵃᵃ(Ax * diffusive_flux_x) + δyᵃᶜᵃ(Ay * diffusive_flux_y) + δzᵃᵃᶜ(Az * diffusive_flux_z)]

which will end up at the location `ccc`.
"""
@inline function ∇_dot_qᶜ(i, j, k, grid, closure::AbstractTurbulenceClosure, c, tracer_index, args...)

    disc = time_discretization(closure)

    return 1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_uᶠᶜᶜ, diffusive_flux_x, disc, closure, c, tracer_index, args...) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_vᶜᶠᶜ, diffusive_flux_y, disc, closure, c, tracer_index, args...) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_wᶜᶜᵃ, diffusive_flux_z, disc, closure, c, tracer_index, args...))
end

#####
##### Gradients of Laplacians
#####

@inline ∂x_∇²h_cᶠᶜᶜ(i, j, k, grid, c) = 1 / Azᶠᶜᵃ(i, j, k, grid) * δxᶠᵃᵃ(i, j, k, Δy_cᶜᶜᵃ, ∇²hᶜᶜᶜ, c)
@inline ∂y_∇²h_cᶜᶠᶜ(i, j, k, grid, c) = 1 / Azᶜᶠᵃ(i, j, k, grid) * δyᵃᶠᵃ(i, j, k, Δx_cᶜᶜᵃ, ∇²hᶜᶜᶜ, c)


text:
module Algorithm

include("../compat.jl")

using OpenGene

export edit_distance,
    hamming_distance,
    overlap,
    simple_merge,
    try_merge

include("distance.jl")
include("overlap.jl")

end


text:
using ASE
using Base.Test

println("-------------------")
println(" Testing ASE.jl")
println("-------------------")

@testset "ASE" begin
   @testset "Build" begin include("test_build.jl"); end
   @testset "Atoms"  begin include("testase.jl"); end
   @testset "Calculators" begin include("test_calculators.jl"); end
end


text:
export Network

"""
# Network
Network structs are used to run simulations. It contains all neurons and connections.
# Arguments
- `neurons::NeuronType`
- `batch_size::Integer=1`
...
- `weight::AbstractArray`: Final size= (#neurons, #neurons)
- `adjacency::AbstractArray`: Final size= (#neurons, #neurons)
- `spikes::AbstractArray`: Final size= (Batch size, #neurons)
- `voltage::AbstractArray`: Final size= (Batch size, #neurons)
- `recovery::AbstractArray`: Membrane recovery variable, only applicable to Izhikevich neurons
- `refractory::AbstractArray`: Refractory period, only applicable to LIF neurons. Final size= (Batch size, #neurons)
...
"""
@with_kw mutable struct Network{T<:NeuronType} <: SpikingNetwork{T}
        neurons::T
        batch_size::Int = 1
        learning_rule::Union{LearningRule,Nothing} = nothing
        weight::AbstractArray = Array{Float64}(undef, (0, 0))
        adjacency::AbstractArray = Array{Bool}(undef, (0, 0))
        spikes::AbstractArray = zeros(Bool, batch_size, 0)
        voltage::AbstractArray = ones(Float32, batch_size, 0)
        recovery::AbstractArray = ones(Float32, batch_size, 0)
        refractory::AbstractArray = zeros(Int32, batch_size, 0)
        learning::Bool = isnothing(learning_rule) ? false : true
        groups::Dict{String,NeuronGroup} = Dict{String,NeuronGroup}()
end

Network(neurons) = Network(neurons=neurons)

function add_group!(network::SpikingNetwork, N::Int; name::Union{String,Nothing}=nothing)
        Group = NeuronGroup(N, size(network.weight, 1)+1:size(network.weight, 1)+N)
        network.weight = pad2D(network.weight, N)
        network.adjacency = pad2D(network.adjacency, N)
        network.spikes = pad1D(network.spikes, N)
        if !isnothing(network.learning_rule)
                add_group!(network.learning_rule, N, network.batch_size)
        end
        _add_neuron_features!(network, N)
        if isnothing(name)
                name = "group_$(length(network.groups)+1)"
        end
        network.groups[name] = Group
        return Group
end

function _add_neuron_features!(network::SpikingNetwork{LIF}, N::Int)
        network.voltage = pad1D(network.voltage, N)
        fill!(network.voltage, network.neurons.v_rest)
        network.refractory = pad1D(network.refractory, N)
end

function _add_neuron_features!(network::SpikingNetwork{Izhikevich}, N::Int)
        network.voltage = pad1D(network.voltage, N)
        if typeof(network.neurons.c) <: AbstractArray
                network.voltage = network.neurons.c[:, 1:size(network.voltage, 2)]
        else
                fill!(network.voltage, network.neurons.c)
        end
        if typeof(network.neurons.b) <: AbstractArray
                network.recovery = network.neurons.b[:, 1:size(network.voltage, 2)] .* network.voltage
        else
                network.recovery = network.neurons.b .* network.voltage
        end
end

"""
Add a connection between two neuron groups in the network
"""
function connect!(
        network::SpikingNetwork,
        source::NeuronGroup,
        target::NeuronGroup,
        weight::AbstractArray;
        adjacency::AbstractArray=ones(Bool, source.n, target.n)
)
        network.weight[source.idx, target.idx] = weight
        network.adjacency[source.idx, target.idx] = adjacency
        return
end

function connect!(
        network::SpikingNetwork,
        source::NeuronGroup,
        target::NeuronGroup,
        connection::Connection,
)
        connect!(network, source, target, connection.weight; adjacency=connection.adjacency)
end

function run!(
        network::Network;
        input_spikes::Union{AbstractMatrix{Bool},Nothing}=nothing,
        input_voltage::Union{AbstractMatrix,Nothing}=nothing
)
        # Evoke spikes
        network.spikes = network.voltage .>= network.neurons.v_thresh
        # External spikes
        if !isnothing(input_spikes)
                network.spikes = network.spikes .| input_spikes
        end
        # Supress neurons #TODO
        _update!(network, input_voltage)
        # Learning process
        if network.learning # Apply the learning rule and update weights
                train!(network)
        end
        return network.spikes, network.voltage
end

function _update!(
        network::Network{LIF},
        input_voltage::Union{AbstractMatrix,Nothing}
)
        # Decay voltages.
        network.voltage = (
                network.neurons.voltage_decay_factor .*
                (network.voltage .- network.neurons.v_rest)
                .+
                network.neurons.v_rest
        )
        # Get current
        current = network.spikes * network.weight  # + network.bias #TODO
        # External voltage:
        if !isnothing(input_voltage)
                input_voltage[network.refractory.>0] .= 0.0
                network.voltage += input_voltage
        end
        # update voltages
        network.voltage += current
        network.voltage[network.refractory.>0] .= network.neurons.v_rest # reset the voltage of the neurons in the refractory period
        network.voltage[network.spikes] .= network.neurons.v_reset  # change the voltage of spiked neurons to v_reset
        # Update refractory timepoints
        network.refractory .-= network.neurons.dt
        network.refractory[network.spikes] .= network.neurons.refractory_period
end

function _update!(
        network::Network{Izhikevich},
        input_voltage::Union{AbstractMatrix,Nothing}
)
        # Get current
        current = network.spikes * network.weight
        # External voltage:
        if !isnothing(input_voltage)
                current += input_voltage
        end
        # Update voltages (spiked neurons)
        network.voltage[network.spikes] .= (network.spikes.*network.neurons.c)[network.spikes]  # change the voltage of spiked neurons to c #TODO: optimize for scalar values
        network.recovery[network.spikes] .+= (network.spikes.*network.neurons.d)[network.spikes]  # add d to the recovery parameter of spiked neurons #TODO: optimize for scalar values
        # Update voltages
        for _ in 1:network.neurons.runge_kutta_order
                network.voltage += (network.neurons.dt / network.neurons.runge_kutta_order) .* (
                        0.04 .* network.voltage .^ 2 + 5 .* network.voltage .+ 140 - network.recovery + current
                )
        end
        network.recovery += network.neurons.dt .* (
                network.neurons.a .* (network.neurons.b .* network.voltage - network.recovery)
        )
end

function _reset!(network::SpikingNetwork)
        fill!(network.spikes, 0)
        if !isnothing(network.learning_rule)
                reset!(network.learning_rule)
        end
        return
end

function reset!(network::SpikingNetwork{LIF})
        fill!(network.voltage, network.neurons.v_rest)
        fill!(network.refractory, 0)
        _reset!(network)
end

function reset!(network::SpikingNetwork{Izhikevich})
        network.voltage .= network.neurons.c
        network.recovery = network.neurons.b .* network.voltage
        _reset!(network)
end

function makeInput(network::SpikingNetwork, time::Integer, inputs::Dict{NeuronGroup}, type=Float64)
        input = zeros(type, time, size(network.weight)[1])
        for group_input in inputs
                input[:, group_input[1].idx] .= group_input[2]
        end
        return input
end

function save(network::SpikingNetwork, filename::AbstractString)
        save_object(filename, network |> cpu)
end

function Base.getindex(network::SpikingNetwork, idx::Union{UnitRange{Int},Vector{Int},Int})
        typeof(network)(
                [
                        begin
                                if typeof(getfield(network, property)) <: AbstractArray
                                        if length(getfield(network, property)) > 0
                                                getfield(network, property)[idx]
                                        else
                                                getfield(network, property)
                                        end
                                else
                                        getfield(network, property)
                                end
                        end for property in propertynames(network)
                ]...
        )
end

function Base.getindex(network::SpikingNetwork, group::NeuronGroup)
        return network[group.idx]
end


text:
module TestFreqResp
using CustomTest
using Base.Test
using ControlSystems

## EVALFR ##
H = [tf(0) tf([3, 0],[1, 1, 10]) ; tf([1, 1],[1, 5]) tf([2],[1, 6])]
G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2],
       [0 3 0 0; -2 0 0 1], [0 0; 1 0])

@test evalfr(H, -6) == [0.0 -0.45; 5.0 Inf]
@test evalfr(H, -5) == [0.0 -0.5; Inf 2.0]
@test evalfr(H, -1) == [0.0 -0.3; 0.0 0.4]
@test_approx_eq evalfr(H, 0) [0.0 0.0; 0.2 1/3]

@test evalfr(G, -6) == [Inf Inf; Inf Inf]
@test evalfr(G, -5) == [Inf Inf; Inf Inf]
@test_approx_eq evalfr(G, -1) [0.0 -0.3; 0.0 0.4]
@test_approx_eq evalfr(G, 0) [0.0 0.0; 0.2 1/3]

## Shortcut notation for evalfr ##
F = tf([1],[1,0.5],-1)
omega = 2
z = 0.5(1+im)
@test F(omega,true)[1] == 1/(exp(-im*2)+0.5)
@test F(z,false)[1] == 1/(z+0.5)
@test_throws ErrorException F(z,true)

## Test bode, nyquist and sigma
sys = [tf([1,-1], [1,1,1]) 0; 0 tf([1],[1,1])]
f(s) = [(s-1)./(s.^2+s+1) 0; 0 1./(1+s)]
ws = logspace(-2,2,50)
resp = Array(Complex128,50,2,2)
for (i,w) in enumerate(ws)
    resp[i,:,:] = f(im*w)
end

@test bode(sys, ws)[1:2] == (abs(resp), rad2deg(angle(resp)))
@test nyquist(sys, ws)[1:2] == (real(resp), imag(resp))
sigs = Array(Float64,50,2)
for i in eachindex(ws)
    sigs[i,:] =  svdvals(resp[i,:,:])
end
@test sigma(sys, ws)[1] == sigs

#Test default freq vector contains at least half a decade more than all features
p, z = 100, 1/1000
sys2 = [tf([1],[1/p,1]) tf([1/z, 2/z, 1],[1])]
mag, phase, ws2 = bode(sys2)
@test maximum(ws2) >= 2max(p,z)
@test minimum(ws2) <= 0.5min(p,z)

#Test default freq vector not too small
p, z = 1, 1.2
sys2 = [tf([1],[1/p,1]) tf([1/z, 2/z, 1],[1])]
mag, mag, ws2 = bode(sys2)
@test maximum(ws2) >= 5max(p,z)
@test minimum(ws2) <= 0.2min(p,z)
@test length(ws2) > 100
end


text:
const queue_map = Dict(
    :compute => QueueCompute(),
    :graphics => QueueGraphics(),
    :transfer => QueueTransfer(),
    :sparse_binding => QueueSparseBinding(),
    :decode => QueueVideoDecode(),
    :encode => QueueVideoEncode(),
)

const render_pass_compatibiltiy_map = Dict(
    :both => [RenderPassInside(), RenderPassOutside()],
    :inside => [RenderPassInside()],
    :outside => [RenderPassOutside()],
)

Base.:(==)(x::S, y::S) where {S<:Spec} = all(name -> getproperty(x, name) == getproperty(y, name), fieldnames(S))

print_parent_info(io::IO, spec::Union{SpecStruct,SpecFunc}, props) =
    println(io, join(string.(vcat(typeof(spec), spec.name, props)), ' '))

children(spec::SpecStruct) = spec.members
children(spec::SpecFunc) = spec.params

print_children(io::IO, spec::Union{SpecStruct,SpecFunc}) = println.(Ref(io), string.("    ", children(spec)))

function Base.show(io::IO, spec::SpecStruct)
    props = string.([spec.type])
    spec.is_returnedonly && push!(props, "returnedonly")
    !isempty(spec.extends) && push!(props, "extends $(spec.extends)")

    print_parent_info(io, spec, props)
    print_children(io, spec)
end

function Base.show(io::IO, spec::SpecFunc)
    props = string.([spec.type, spec.return_type])
    !isempty(spec.render_pass_compatibility) && push!(
        props,
        string("to be executed ", join(map(spec.render_pass_compatibility) do compat
            @match compat begin
                ::RenderPassInside => "inside"
                ::RenderPassOutside => "outside"
            end
        end, " and "), " render passes"),
    )
    !isempty(spec.queue_compatibility) &&
        push!(props, string(" (compatible with $(join(string.(spec.queue_compatibility), ", ")) queues)"))

    print_parent_info(io, spec, props)
    print_children(io, spec)
end

function Base.show(io::IO, spec::Union{SpecStructMember,SpecFuncParam})
    props = string.([spec.requirement])
    spec.is_constant && push!(props, "constant")
    spec.is_externsync && push!(props, "externsync")
    is_arr(spec) && push!(props, "with length $(spec.len)")
    !isempty(spec.arglen) && push!(props, string("length of ", join(spec.arglen, ' ')))

    print(io, join(string.(vcat(typeof(spec), spec.type, spec.name, props)), ' '))
end

function SpecStructMember(node::Node)
    id = extract_identifier(node)
    SpecStructMember(
        parent_name(node),
        id == :module ? :_module : id,
        extract_type(node),
        is_constant(node),
        externsync(node),
        PARAM_REQUIREMENT(node),
        len(node),
        arglen(node, neighbor_type = "member"),
        !parse(Bool, getattr(node, "noautovalidity", default="false", symbol=false)),
    )
end

function SpecStruct(node::Node)
    name_str = node["name"]
    returnedonly = haskey(node, "returnedonly")
    type = @match returnedonly begin
        true => PROPERTY
        _ && if occursin("CreateInfo", name_str)
        end => CREATE_INFO
        _ && if occursin("AllocateInfo", name_str)
        end => ALLOCATE_INFO
        _ && if occursin("Info", name_str)
        end => GENERIC_INFO
        _ => DATA
    end
    extends = @match struct_csv = getattr(node, "structextends", symbol = false) begin
        ::String => Symbol.(split(struct_csv, ','))
        ::Nothing => []
    end
    SpecStruct(
        Symbol(name_str),
        type,
        returnedonly,
        extends,
        StructVector(SpecStructMember.(findall("./member", node))),
    )
end

function SpecUnion(node::Node)
    members = findall("./member", node)
    selectors = getattr.(members, "selection")
    SpecUnion(
        getattr(node, "name"),
        extract_type.(members),
        filter(!isnothing, selectors),
        haskey(node, "returnedonly"),
    )
end

SpecFuncParam(node::Node) = SpecFuncParam(
    parent_name(node),
    extract_identifier(node),
    extract_type(node),
    is_constant(node),
    externsync(node),
    PARAM_REQUIREMENT(node),
    len(node),
    arglen(node),
    !parse(Bool, getattr(node, "noautovalidity", default="false", symbol=false)),
)

function SpecFunc(node::Node)
    name = command_name(node)
    queues = @match getattr(node, "queues", symbol = false) begin
        qs::String => [queue_map[Symbol(q)] for q ∈ split(qs, ',')]
        ::Nothing => []
    end
    rp_reqs = @match getattr(node, "renderpass") begin
        x::Symbol => render_pass_compatibiltiy_map[x]
        ::Nothing => []
    end
    ctype = @match findfirst(startswith.(string(name), ["vkCreate", "vkDestroy", "vkAllocate", "vkFree", "vkCmd"])) begin
        i::Integer => FUNC_TYPE(i)
        if any(startswith.(string(name), ["vkGet", "vkEnumerate"]))
        end => QUERY
        _ => OTHER
    end
    return_type = extract_type(findfirst("./proto", node))
    codes(type) = Symbol.(filter(!isempty, split(getattr(node, type; default = "", symbol = false), ',')))
    SpecFunc(
        name,
        ctype,
        return_type,
        rp_reqs,
        queues,
        StructVector(SpecFuncParam.(findall("./param", node))),
        codes("successcodes"),
        codes("errorcodes"),
    )
end

SpecEnum(node::Node) =
    SpecEnum(Symbol(node["name"]), StructVector(SpecConstant.(findall("./enum[@name and not(@alias)]", node))))

SpecBit(node::Node) = SpecBit(
    Symbol(node["name"]),
    parse(Int, something(getattr(node, "value", symbol = false), getattr(node, "bitpos", symbol = false))),
)

SpecBitmask(node::Node) =
    SpecBitmask(Symbol(node["name"]), StructVector(SpecBit.(findall("./enum[not(@alias)]", node))), parse(Int, getattr(node, "bitwidth", symbol = false, default = "32")))

function SpecFlag(node::Node)
    name = Symbol(findfirst("./name", node).content)
    typealias = Symbol(findfirst("./type", node).content)
    bitmask = if haskey(node, "requires")
        bitmask_name = Symbol(node["requires"])
        if bitmask_name ∈ disabled_specs
            nothing
        else
            spec_bitmasks[findfirst(==(bitmask_name), spec_bitmasks.name)]
        end
    else
        nothing
    end

    SpecFlag(name, typealias, bitmask)
end

function SpecHandle(node::Node)
    is_dispatchable = findfirst("./type", node).content == "VK_DEFINE_HANDLE"
    name = Symbol(findfirst("./name", node).content)
    SpecHandle(name, getattr(node, "parent"), is_dispatchable)
end

function SpecConstant(node::Node)
    name = Symbol(haskey(node, "name") ? node["name"] : findfirst("./name", node).content)
    value = if haskey(node, "offset")
        ext_value =
            parse(
                Int,
                something(
                    getattr(node, "extnumber", symbol = false),
                    getattr(node.parentnode.parentnode, "number", symbol = false),
                ),
            ) - 1
        offset = parse(Int, node["offset"])
        sign = (1 - 2 * (getattr(node, "dir", default = "", symbol = false) == "-"))
        sign * Int(1e9 + ext_value * 1e3 + offset)
    elseif haskey(node, "value")
        @match node["value"] begin
            "(~0U)" => :(typemax(UInt32))
            "(~0ULL)" => :(typemax(UInt64))
            "(~0U-1)" => :(typemax(UInt32) - 1)
            "(~0U-2)" => :(typemax(UInt32) - 2)
            "1000.0f" => :(1000.0f0)
            str::String && if contains(str, """)
            end => replace(str, """ => "")
            str => Meta.parse(str)
        end
    elseif haskey(node, "category")
        @match cat = node["category"] begin
            ("basetype" || "bitmask") => @match type = findfirst("./type", node) begin
                ::Nothing && if cat == "basetype"
                end => :Cvoid
                ::Node => translate_c_type(Symbol(type.content))
            end
        end
    else
        error("Unknown constant specification for node $node")
    end
    SpecConstant(name, value)
end

function SpecAlias(node::Node)
    name = Symbol(node["name"])
    alias = follow_alias(name)
    spec_names = getproperty.(spec_all_noalias, :name)
    alias_spec = spec_all_noalias[findfirst(==(alias), spec_names)]
    SpecAlias(name, alias_spec)
end

function CreateFunc(spec::SpecFunc)
    created_param = last(spec.params)
    handle = spec_handles[findfirst(==(follow_alias(innermost_type(created_param.type))), spec_handles.name)]
    create_info_params =
        getindex.(Ref(spec.params), findall(in(spec_create_info_structs.name), innermost_type.(spec.params.type)))
    @assert length(create_info_params) <= 1 "Found $(length(create_info_params)) create info types from the parameters of $spec:\n    $create_info_params"
    if length(create_info_params) == 0
        CreateFunc(spec, handle, nothing, nothing, false)
    else
        create_info_param = first(create_info_params)
        create_info_struct = spec_create_info_structs[findfirst(
            ==(innermost_type(create_info_param.type)),
            spec_create_info_structs.name,
        )]
        batch = is_arr(created_param)
        CreateFunc(spec, handle, create_info_struct, create_info_param, batch)
    end
end

function find_destroyed_param(spec::SpecFunc)
    idx = findlast(spec.params.is_externsync)
    @match idx begin
        ::Integer => spec.params[idx]
        ::Nothing => @match idx = findfirst(in(spec_handles.name), innermost_type.(spec.params.type)) begin
            ::Integer => spec.params[idx + 1]
            ::Nothing => error("Failed to retrieve the parameter to be destroyed:\n $spec")
        end
    end
end

function DestroyFunc(spec::SpecFunc)
    destroyed_param = find_destroyed_param(spec)
    handle = spec_handles[findfirst(==(innermost_type(destroyed_param.type)), spec_handles.name)]
    DestroyFunc(spec, handle, destroyed_param, is_arr(destroyed_param))
end

function queue_compatibility(node::Node)
    @when let _ = node, if haskey(node, "queues")
        end
        queue.(split(node["queues"], ','))
    end
end

externsync(node::Node) = haskey(node, "externsync") && node["externsync"] ≠ "false"

function len(node::Node)
    @match node begin
        _ && if haskey(node, "altlen")
        end => Meta.parse(node["altlen"])
        _ => @match getattr(node, "len", symbol = false) begin
            val::Nothing => nothing
            val => begin
                val_arr = filter(≠("null-terminated"), split(val, ','))
                if length(val_arr) > 1
                    if length(last(val_arr)) == 1 && isdigit(first(last(val_arr))) # array of pointers, length is unaffected
                        pop!(val_arr)
                    else
                        error("Failed to parse 'len' parameter '$val' for $(node.parentnode["name"]).")
                    end
                end
                isempty(val_arr) ? nothing : Symbol(first(val_arr))
            end
        end
    end
end

function arglen(node::Node; neighbor_type = "param")
    neighbors = findall("../$neighbor_type", node)
    arg_name = name(node)
    map(name, filter(x -> len(x) == arg_name, neighbors))
end

is_arr(spec::Union{SpecStructMember,SpecFuncParam}) = has_length(spec) && innermost_type(spec.type) ≠ :Cvoid
is_length(spec::Union{SpecStructMember,SpecFuncParam}) = !isempty(spec.arglen) && !endswith(string(spec.name), "Size")
is_size(spec::Union{SpecStructMember,SpecFuncParam}) = !isempty(spec.arglen) && endswith(string(spec.name), "Size")
has_length(spec::Union{SpecStructMember,SpecFuncParam}) = !isnothing(spec.len)
has_computable_length(spec::Union{SpecStructMember,SpecFuncParam}) =
    !spec.is_constant && spec.requirement == POINTER_REQUIRED && is_arr(spec)
is_data(spec::Union{SpecStructMember,SpecFuncParam}) = has_length(spec) && spec.type == :(Ptr{Cvoid})
is_version(spec::Union{SpecStructMember,SpecFuncParam}) =
    !isnothing(match(r"($v|V)ersion", string(spec.name))) && (
        follow_constant(spec.type) == :UInt32 ||
        is_ptr(spec.type) && !is_arr(spec) && !spec.is_constant && follow_constant(ptr_type(spec.type)) == :UInt32
    )

"""
Iterate through function or struct specification fields from a list of fields.
`list` is a sequence of fields to get through from `root`.
"""
struct FieldIterator
    root::Any
    list::Any
end

function Base.iterate(f::FieldIterator)
    spec = field(f.root, popfirst!(f.list))
    root = struct_by_name(innermost_type(spec.type))
    if isnothing(root)
        isempty(f.list) || error("Failed to retrieve a struct from $spec to continue the list $(f.list)")
        (spec, nothing)
    else
        (spec, FieldIterator(root, f.list))
    end
end

Base.iterate(_, f::FieldIterator) = iterate(f)
Base.iterate(f::FieldIterator, ::Nothing) = nothing
Base.length(f::FieldIterator) = length(f.list)

function field(spec, name)
    index = findfirst(==(name), children(spec).name)
    !isnothing(index) || error("Failed to retrieve field $name in $spec")
    children(spec)[index]
end

function field_struct(spec::Union{SpecStructMember,SpecFuncParam}, name)
    field_spec = field(parent_spec(spec), name)
    @match t = innermost_type(field_spec) begin
        GuardBy(in(spec_structs.name)) => struct_by_name(t)
        _ => field
    end
end

function length_chain(spec::Union{SpecStructMember,SpecFuncParam}, chain)
    parts = Symbol.(split(string(chain), "->"))
    collect(FieldIterator(parent_spec(spec), parts))
end

"""
Some specifications are disabled in the Vulkan headers (see https://github.com/KhronosGroup/Vulkan-Docs/issues/1225).
"""
const disabled_specs = map(x -> Symbol(x["name"]), findall("//extension[@supported = 'disabled']//*[@name]", xroot))

enabled_specs(specs) = filter(x -> x.name ∉ disabled_specs, specs)

"""
Specification constants, usually defined in C with #define.
"""
const spec_constants = enabled_specs(
    StructVector(
        SpecConstant.(
            vcat(
                findall("//enums[@name = 'API Constants']/*[@value and @name]", xroot),
                findall("//extension/require/enum[not(@extends) and not(@alias) and @value]", xroot),
                findall("/registry/types/type[@category = 'basetype' or @category = 'bitmask' and not(@alias)]", xroot),
            ),
        ),
    ),
)

"""
Specification enumerations, excluding bitmasks.
"""
const spec_enums = enabled_specs(StructVector(SpecEnum.(findall("//enums[@type = 'enum' and not(@alias)]", xroot))))

"""
Specification bitmask enumerations.
"""
const spec_bitmasks =
    enabled_specs(StructVector(SpecBitmask.(findall("//enums[@type = 'bitmask' and not(@alias)]", xroot))))

"""
Flag types.
"""
const spec_flags =
    enabled_specs(StructVector(SpecFlag.(findall("//type[@category = 'bitmask' and not(@alias)]", xroot))))

"""
Specification functions.
"""
const spec_funcs = enabled_specs(StructVector(SpecFunc.(findall("//command[not(@name)]", xroot))))

"""
Specification structures.
"""
const spec_structs =
    enabled_specs(StructVector(SpecStruct.(findall("//type[@category = 'struct' and not(@alias)]", xroot))))

"""
Specification for union types.
"""
const spec_unions =
    enabled_specs(StructVector(SpecUnion.(findall("//type[@category = 'union' and not(@alias)]", xroot))))

"""
Specification handle types.
"""
const spec_handles =
    enabled_specs(StructVector(SpecHandle.(findall("//type[@category = 'handle' and not(@alias)]", xroot))))

"""
All enumerations, regardless of them being bitmasks or regular enum values.
"""
const spec_all_semantic_enums = [spec_enums..., spec_bitmasks...]

# extend all types with core and extension values
let nodes = findall("//*[@extends and not(@alias)]", xroot)
    database = vcat(spec_all_semantic_enums, spec_structs)
    foreach(nodes) do node
        spec = database[findfirst(x -> x.name == Symbol(node["extends"]), database)]
        @switch spec begin
            @case ::SpecBitmask
            push!(spec.bits, SpecBit(node))
            @case ::SpecEnum
            push!(spec.enums, SpecConstant(node))
        end
    end
end

"""
All specifications except aliases.
"""
const spec_all_noalias = [
    spec_funcs...,
    spec_structs...,
    spec_handles...,
    spec_constants...,
    spec_enums...,
    (spec_enums.enums...)...,
    spec_bitmasks...,
    (spec_bitmasks.bits...)...,
]

"""
All specification aliases.
"""
const spec_aliases = StructVector(
    SpecAlias.(
        filter(
            x -> (par = x.parentnode.parentnode;
            par.name ≠ "extension" || getattr(par, "supported") ≠ :disabled),
            findall("//*[@alias and @name]", xroot),
        ),
    ),
)
const spec_func_params = collect(Iterators.flatten(spec_funcs.params))
const spec_struct_members = collect(Iterators.flatten(spec_structs.members))
const spec_create_info_structs = filter(x -> x.type ∈ [CREATE_INFO, ALLOCATE_INFO], spec_structs)

function spec_by_field(specs, field, value)
    specs[findall(==(value), getproperty(specs, field))]
end

function spec_by_name(specs, name)
    specs = spec_by_field(specs, :name, name)
    if !isempty(specs)
        length(specs) == 1 || error("Non-uniquely identified spec '$name': $specs")
        first(specs)
    else
        nothing
    end
end

func_by_name(name) = spec_by_name(spec_funcs, name)

struct_by_name(name) = spec_by_name(spec_structs, name)

union_by_name(name) = spec_by_name(spec_unions, name)

handle_by_name(name) = spec_by_name(spec_handles, name)

bitmask_by_name(name) = spec_by_name(spec_bitmasks, name)

flag_by_name(name) = spec_by_name(spec_flags, name)

enum_by_name(name) = spec_by_name(spec_enums, name)

constant_by_name(name) = spec_by_name(spec_constants, name)

Base.parent(spec::SpecFuncParam) = spec.func
Base.parent(spec::SpecStructMember) = spec.parent
Base.parent(spec::SpecHandle) = spec.parent

parent_spec(spec::SpecFuncParam) = func_by_name(parent(spec))
parent_spec(spec::SpecStructMember) = struct_by_name(parent(spec))
parent_spec(spec::SpecHandle) = handle_by_name(parent(spec))

function len(spec::Union{SpecFuncParam,SpecStructMember})
    params = children(parent_spec(spec))
    params[findfirst(x -> x.name == spec.len, params)]
end

function arglen(spec::Union{SpecFuncParam,SpecStructMember})
    params = children(parent_spec(spec))
    params[findall(x -> x.name ∈ spec.arglen, params)]
end

const spec_create_funcs = StructVector(CreateFunc.(filter(x -> x.type ∈ [CREATE, ALLOCATE], spec_funcs)))
const spec_destroy_funcs = StructVector(DestroyFunc.(filter(x -> x.type ∈ [DESTROY, FREE], spec_funcs)))

function wrappable_constructors(handle::SpecHandle)::Vector{CreateFunc}
    # don't wrap VkSurfaceKHR, almost all signatures conflict with one another with create info parameters exposed
    handle.name == :VkSurfaceKHR && return []
    all_constructors_nobatch = filter(x -> x.handle == handle && !x.batch, spec_create_funcs)
    if length(unique(all_constructors_nobatch.create_info_struct)) == length(all_constructors_nobatch)
        all_constructors_nobatch
    else
        []
    end
end

const spec_handles_with_wrappable_constructors = filter(x -> !isempty(wrappable_constructors(x)), spec_handles)

is_destructible(spec::SpecHandle) = spec ∈ spec_destroy_funcs.handle

default(::SpecHandle) = :C_NULL
default(spec::Union{SpecStructMember,SpecFuncParam}) = @match spec.requirement begin
    if spec.type ∈ spec_handles.name
    end => default(handle_by_name(spec.type))
    &POINTER_OPTIONAL || &POINTER_REQUIRED || if is_ptr(spec.type) || spec.type == :Cstring
    end => :C_NULL
    &OPTIONAL || &REQUIRED => 0
end

create_func(func::SpecFunc) = first(spec_by_field(spec_create_funcs, :func, func))
create_func(name) = first(spec_by_field(spec_create_funcs, :func, func_by_name(name)))

create_funcs(handle::SpecHandle) = spec_by_field(spec_create_funcs, :handle, handle)

function create_func_no_batch(handle::SpecHandle)
    fs = create_funcs(handle)
    fs[findfirst(x -> x.handle == handle && !x.batch, fs)]
end

destroy_func(func::SpecFunc) = first(spec_by_field(spec_destroy_funcs, :func, func))
destroy_func(name) = first(spec_by_field(spec_destroy_funcs, :func, func_by_name(name)))

destroy_funcs(handle::SpecHandle) = spec_by_field(spec_destroy_funcs, :handle, handle)

function follow_constant(spec::SpecConstant)
    @match val = spec.value begin
        ::Symbol && if val ∈ spec_constants.name
        end => follow_constant(constant_by_name(val))
        val => val
    end
end

function follow_constant(name)
    constant = constant_by_name(name)
    isnothing(constant) ? name : follow_constant(constant)
end


text:
module StoppingTutorials

deps_jl = joinpath(@__DIR__,"..", "deps", "deps.jl")

if !isfile(deps_jl)
  s = """
  Package StoppingTutorials not installed properly.
  Run Pkg.build(\"StoppingTutorials\"), restart Julia and try again
  """
  error(s)
end

include(deps_jl)

end


text:
"""
    AbstractBernsteinForm{T}

Abstract supertype for all Bernstein forms types.

### Notes

This type is parametric in the numeric type `T`.

```jldoctest
julia> subtypes(AbstractBernsteinForm)

```
"""
abstract type AbstractBernsteinForm{T} end

"""
    numvars(::AbstractBernsteinForm)

Return the number of variables of the given Bernstein form.

### Input

### Output
"""
function numvars(::AbstractBernsteinForm) end

"""
    polynomial(::AbstractBernsteinForm)

Return the polynomial represented by the given Bernstein form.

### Input

### Output

"""
function polynomial(::AbstractBernsteinForm) end

"""
    domain(::AbstractBernsteinForm)
"""
function domain(::AbstractBernsteinForm) end

"""
    enclose(::AbstractBernsteinForm, ::IntervalOrIntervalBox)

### Input

### Output
"""
function enclose(::AbstractBernsteinForm, ::IntervalOrIntervalBox) end

"""
    getcoeff(::AbstractBernsteinForm, ids::NTuple{D, Int}) where {D}

Return the coefficient of the Bernstein expansion corresponding to the tuple `ids`.

### Input

### Output
"""
function getcoeff(::AbstractBernsteinForm, ids::NTuple{D, Int}) where {D} end

max_stars_repo_name: JuliaReach/BernsteinExpansions

max_issues_repo_name: JuliaReach/BernsteinExpansions

max_forks_repo_name: JuliaReach/BernsteinExpansions


text:
struct DeivVecTag end

# J(f(x))*v
function auto_jacvec!(
    dy,
    f,
    x,
    v,
    cache1 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x))))),
    cache2 = similar(cache1),
)
    cache1 .= Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x)))))
    f(cache2, cache1)
    vecdy = _vec(dy)
    vecdy .= partials.(_vec(cache2), 1)
end

_vec(v) = vec(v)
_vec(v::AbstractVector) = v

function auto_jacvec(f, x, v)
    vv = reshape(v, axes(x))
    y = ForwardDiff.Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(vv)))
    vec(partials.(vec(f(y)), 1))
end

function num_jacvec!(
    dy,
    f,
    x,
    v,
    cache1 = similar(v),
    cache2 = similar(v);
    compute_f0 = true,
)
    vv = reshape(v, axes(x))
    compute_f0 && (f(cache1, x))
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ * vv
    f(cache2, x)
    @. x -= ϵ * vv
    vecdy = _vec(dy)
    veccache1 = _vec(cache1)
    veccache2 = _vec(cache2)
    @. vecdy = (veccache2 - veccache1) / ϵ
end

function num_jacvec(f, x, v, f0 = nothing)
    vv = reshape(v, axes(x))
    f0 === nothing ? _f0 = f(x) : _f0 = f0
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(minimum(x)))
    vec((f(x .+ ϵ .* vv) .- _f0) ./ ϵ)
end

function num_hesvec!(
    dy,
    f,
    x,
    v,
    cache1 = similar(v),
    cache2 = similar(v),
    cache3 = similar(v),
)
    cache = FiniteDiff.GradientCache(v[1], cache1, Val{:central})
    g = let f = f, cache = cache
        (dx, x) -> FiniteDiff.finite_difference_gradient!(dx, f, x, cache)
    end
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ * v
    g(cache2, x)
    @. x -= 2ϵ * v
    g(cache3, x)
    @. dy = (cache2 - cache3) / (2ϵ)
end

function num_hesvec(f, x, v)
    g = (x) -> FiniteDiff.finite_difference_gradient(f, x)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    x += ϵ * v
    gxp = g(x)
    x -= 2ϵ * v
    gxm = g(x)
    (gxp - gxm) / (2ϵ)
end

function numauto_hesvec!(
    dy,
    f,
    x,
    v,
    cache = ForwardDiff.GradientConfig(f, v),
    cache1 = similar(v),
    cache2 = similar(v),
)
    g = let f = f, x = x, cache = cache
        g = (dx, x) -> ForwardDiff.gradient!(dx, f, x, cache)
    end
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ * v
    g(cache1, x)
    @. x -= 2ϵ * v
    g(cache2, x)
    @. dy = (cache1 - cache2) / (2ϵ)
end

function numauto_hesvec(f, x, v)
    g = (x) -> ForwardDiff.gradient(f, x)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    x += ϵ * v
    gxp = g(x)
    x -= 2ϵ * v
    gxm = g(x)
    (gxp - gxm) / (2ϵ)
end

function autonum_hesvec!(
    dy,
    f,
    x,
    v,
    cache1 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x))))),
    cache2 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x))))),
)
    cache = FiniteDiff.GradientCache(v[1], cache1, Val{:central})
    g = (dx, x) -> FiniteDiff.finite_difference_gradient!(dx, f, x, cache)
    cache1 .= Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x)))))
    g(cache2, cache1)
    dy .= partials.(cache2, 1)
end

function autonum_hesvec(f, x, v)
    g = (x) -> FiniteDiff.finite_difference_gradient(f, x)
    partials.(g(Dual{DeivVecTag}.(x, v)), 1)
end

function num_hesvecgrad!(dy, g, x, v, cache2 = similar(v), cache3 = similar(v))
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ * v
    g(cache2, x)
    @. x -= 2ϵ * v
    g(cache3, x)
    @. dy = (cache2 - cache3) / (2ϵ)
end

function num_hesvecgrad(g, x, v)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    x += ϵ * v
    gxp = g(x)
    x -= 2ϵ * v
    gxm = g(x)
    (gxp - gxm) / (2ϵ)
end

function auto_hesvecgrad!(
    dy,
    g,
    x,
    v,
    cache2 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x))))),
    cache3 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x))))),
)
    cache2 .= Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x)))))
    g(cache3, cache2)
    dy .= partials.(cache3, 1)
end

function auto_hesvecgrad(g, x, v)
    y = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x)))))
    partials.(g(y), 1)
end

### Operator Forms

struct JacVec{F,T1,T2,xType}
    f::F
    cache1::T1
    cache2::T2
    x::xType
    autodiff::Bool
end

function JacVec(f, x::AbstractArray; autodiff = true)
    if autodiff
        cache1 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(x)))
        cache2 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(x)))
    else
        cache1 = similar(x)
        cache2 = similar(x)
    end
    JacVec(f, cache1, cache2, x, autodiff)
end

Base.eltype(L::JacVec) = eltype(L.x)
Base.size(L::JacVec) = (length(L.cache1), length(L.cache1))
Base.size(L::JacVec, i::Int) = length(L.cache1)
Base.:*(L::JacVec, v::AbstractVector) =
    L.autodiff ? auto_jacvec(_x -> L.f(_x), L.x, v) :
    num_jacvec(_x -> L.f(_x), L.x, v)

function LinearAlgebra.mul!(dy::AbstractVector, L::JacVec, v::AbstractVector)
    if L.autodiff
        auto_jacvec!(dy, (_y, _x) -> L.f(_y, _x), L.x, v, L.cache1, L.cache2)
    else
        num_jacvec!(dy, (_y, _x) -> L.f(_y, _x), L.x, v, L.cache1, L.cache2)
    end
end

struct HesVec{F,T1,T2,xType}
    f::F
    cache1::T1
    cache2::T2
    cache3::T2
    x::xType
    autodiff::Bool
end

function HesVec(f, x::AbstractArray; autodiff = true)
    if autodiff
        cache1 = ForwardDiff.GradientConfig(f, x)
        cache2 = similar(x)
        cache3 = similar(x)
    else
        cache1 = similar(x)
        cache2 = similar(x)
        cache3 = similar(x)
    end
    HesVec(f, cache1, cache2, cache3, x, autodiff)
end

Base.size(L::HesVec) = (length(L.cache2), length(L.cache2))
Base.size(L::HesVec, i::Int) = length(L.cache2)
Base.:*(L::HesVec, v::AbstractVector) =
    L.autodiff ? numauto_hesvec(L.f, L.x, v) : num_hesvec(L.f, L.x, v)

function LinearAlgebra.mul!(dy::AbstractVector, L::HesVec, v::AbstractVector)
    if L.autodiff
        numauto_hesvec!(dy, L.f, L.x, v, L.cache1, L.cache2, L.cache3)
    else
        num_hesvec!(dy, L.f, L.x, v, L.cache1, L.cache2, L.cache3)
    end
end

struct HesVecGrad{G,T1,T2,uType}
    g::G
    cache1::T1
    cache2::T2
    x::uType
    autodiff::Bool
end

function HesVecGrad(g, x::AbstractArray; autodiff = false)
    if autodiff
        cache1 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(x)))
        cache2 = Dual{typeof(ForwardDiff.Tag(DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(x)))
    else
        cache1 = similar(x)
        cache2 = similar(x)
    end
    HesVecGrad(g, cache1, cache2, x, autodiff)
end

Base.size(L::HesVecGrad) = (length(L.cache2), length(L.cache2))
Base.size(L::HesVecGrad, i::Int) = length(L.cache2)
Base.:*(L::HesVecGrad, v::AbstractVector) =
    L.autodiff ? auto_hesvecgrad(L.g, L.x, v) : num_hesvecgrad(L.g, L.x, v)

function LinearAlgebra.mul!(
    dy::AbstractVector,
    L::HesVecGrad,
    v::AbstractVector,
)
    if L.autodiff
        auto_hesvecgrad!(dy, L.g, L.x, v, L.cache1, L.cache2)
    else
        num_hesvecgrad!(dy, L.g, L.x, v, L.cache1, L.cache2)
    end
end


text:
# fix eta k cchoice
# fix quadstep?

struct HZAW{T} <: LineSearch
  decrease::T
  curvature::T
  θ::T # if 0.5 -> bisection in U3
  γ::T
end
function HZAW(; decrease=0.1, curvature=0.9, theta=0.5, gamma=2/3)
  if !(0 < decrease ≤ curvature)
    println("Decrease constant must be positive and smaller than the curvature condition.")
  end
  if !(curvature < 1)
    println("Curvature constant must be smaller than one.")
  end
  HZAW(decrease, curvature, theta, gamma)
end
function find_steplength(hzl::HZAW, φ, c)
  # c = initial(k) but this is done outisde
  T = typeof(φ.φ0)
  ϵk = 1e-6
  δ = T(hzl.decrease)
  σ = T(hzl.curvature)
  φ0, dφ0 = φ.φ0, φ.dφ0
  φc, dφc = φ0, dφ0
  # initial convergence
  # Wolfe conditions
  if δ*dφ0 ≥ (φc-φ0)/c && dφc ≥ σ*dφ0
    return c, φc, true
  end
  # Approximate Wolfe conditions
  if (2*δ-1)*dφ0 ≥ dφc ≥ σ*dφ0 && φc ≤ φ0 + ϵk
    return c, φc, true
  end

  a0, b0 = bracket(hzl, c, φ)
  j = 0
  aj, bj  = a0, b0
  while j < 100
    a, b = secant²(hzl, φ, aj, bj)
    if b-a>hzl.γ*(bj-aj)
      c = (a+b)/2
      φc, dφc = φ(c, true)
      a, b = update(hzl, a, b, c, φ, φc, dφc, ϵk)
    end
   aj, bj = a, b
    j += 1
    if δ*dφ0 ≥ (φc-φ0)/c && dφc ≥ σ*dφ0
      return c, φc, true
    end
    if (2*δ-1)*dφ0 ≥ dφc ≥ σ*dφ0 && φc ≤ φ0 + ϵk
      return c, φc, true
    end
  end
  return c, φc, true
end

"""
   _U3

Used to take step U3 of the updating procedure [HZ, p.123]. The other steps
are in update, but this step is separated out to be able to use it in
step B2 of bracket.
"""
function _U3(hzl::HZAW, φ, a::T, b::T, c::T, φ0, φc, dφc, ϵk) where T
  _a, _b = a, c
  # a)
  searching = true
  j = 1
  while searching && j < 50
    d = (1-hzl.θ)*_a+ hzl.θ*_b
    φd, dφd = φ(d, true)
    if dφd ≥ T(0)
      _b = d
      searching = false
    else
      if dφd < T(0) && φd < φ0 + ϵk
        _a = d
      else # φ(d) ≥ φ0 + ϵk
        _b = d
      end
    end
    j += 1
  end
  _a, _b
end

function update(hzl::HZ,
                 a::T, b::T, c::T,
                 φ, φc, dφc,
                 ϵk) where {HZ<:HZAW, T}

  φ0 = φ.φ0
  #== U0 ==#
  if c < a || c > b
    return a, b, (a=false, b=false)
  end
  #== U1 ==#
  if dφc ≥ T(0)
    return a, c, (a=false, b=true)
  end
  #== U2 ==#
  if dφc < T(0) && φc ≤ φ0 + ϵk
    return c, b, (a=true, b=false)
  end
  #== U3 ==#
  a, b = _U3(hzl, φ, a, b, c, φ0, φc, dφc, ϵk)
  return a, b, (a=a==c, b=b==c)
end

"""
  bracket

Find an interval satisfying the opposite slope condition [OSC] starting from
[0, c].
"""
function bracket(hzl::HZAW, c::T, φ) where T
  # Note, we know that dφ(0) < 0 since we're accepted that the current step is in a
  # direction of descent.
  φ0 = φ.φ0
  ρ = T(5)
  ϵk = T(1)/10^6
  #== B0 ==#
  cj = c
  φc, dφc = (cj, true)
  # we only want to store a number, so we don't store all iterates
  ci, φi = T(0), φ0

  maxj = 100
  for j = 0:maxj
    #== B1: φ is increasing at c, set b to cj as this is an upper bound ==#
    if dφc ≥ T(0)
      a, b = ci, cj
      return a, b
    end
    #== B2 : φ is decreasing at cj but function value is sufficiently larger than
    # φ0, use U3 to update. ==#
    if dφc < T(0) && φc > φ0 + ϵk
      a, b = _U3(hzl, a, b, c, φ, φ0, φc, dφc, ϵk)
      return a, b
    end
    #== B3 ==#
    ci_candidate = cj
    φci_candidate = φc
    cj = ρ*cj
    φc, dφc = φ(cj, true)
    # update ci instead of keeping all c's
    if φci_candidate ≤ φ0 + ϵk
      ci = ci_candidate
      φci = φci_candidate
    end
  end
end

function secant(hzl::HZAW, a, dφa, b, dφb)
  (a*dφb + b*dφa)/(dφb-dφa)
end
function secant²(hzl::HZAW, φ, a, b)
  ϵk = 1e-6
  #== S1 ==#
  φa, dφa = φ(a, true)
  φb, dφb = φ(b, true)
  c = secant(hzl, a, dφa, b, dφb)
  φc, dφc = φ(c, true)
  A, B, updates = update(hzl, a, b, c, φ, φc, dφc, ϵk)
  if updates.b
    #== S2: c is the upper bound ==#
    φB, dφB = φ(B, true)
    _c = secant(hzl, b, dφb, B, dφB)
  elseif updates.a
    #== S3: c is the lower bound ==#
    φA, dφA = φ(A, true)
    _c = secant(hzl, a, dφa, A, dφA)
  end

  if any(updates)
    #== S4.if: c was upper or lower bound ==#
    φ_c, dφ_c = φ(_c, true)
    _a, _b = update(hzl, A, B, _c, φ, φ_c, dφ_c, ϵk)
    return _a, _b
  else
    #== S4.otherwise: c was neither ==#
    return A, B
  end
end


text:
module EasyHamiltonian
    import ITensors.AutoMPO

    include("expression.jl")
    include("hamiltonian_type.jl")
    include("translate.jl")

    # Model zoo
    include("./hamiltonian_zoo/SpinfulFermionicHamiltonian.jl")
end


text:
@testset "point_light_attributes" begin
    intensity = RayTracer.color(1,1,1)
    position = RayTracer.point(0,0,0)
    light = RayTracer.point_light(position, intensity)
    @test light.position ≈ position
    @test light.intensity ≈ intensity
end


text:
#### dims2indices ####

"""
    dims2indices(dim::Dimension, I) => NTuple{Union{Colon,AbstractArray,Int}}

Convert a `Dimension` or `Selector` `I` to indices of `Int`, `AbstractArray` or `Colon`.
"""
@inline dims2indices(dim::Dimension, I::StandardIndices) = I
@inline dims2indices(dim::Dimension, I) = _dims2indices(dim, I)

@inline dims2indices(x, I) = dims2indices(dims(x), I)
@inline dims2indices(::Nothing, I) = _dimsnotdefinederror()
@inline dims2indices(dims::DimTuple, I) = dims2indices(dims, (I,))
# Standard array indices are simply returned
@inline dims2indices(dims::DimTuple, I::Tuple{Vararg{<:StandardIndices}}) = I
@inline dims2indices(dims::DimTuple, I::Tuple{<:CartesianIndex}) = I
@inline dims2indices(dims::DimTuple, sel::Tuple) =
    LookupArrays.selectindices(lookup(dims), sel)
@inline dims2indices(dims::DimTuple, ::Tuple{}) = ()
# Otherwise attempt to convert dims to indices
@inline function dims2indices(dims::DimTuple, I::DimTuple)
    extradims = otherdims(I, dims)
    length(extradims) > 0 && _warnextradims(extradims)
    _dims2indices(lookup(dims), dims, sortdims(I, dims))
end

# Handle tuples with @generated
@inline _dims2indices(::Tuple{}, dims::Tuple{}, ::Tuple{}) = ()
@generated function _dims2indices(lookups::Tuple, dims::Tuple, I::Tuple)
    # We separate out Aligned and Unaligned lookups as
    # Unaligned must be selected in groups e.g. X and Y together.
    unalligned = Expr(:tuple)
    uaI = Expr(:tuple)
    alligned = Expr(:tuple)
    dimmerge = Expr(:tuple)
    a_count = ua_count = 0
    for (i, lkup) in enumerate(lookups.parameters)
        if lkup <: Unaligned
            ua_count += 1
            push!(unalligned.args, :(dims[$i]))
            push!(uaI.args, :(I[$i]))
            push!(dimmerge.args, :(uadims[$ua_count]))
        else
            a_count += 1
            push!(alligned.args, :(_dims2indices(dims[$i], I[$i])))
            # Update the merged tuple
            push!(dimmerge.args, :(adims[$a_count]))
        end
    end

    if length(unalligned.args) > 1
        # Output the `dimmerge` tuple, which will
        # combine uadimsand adims in the right order
        quote
             adims = $alligned
             # Unaligned dims have to be run together as a set
             uadims = unalligned_dims2indices($unalligned, map(_unwrapdim, $uaI))
             $dimmerge
        end
    else
        alligned
    end
end

@inline function unalligned_dims2indices(dims::DimTuple, sel::Tuple)
    map(sel) do s
        s isa Selector && throw(ArgumentError("Unalligned dims: use selectors for all $(join(map(string ∘ dim2key, dims), ", ")) dims, or none of them"))
        isnothing(s) ? Colon() : s
    end
end
@inline function unalligned_dims2indices(dims::DimTuple, sel::Tuple{<:Selector,Vararg{<:Selector}})
    LookupArrays.select_unalligned_indices(lookup(dims), sel)
end

_unwrapdim(dim::Dimension) = val(dim)
_unwrapdim(x) = x

# Single dim methods
# A Dimension type always means Colon(), as if it was constructed with the default value.
@inline _dims2indices(dim::Dimension, ::Type{<:Dimension}) = Colon()
# Nothing means nothing was passed for this dimension
@inline _dims2indices(dim::Dimension, ::Nothing) = Colon()
# Simply unwrap dimensions
@inline _dims2indices(dim::Dimension, seldim::Dimension) =
    LookupArrays.selectindices(val(dim), val(seldim))