text:
export  IntegratedFunction,
        active_dim,
        getbasis,
        grad_xd,
        grad_coeff_grad_xd,
        hess_coeff_grad_xd,
        integrate_xd,
        evaluate!,
        evaluate,
        grad_x!,
        grad_x,
        reduced_grad_x!,
        reduced_grad_x,
        hess_x!,
        hess_x,
        reduced_hess_x!,
        reduced_hess_x,
        grad_coeff_integrate_xd,
        hess_coeff_integrate_xd,
        grad_coeff,
        hess_coeff,
        evalgrad_coeff!

"""
$(TYPEDEF)

An immutable structure to hold the integrated function
with `g` a `Rectifier` and `f` an `ExpandedFunction`.

## Fields

$(TYPEDFIELDS)
"""
struct IntegratedFunction
    m::Int64
    Nψ::Int64
    Nx::Int64
    g::Rectifier
    f::ExpandedFunction
end

function IntegratedFunction(f::ExpandedFunction; rectifier::String="softplus")
    return IntegratedFunction(f.m, f.Nψ, f.Nx, Rectifier(rectifier), f)
end

"""
$(TYPEDSIGNATURES)

Extracts the kind of basis of the integrated function `R`.
"""
getbasis(R::IntegratedFunction) = getbasis(R.f)

"""
$(TYPEDSIGNATURES)

Extracts the active dimensions of the integrated function `R`.
"""
active_dim(R::IntegratedFunction) = R.f.dim

"""
$(TYPEDSIGNATURES)

Computes ∫_0^xₖ g(∂ₖ f(x_{1:k-1}, t))dt for the ensemble matrix `X`.
"""
function integrate_xd(R::IntegratedFunction, X::Array{Float64,2})
    NxX, Ne = size(X)
    ψoff = evaluate_offdiagbasis(R.f, X)
    Nx = R.Nx
    xk = deepcopy(X[Nx, :])
    cache = zeros(Ne)
    # ψoffdxdψ = zeros(Ne, Nψ)

    function integrand!(v::Vector{Float64}, t::Float64)
        # ψoffdxdψ .= repeated_grad_xk_basis(R.f,  0.5*(t+1)*xk)
        # ψoffdxdψ .*= ψoff
        # map!(R.g, v, ψoffdxdψ*R.f.coeff)
        evaluate!(v, R.g, (repeated_grad_xk_basis(R.f,  t*xk) .* ψoff)*R.f.coeff)
    end

    return xk .* quadgk!(integrand!, cache, 0.0, 1.0)[1]
end

# Compute g(∂ₖf(x_{1:k}))
"""
$(TYPEDSIGNATURES)

Computes g(∂ₖ f(x_{1:k})) for the ensemble matrix `X`.
"""
function grad_xd(R::IntegratedFunction, X)
    dψ = grad_xd(R.f, X)
    evaluate!(dψ, R.g, dψ)
    # @show size(dψ)
    # gdψ = R.g.(dψ)
    return dψ
end

# Compute ∂_c( g(∂ₖf(x_{1:k}) ) ) = ∂_c∂ₖf(x_{1:k}) × g′(∂ₖf(x_{1:k}))
"""
$(TYPEDSIGNATURES)

Computes ∂_c( g(∂ₖf(x_{1:k})) ) = ∂_c∂ₖf(x_{1:k}) × g′(∂ₖf(x_{1:k})) for the ensemble matrix `X`.
"""
function grad_coeff_grad_xd(R::IntegratedFunction, X::Array{Float64,2})
    dψ = grad_xd(R.f, X)
    dcdψ = grad_coeff_grad_xd(R.f, X)
    return vgrad_x(R.g, dψ) .* dcdψ
end

# Compute ∂²_c( g(∂ₖf(x_{1:k}) ) ) = ∂²_c∂ₖf(x_{1:k}) × g′(∂ₖf(x_{1:k})) + ∂_c∂ₖf(x_{1:k}) × ∂_c∂ₖf(x_{1:k}) g″(∂ₖf(x_{1:k}))
"""
$(TYPEDSIGNATURES)

Computes ∂²_c( g(∂ₖf(x_{1:k}) ) ) = ∂²_c∂ₖf(x_{1:k}) × g′(∂ₖf(x_{1:k})) + ∂_c∂ₖf(x_{1:k}) × ∂_c∂ₖf(x_{1:k}) g″(∂ₖf(x_{1:k})) for the ensemble matrix `X`.
"""
function hess_coeff_grad_xd(R::IntegratedFunction, X::Array{Float64,2})
    # The second term can be dropped for improve performance
    Nψ = R.Nψ
    NxX, Ne = size(X)
    dψ    = grad_xd(R.f, X)
    dcdψ  = grad_coeff_grad_xd(R.f, X)
    # d2cdψ = hess_coeff_grad_xd(R.f, X)

    dcdψouter = zeros(Ne, Nψ, Nψ)
    @inbounds for i=1:Nψ
        for j=1:Nψ
            dcdψouter[:,i,j] = dcdψ[:,i] .* dcdψ[:, j]
        end
    end
    return vhess_x(R.g, dψ) .* dcdψouter # + dψ .* d2cdψ
end

## integrate_xd
# Compute ∫_0^{x_k} g(∂ₖf(x_{1:k-1},t)) dt

# function repeated_grad_xk_basis(f::ExpandedFunction, x, idx::Array{Int64,2})
#     # Compute the k=th order deriviative of an expanded function along the direction grad_dim
#     N = size(x,1)
#     Nx = f.Nx
#     # ∂ᵏf/∂x_{grad_dim} = ψ
#     k = 1
#     grad_dim = Nx
#     dims = Nx
#
#     midxj = idx[:,Nx]
#     maxj = maximum(midxj)
#     #   Compute the kth derivative along grad_dim
#     dkψj = vander(f.B.B, maxj, k, x)
#     return dkψj[:, midxj .+ 1]
# end

# function evaluate!(out::Array{Float64,1}, R::IntegratedFunction{m, Nψ, Nx}, X::Array{Float64,2}) where {m, Nψ, Nx}
"""
$(TYPEDSIGNATURES)

Evaluates in-place the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
"""
function evaluate!(out, R::IntegratedFunction, X)
    NxX, Ne = size(X)
    Nx = R.Nx
    @assert NxX == Nx "Wrong dimension of the sample X"
    ψoff = evaluate_offdiagbasis(R.f, X)
    ψdiag = repeated_evaluate_basis(R.f, zeros(Ne))
    xlast = view(X,NxX,:)
    cache = zeros(Ne)

    @assert size(out,1) == Ne

    function integrand!(v::Vector{Float64}, t::Float64)
        v .= (repeated_grad_xk_basis(R.f,  t*xlast) .* ψoff) * R.f.coeff
        evaluate!(v, R.g, v)

        # v .= R.g((repeated_grad_xk_basis(R.f,  t*xk) .* ψoff)*R.f.coeff)
    end

     out .= (ψoff .* ψdiag)*R.f.coeff + xlast .* quadgk!(integrand!, cache, 0.0, 1.0)[1]

     return out
 end

"""
$(TYPEDSIGNATURES)

Evaluates the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
"""
evaluate(R::IntegratedFunction, X::Array{Float64,2}) = evaluate!(zeros(size(X,2)), R, X)
evaluate(R::IntegratedFunction, x::Array{Float64,1}) = evaluate!(zeros(1), R, reshape(x,(size(x,1), 1)))[1]

## Gradient of an IntegratedFunction

# function grad_x!(out, R::IntegratedFunction, X)
#     NxX, Ne = size(X)
#     Nx = R.Nx
#     @assert NxX == Nx "Wrong dimension of the sample"
#     @assert size(out) == (Ne, Nx) "Dimensions of the output and the samples don't match"
#
#     x0 = zeros(Ne)
#     xlast = view(X,Nx,:)
#     ψk0  = repeated_evaluate_basis(R.f, x0)
#     ψoff = evaluate_offdiagbasis(R.f, X)
#     dxkψ = zero(ψk0)
#
#     coeff = R.f.coeff
#
#     # Cache for the integration
#     cache = zeros((Nx-1)*Ne)
#     cacheg = zeros(Ne)
#
#     # Compute the basis for each component
#     ψbasis = zeros(Ne, R.Nψ, Nx-1)
#     @inbounds for i=1:Nx-1
#         ψbasis_i = view(ψbasis, :, :, i)
#         ψbasis_i .= evaluate_basis(R.f, X, [i], R.f.idx)
#     end
#
#     # Compute ψ1 ⊗ ψ2 ⊗ ψi′ ⊗ … ⊗ ψk-1
#     dxψbasis = zero(ψbasis)
#     fill!(dxψbasis, 1.0)
#     @inbounds for i=1:Nx-1
#         for j=1:Nx-1
#             if i==j
#             dxψbasis[:,:,i] .*= grad_xk_basis(R.f, X, 1, j, j, R.f.idx)
#             else
#             dxψbasis[:,:,i] .*= ψbasis[:,:,j]
#             end
#         end
#     end
#
#
#     # This integrand computes ψ1 ⊗ … ⊗ ψi′ ⊗ … ⊗ ψk′(t) × c g′(ψ1 ⊗ … ⊗ ψk′(t)×c)
#     function integrand!(v::Vector{Float64}, t::Float64)
#     dxkψ .= repeated_grad_xk_basis(R.f,  t*xlast)
#     @avx @. cacheg = (dxkψ * ψoff) *ˡ coeff
#     grad_x!(cacheg, R.g, cacheg)
#
#         @inbounds for i=1:Nx-1
#             vi = view(v, (i-1)*Ne+1:i*Ne)
#             dxψbasisi = view(dxψbasis,:,:,i)
#             vi .= ((dxψbasisi .* dxkψ) * coeff) .* cacheg
#         end
#     end
#     quadgk!(integrand!, cache, 0.0, 1.0; rtol = 1e-3)
#
#     # Multiply integral by xlast (change of variable in the integration)
#     @inbounds for i=1:Nx-1
#         @. cache[(i-1)*Ne+1:i*Ne] *= xlast
#     end
#
#     # Add ψ1 ⊗ … ⊗ ψi′ ⊗ … ⊗ ψk(0) × c
#     @inbounds for i=1:Nx-1
#         dxψbasisi = view(dxψbasis,:,:,i)
#         view(out,:,i) .+= (dxψbasisi .* ψk0)*coeff + view(cache,(i-1)*Ne+1:i*Ne)
#     end
#
#     # Fill the last column
#     lastcol = view(out, :, Nx)
#     evaluate!(lastcol, R.g, (repeated_grad_xk_basis(R.f,  xlast) .* ψoff)*coeff)
#     return out
# end

"""
$(TYPEDSIGNATURES)

Evaluates in-place the gradient (with respect to the state) the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates ∂ₓ f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
"""
function grad_x!(out, R::IntegratedFunction, X)
    NxX, Ne = size(X)
    Nx = R.Nx
    @assert NxX == Nx "Wrong dimension of the sample"
    @assert size(out) == (Ne, Nx) "Dimensions of the output and the samples don't match"

    x0 = zeros(Ne)
    xlast = view(X,Nx,:)
    ψk0  = repeated_evaluate_basis(R.f, x0)
    ψoff = evaluate_offdiagbasis(R.f, X)
    dxkψ = zero(ψk0)

    coeff = R.f.coeff

    # Define active and off diagonal active dimension
    dim = active_dim(R)
    dimoff = dim[dim .< Nx]

    # Cache for the integration
    cache = zeros(length(dimoff)*Ne)
    cacheg = zeros(Ne)

    # Compute the basis for each component
    ψbasis = zeros(Ne, R.Nψ, length(dimoff))
    @inbounds for (i, dimi) in enumerate(dimoff)
        ψbasis_i = view(ψbasis, :, :, i)
        ψbasis_i .= evaluate_basis(R.f, X, [dimi], R.f.idx)
    end

    # Compute ψ1 ⊗ ψ2 ⊗ ψi′ ⊗ … ⊗ ψk-1
    dxψbasis = zero(ψbasis)
    fill!(dxψbasis, 1.0)
    @inbounds for (i, dimi) in enumerate(dimoff)
        for (j, dimj) in enumerate(dimoff)
            if i==j
            dxψbasis[:,:,i] .*= grad_xk_basis(R.f, X, 1, dimj, dimj, R.f.idx)
            else
            dxψbasis[:,:,i] .*= ψbasis[:,:,j]
            end
        end
    end

    # This integrand computes ψ1 ⊗ … ⊗ ψi′ ⊗ … ⊗ ψk′(t) × c g′(ψ1 ⊗ … ⊗ ψk′(t)×c)
    function integrand!(v::Vector{Float64}, t::Float64)
    dxkψ .= repeated_grad_xk_basis(R.f,  t*xlast)
    @avx @. cacheg = (dxkψ * ψoff) *ˡ coeff
    grad_x!(cacheg, R.g, cacheg)

        # @inbounds for i=1:Nx-1
        #     vi = view(v, (i-1)*Ne+1:i*Ne)
        #     dxψbasisi = view(dxψbasis,:,:,i)
        #     vi .= ((dxψbasisi .* dxkψ) * coeff) .* cacheg
        # end

        @inbounds for (i, dimi) in enumerate(dimoff)
            vi = view(v, (i-1)*Ne+1:i*Ne)
            dxψbasisi = view(dxψbasis,:,:,i)
            vi .= ((dxψbasisi .* dxkψ) * coeff) .* cacheg
        end
    end
    quadgk!(integrand!, cache, 0.0, 1.0)#; rtol = 1e-3)

    # Multiply integral by xlast (change of variable in the integration)
    @inbounds for (i, dimi) in enumerate(dimoff)
        @. cache[(i-1)*Ne+1:i*Ne] *= xlast
    end

    # Add ψ1 ⊗ … ⊗ ψi′ ⊗ … ⊗ ψk(0) × c
    @inbounds for (i, dimi) in enumerate(dimoff)
        dxψbasisi = view(dxψbasis,:,:,i)
        view(out,:,dimi) .+= (dxψbasisi .* ψk0)*coeff + view(cache,(i-1)*Ne+1:i*Ne)
    end

    # Fill the last column
    lastcol = view(out, :, Nx)
    evaluate!(lastcol, R.g, (repeated_grad_xk_basis(R.f,  xlast) .* ψoff)*coeff)
    return out
end

"""
$(TYPEDSIGNATURES)

Evaluates the gradient (with respect to the state) the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates ∂ₓ f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
"""
grad_x(R::IntegratedFunction, X) = grad_x!(zeros(size(X')), R, X)

"""
$(TYPEDSIGNATURES)

Evaluates in-place the gradient (with respect to the state) the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates ∂ₓ f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
Note that the gradient is only taken with respect to the active dimensions in the feature expansion.
"""
function reduced_grad_x!(out, R::IntegratedFunction, X)
    NxX, Ne = size(X)
    Nx = R.Nx

    dim = active_dim(R)
    dimoff = dim[dim .< Nx]

    @assert NxX == Nx "Wrong dimension of the sample"
    @assert size(out) == (Ne, length(dim)) "Dimensions of the output and the samples don't match"

    x0 = zeros(Ne)
    xlast = view(X,Nx,:)
    ψk0  = repeated_evaluate_basis(R.f, x0)
    ψoff = evaluate_offdiagbasis(R.f, X)
    dxkψ = zero(ψk0)

    coeff = R.f.coeff

    # Define active and off diagonal active dimension
    dim = active_dim(R)
    dimoff = dim[dim .< Nx]

    # Cache for the integration
    cache = zeros(length(dimoff)*Ne)
    cacheg = zeros(Ne)

    # Compute the basis for each component
    ψbasis = zeros(Ne, R.Nψ, length(dimoff))
    @inbounds for (i, dimi) in enumerate(dimoff)
        ψbasis_i = view(ψbasis, :, :, i)
        ψbasis_i .= evaluate_basis(R.f, X, [dimi], R.f.idx)
    end

    # Compute ψ1 ⊗ ψ2 ⊗ ψi′ ⊗ … ⊗ ψk-1
    dxψbasis = zero(ψbasis)
    fill!(dxψbasis, 1.0)
    @inbounds for (i, dimi) in enumerate(dimoff)
        for (j, dimj) in enumerate(dimoff)
            if i==j
            dxψbasis[:,:,i] .*= grad_xk_basis(R.f, X, 1, dimj, dimj, R.f.idx)
            else
            dxψbasis[:,:,i] .*= ψbasis[:,:,j]
            end
        end
    end

    # This integrand computes ψ1 ⊗ … ⊗ ψi′ ⊗ … ⊗ ψk′(t) × c g′(ψ1 ⊗ … ⊗ ψk′(t)×c)
    function integrand!(v::Vector{Float64}, t::Float64)
    dxkψ .= repeated_grad_xk_basis(R.f,  t*xlast)
    @avx @. cacheg = (dxkψ * ψoff) *ˡ coeff
    grad_x!(cacheg, R.g, cacheg)

        # @inbounds for i=1:Nx-1
        #     vi = view(v, (i-1)*Ne+1:i*Ne)
        #     dxψbasisi = view(dxψbasis,:,:,i)
        #     vi .= ((dxψbasisi .* dxkψ) * coeff) .* cacheg
        # end

        @inbounds for (i, dimi) in enumerate(dimoff)
            vi = view(v, (i-1)*Ne+1:i*Ne)
            dxψbasisi = view(dxψbasis,:,:,i)
            vi .= ((dxψbasisi .* dxkψ) * coeff) .* cacheg
        end
    end
    quadgk!(integrand!, cache, 0.0, 1.0)#; rtol = 1e-3)

    # Multiply integral by xlast (change of variable in the integration)
    @inbounds for (i, dimi) in enumerate(dimoff)
        @. cache[(i-1)*Ne+1:i*Ne] *= xlast
    end

    # Add ψ1 ⊗ … ⊗ ψi′ ⊗ … ⊗ ψk(0) × c
    @inbounds for (i, dimi) in enumerate(dimoff)
        dxψbasisi = view(dxψbasis,:,:,i)
        view(out,:,i) .+= (dxψbasisi .* ψk0)*coeff + view(cache,(i-1)*Ne+1:i*Ne)
    end

    if dim[end] == Nx
        # Fill the last column
        lastcol = view(out, :, length(dim))
        evaluate!(lastcol, R.g, (repeated_grad_xk_basis(R.f,  xlast) .* ψoff)*coeff)
    end
    return out
end

"""
$(TYPEDSIGNATURES)

Evaluates the gradient (with respect to the state) the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates ∂ₓ f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
Note that the gradient is only taken with respect to the active dimensions in the feature expansion.
"""
reduced_grad_x(R::IntegratedFunction, X) = reduced_grad_x!(zeros(size(X,2), length(active_dim(R))), R, X)

## Hessian of an IntegratedFunction

"""
$(TYPEDSIGNATURES)

Evaluates in-place the hessian (with respect to the state) the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates ∂²ₓ f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
"""
function hess_x!(out, R::IntegratedFunction, X)
    NxX, Ne = size(X)
    Nx = R.Nx
    Nψ = R.Nψ
    @assert NxX == Nx "Wrong dimension of the sample"
    @assert size(out) == (Ne, Nx, Nx) "Dimensions of the output and the samples don't match"

    x0 = zeros(Ne)
    xlast = view(X,Nx,:)
    ψk0  = repeated_evaluate_basis(R.f, x0)
    ψoff = evaluate_offdiagbasis(R.f, X)
    dxkψ = zero(ψk0)

    dgψ = zeros(Ne)

    coeff = R.f.coeff

    # Define active and off diagonal active dimension
    dim = active_dim(R)
    dimoff = dim[dim .< Nx]

    if Nx>1

        # Cache for the integration
        cache = zeros(length(dimoff)*Ne)
        cachedg  = zeros(Ne)
        cached2g = zeros(Ne)

        # Compute the basis for each component ψi(xi)
        ψbasis = zeros(Ne, R.Nψ, length(dimoff))
        @inbounds for (i, dimi) in enumerate(dimoff)
            ψbasis_i = view(ψbasis, :, :, i)
            ψbasis_i .= evaluate_basis(R.f, X, [dimi], R.f.idx)
        end

        dxψbasis = zero(ψbasis)
        d2xψbasis = zero(ψbasis)
        fill!(dxψbasis, 1.0)
        fill!(d2xψbasis, 1.0)

        # Compute ψ1 ⊗ ψ2 ⊗ ψi′⊗ … ⊗ ψk-1 && ψ1 ⊗ ψ2 ⊗ ψi″⊗ … ⊗ ψk-1
        @inbounds for (i, dimi) in enumerate(dimoff)
            for (j, dimj) in enumerate(dimoff)
                if i==j
                dxψbasis[:,:,i] .*= grad_xk_basis(R.f, X, 1, dimj, dimj, R.f.idx)
                d2xψbasis[:,:,i] .*= grad_xk_basis(R.f, X, 2, dimj, dimj, R.f.idx)

                else
                ψbasisj = view(ψbasis, :, :, j)
                @avx @. dxψbasis[:,:,i] *= ψbasisj#[:,:,j]
                @avx @. d2xψbasis[:,:,i] .*= ψbasisj#[:,:,j]
                end
            end
        end
    end

    ##########################################################################################
    # Compute ∂^2_ij R(x1:k) = ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψj′ ⊗ … ⊗ ψk(0) c + ∫_0^x_k [ψ1 ⊗ … ⊗ ψi′ ⊗ … ⊗ ψj′ ⊗ … ⊗ ψk′(t)c
    # g′(ψ1 ⊗ … ⊗ ψk′(t)c) +  ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψk′(t) c ⨂ ψ1 ⊗ … ⊗ ψj′⊗ … ⊗ ψk′(t) c
    # g″(ψ1 ⊗ … ⊗ ψk′(t)c) dt for i,j ∈[1,k-1] and i<j

    if Nx>2
        cacheij = zeros(ceil(Int64, (length(dimoff))*(length(dimoff)-1)*Ne/2))
        dxijψbasis = zeros(Ne, Nψ, ceil(Int64, (length(dimoff))*(length(dimoff)-1)/2))
        fill!(dxijψbasis, 1.0)
        # Compute ψ1 ⊗ ψ2 ⊗ ψi′⊗ … ⊗ ψj′⊗ … ⊗ ψk-1
        count = 0
        @inbounds for i=1:length(dimoff)
            for j=i+1:length(dimoff)
                count += 1
                dxijψbasis[:,:,count] .*= grad_xk_basis(R.f, X, 1, [dimoff[i]; dimoff[j]], [dimoff[i]; dimoff[j]], R.f.idx)
                for k=1:length(dimoff)
                    if dimoff[k] != dimoff[i] && dimoff[k] != dimoff[j]
                        ψbasisk = view(ψbasis, :, :, k)
                        @avx @. dxijψbasis[:,:,count] .*= ψbasisk
                    end
                end
            end
        end

        function integrandij!(v::Vector{Float64}, t::Float64)
        dxkψ .= repeated_grad_xk_basis(R.f,  t*xlast)
        @avx @. cachedg = (dxkψ * ψoff) *ˡ coeff
        hess_x!(cached2g, R.g, cachedg)
        grad_x!(cachedg, R.g, cachedg)
            count = 0
            @inbounds for i=1:length(dimoff)
                for j=i+1:length(dimoff)
                    count +=1
                    vij= view(v, (count-1)*Ne+1:count*Ne)
                    dxijψbasis_count = view(dxijψbasis,:,:,count)
                    dxψbasisi = view(dxψbasis,:,:,i)
                    dxψbasisj = view(dxψbasis,:,:,j)
                    vij .= ((dxijψbasis_count .* dxkψ) * coeff) .* cachedg + (((dxψbasisi .* dxkψ) * coeff) .* ((dxψbasisj .* dxkψ) * coeff)) .* cached2g
                end
            end
        end

        quadgk!(integrandij!, cacheij, 0.0, 1.0)#; rtol = 1e-3 )

        # Multiply integral by xlast (change of variable in the integration)
        @inbounds for i=1:length(dimoff)
            @. cacheij[(i-1)*Ne+1:i*Ne] *= xlast
        end
        count = 0
        @inbounds for i=1:length(dimoff)
            for j=i+1:length(dimoff)
                count += 1
                colij = view(out, :, dimoff[i], dimoff[j])
                colji = view(out, :, dimoff[j], dimoff[i])
                cacheij_count = view(cacheij, (count-1)*Ne+1:count*Ne)
                dxijψbasis_count = view(dxijψbasis,:,:,count)
                @avx @. colij = (dxijψbasis_count * ψk0) *ˡ coeff
                colij .+= cacheij_count
                colji .= colij
            end
        end
    end
    ##########################################################################################
    # Compute ∂^2_ii R(x1:k) = ψ1 ⊗ … ⊗ ψi″⊗ … ⊗ ψk(0) c + ∫_0^x_k [ψ1 ⊗ … ⊗ ψi″ ⊗ … ⊗ ψk′(t)c
    # g′(ψ1 ⊗ … ⊗ ψk′(t)c) +  ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψk′(t) c ⨂ ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψk′(t) c
    # g″(ψ1 ⊗ … ⊗ ψk′(t)c) dt for i ∈[1,k-1]
    if Nx>1
        # Compute integral term
        cache = zeros(length(dimoff)*Ne)

        function integrandii!(v::Vector{Float64}, t::Float64)
        dxkψ .= repeated_grad_xk_basis(R.f,  t*xlast)
        @avx @. cachedg = (dxkψ * ψoff) *ˡ coeff
        hess_x!(cached2g, R.g, cachedg)
        grad_x!(cachedg, R.g, cachedg)

            @inbounds for i=1:length(dimoff)
                vi = view(v, (i-1)*Ne+1:i*Ne)
                dxψbasisi = view(dxψbasis,:,:,i)
                d2xψbasisi = view(d2xψbasis,:,:,i)
                # vi .= ((d2xψbasisi .* dxkψ) * coeff) .* cachedg + (((dxψbasisi .* dxkψ) * coeff) .^2) .* cached2g
                @avx @. vi = ((d2xψbasisi * dxkψ) *ˡ coeff) * cachedg + (((dxψbasisi * dxkψ) *ˡ coeff) ^2) * cached2g
            end
        end

        quadgk!(integrandii!, cache, 0.0, 1.0)#; rtol = 1e-3)

        # Multiply integral by xlast (change of variable in the integration)
        @inbounds for i=1:length(dimoff)
            @. cache[(i-1)*Ne+1:i*Ne] *= xlast
        end

        @inbounds for i=1:length(dimoff)
            colii = view(out, :, dimoff[i], dimoff[i])
            cachei = view(cache, (i-1)*Ne+1:i*Ne)
            d2xψbasisi = view(d2xψbasis,:,:,i)
            @avx @. colii = (d2xψbasisi * ψk0) *ˡ coeff
            colii .+= cachei
        end
    end
    #############################################################################
    # Compute ∂^2_ik R(x1:k) = ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψk′(xk) c g′(ψ1 ⊗ … ⊗ ψk′(xk) c)
    if Nx>1
        dxkψk = repeated_grad_xk_basis(R.f,  xlast)
        grad_x!(dgψ, R.g, (ψoff .* dxkψk)* coeff)

        @inbounds for i=1:length(dimoff)
            colik = view(out,:,dimoff[i],Nx)
            colki = view(out,:,Nx,dimoff[i])
            dxψbasisi = view(dxψbasis,:,:,i)
            @avx @. colik = ((dxψbasisi * dxkψk) *ˡ coeff) * dgψ
            colki .= colik
        end
    end
    #############################################################################
    # Compute ∂^2_k R^k(x1:k) = ψ1 ⊗ … ⊗ ψk″(xk) c g′(ψ1 ⊗ … ⊗ ψk′(xk) c)

    d2xkψk = grad_xk_basis(R.f, X, 2, Nx, Nx, R.f.idx)
    colkk = view(out, :, Nx, Nx)
    if Nx>1
        @avx @. colkk = ((ψoff * d2xkψk) *ˡ coeff) * dgψ
    else
        dxkψk = repeated_grad_xk_basis(R.f,  xlast)
        grad_x!(dgψ, R.g, (dxkψk)* coeff)
        @avx @. colkk = (d2xkψk *ˡ coeff) * dgψ
    end
    return out
end

"""
$(TYPEDSIGNATURES)

Evaluates the hessian (with respect to the state) the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates ∂²ₓ f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
"""
hess_x(R::IntegratedFunction, X) = hess_x!(zeros(size(X,2), size(X,1), size(X,1)), R, X)

# This version outputs a result of size(Ne, active_dim(R), active_dim(R))
"""
$(TYPEDSIGNATURES)

Evaluates in-place the hessian (with respect to the state) the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates ∂²ₓ f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
Note that the hessian is only taken with respect to the active dimensions in the feature expansion.
"""
function reduced_hess_x!(out, R::IntegratedFunction, X)
    NxX, Ne = size(X)
    Nx = R.Nx
    Nψ = R.Nψ

    dim = active_dim(R)
    dimoff = dim[dim .< Nx]

    @assert NxX == Nx "Wrong dimension of the sample"
    @assert size(out) == (Ne, length(dim), length(dim)) "Dimensions of the output and the samples don't match"

    x0 = zeros(Ne)
    xlast = view(X,Nx,:)
    ψk0  = repeated_evaluate_basis(R.f, x0)
    ψoff = evaluate_offdiagbasis(R.f, X)
    dxkψ = zero(ψk0)

    dgψ = zeros(Ne)

    coeff = R.f.coeff

    # Define active and off diagonal active dimension

    if Nx>1

        # Cache for the integration
        cache = zeros(length(dimoff)*Ne)
        cachedg  = zeros(Ne)
        cached2g = zeros(Ne)

        # Compute the basis for each component ψi(xi)
        ψbasis = zeros(Ne, R.Nψ, length(dimoff))
        @inbounds for (i, dimi) in enumerate(dimoff)
            ψbasis_i = view(ψbasis, :, :, i)
            ψbasis_i .= evaluate_basis(R.f, X, [dimi], R.f.idx)
        end

        dxψbasis = zero(ψbasis)
        d2xψbasis = zero(ψbasis)
        fill!(dxψbasis, 1.0)
        fill!(d2xψbasis, 1.0)

        # Compute ψ1 ⊗ ψ2 ⊗ ψi′⊗ … ⊗ ψk-1 && ψ1 ⊗ ψ2 ⊗ ψi″⊗ … ⊗ ψk-1
        @inbounds for (i, dimi) in enumerate(dimoff)
            for (j, dimj) in enumerate(dimoff)
                if i==j
                dxψbasis[:,:,i] .*= grad_xk_basis(R.f, X, 1, dimj, dimj, R.f.idx)
                d2xψbasis[:,:,i] .*= grad_xk_basis(R.f, X, 2, dimj, dimj, R.f.idx)

                else
                ψbasisj = view(ψbasis, :, :, j)
                @avx @. dxψbasis[:,:,i] *= ψbasisj#[:,:,j]
                @avx @. d2xψbasis[:,:,i] .*= ψbasisj#[:,:,j]
                end
            end
        end
    end

    ##########################################################################################
    # Compute ∂^2_ij R(x1:k) = ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψj′ ⊗ … ⊗ ψk(0) c + ∫_0^x_k [ψ1 ⊗ … ⊗ ψi′ ⊗ … ⊗ ψj′ ⊗ … ⊗ ψk′(t)c
    # g′(ψ1 ⊗ … ⊗ ψk′(t)c) +  ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψk′(t) c ⨂ ψ1 ⊗ … ⊗ ψj′⊗ … ⊗ ψk′(t) c
    # g″(ψ1 ⊗ … ⊗ ψk′(t)c) dt for i,j ∈[1,k-1] and i<j

    if Nx>2
        cacheij = zeros(ceil(Int64, (length(dimoff))*(length(dimoff)-1)*Ne/2))
        dxijψbasis = zeros(Ne, Nψ, ceil(Int64, (length(dimoff))*(length(dimoff)-1)/2))
        fill!(dxijψbasis, 1.0)
        # Compute ψ1 ⊗ ψ2 ⊗ ψi′⊗ … ⊗ ψj′⊗ … ⊗ ψk-1
        count = 0
        @inbounds for i=1:length(dimoff)
            for j=i+1:length(dimoff)
                count += 1
                dxijψbasis[:,:,count] .*= grad_xk_basis(R.f, X, 1, [dimoff[i]; dimoff[j]], [dimoff[i]; dimoff[j]], R.f.idx)
                for k=1:length(dimoff)
                    if dimoff[k] != dimoff[i] && dimoff[k] != dimoff[j]
                        ψbasisk = view(ψbasis, :, :, k)
                        @avx @. dxijψbasis[:,:,count] .*= ψbasisk
                    end
                end
            end
        end

        function integrandij!(v::Vector{Float64}, t::Float64)
        dxkψ .= repeated_grad_xk_basis(R.f,  t*xlast)
        @avx @. cachedg = (dxkψ * ψoff) *ˡ coeff
        hess_x!(cached2g, R.g, cachedg)
        grad_x!(cachedg, R.g, cachedg)
            count = 0
            @inbounds for i=1:length(dimoff)
                for j=i+1:length(dimoff)
                    count +=1
                    vij= view(v, (count-1)*Ne+1:count*Ne)
                    dxijψbasis_count = view(dxijψbasis,:,:,count)
                    dxψbasisi = view(dxψbasis,:,:,i)
                    dxψbasisj = view(dxψbasis,:,:,j)
                    vij .= ((dxijψbasis_count .* dxkψ) * coeff) .* cachedg + (((dxψbasisi .* dxkψ) * coeff) .* ((dxψbasisj .* dxkψ) * coeff)) .* cached2g
                    # @avx @. vij = ((dxijψbasis_count * dxkψ) *ˡ coeff) * cachedg + (((dxψbasisi * dxkψ) *ˡ coeff) * ((dxψbasisj * dxkψ) *ˡ coeff)) * cached2g
                end
            end
        end

        quadgk!(integrandij!, cacheij, 0.0, 1.0)#; rtol = 1e-3 )

        # Multiply integral by xlast (change of variable in the integration)
        @inbounds for i=1:length(dimoff)
            @. cacheij[(i-1)*Ne+1:i*Ne] *= xlast
        end
        count = 0
        @inbounds for i=1:length(dimoff)
            for j=i+1:length(dimoff)
                count += 1
                # colij = view(out, :, dimoff[i], dimoff[j])
                # colji = view(out, :, dimoff[j], dimoff[i])
                colij = view(out, :, i, j)
                colji = view(out, :, j, i)
                cacheij_count = view(cacheij, (count-1)*Ne+1:count*Ne)
                dxijψbasis_count = view(dxijψbasis,:,:,count)
                @avx @. colij = (dxijψbasis_count * ψk0) *ˡ coeff
                colij .+= cacheij_count
                colji .= colij
            end
        end
    end
    ##########################################################################################
    # Compute ∂^2_ii R(x1:k) = ψ1 ⊗ … ⊗ ψi″⊗ … ⊗ ψk(0) c + ∫_0^x_k [ψ1 ⊗ … ⊗ ψi″ ⊗ … ⊗ ψk′(t)c
    # g′(ψ1 ⊗ … ⊗ ψk′(t)c) +  ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψk′(t) c ⨂ ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψk′(t) c
    # g″(ψ1 ⊗ … ⊗ ψk′(t)c) dt for i ∈[1,k-1]
    if Nx>1
        # Compute integral term
        cache = zeros(length(dimoff)*Ne)

        function integrandii!(v::Vector{Float64}, t::Float64)
        dxkψ .= repeated_grad_xk_basis(R.f,  t*xlast)
        @avx @. cachedg = (dxkψ * ψoff) *ˡ coeff
        hess_x!(cached2g, R.g, cachedg)
        grad_x!(cachedg, R.g, cachedg)

            @inbounds for i=1:length(dimoff)
                vi = view(v, (i-1)*Ne+1:i*Ne)
                dxψbasisi = view(dxψbasis,:,:,i)
                d2xψbasisi = view(d2xψbasis,:,:,i)
                # vi .= ((d2xψbasisi .* dxkψ) * coeff) .* cachedg + (((dxψbasisi .* dxkψ) * coeff) .^2) .* cached2g
                @avx @. vi = ((d2xψbasisi * dxkψ) *ˡ coeff) * cachedg + (((dxψbasisi * dxkψ) *ˡ coeff) ^2) * cached2g
            end
        end

        quadgk!(integrandii!, cache, 0.0, 1.0)#; rtol = 1e-3)

        # Multiply integral by xlast (change of variable in the integration)
        @inbounds for i=1:length(dimoff)
            @. cache[(i-1)*Ne+1:i*Ne] *= xlast
        end

        @inbounds for i=1:length(dimoff)
            # colii = view(out, :, dimoff[i], dimoff[i])
            colii = view(out, :, i, i)
            cachei = view(cache, (i-1)*Ne+1:i*Ne)
            d2xψbasisi = view(d2xψbasis,:,:,i)
            @avx @. colii = (d2xψbasisi * ψk0) *ˡ coeff
            colii .+= cachei
        end
    end

    if dim[end] == Nx
        #############################################################################
        # Compute ∂^2_ik R(x1:k) = ψ1 ⊗ … ⊗ ψi′⊗ … ⊗ ψk′(xk) c g′(ψ1 ⊗ … ⊗ ψk′(xk) c)
        if Nx>1
            dxkψk = repeated_grad_xk_basis(R.f,  xlast)
            grad_x!(dgψ, R.g, (ψoff .* dxkψk)* coeff)

            @inbounds for i=1:length(dimoff)
                # colik = view(out,:,dimoff[i],Nx)
                # colki = view(out,:,Nx,dimoff[i])
                colik = view(out,:,i,length(dim))
                colki = view(out,:,length(dim),i)
                dxψbasisi = view(dxψbasis,:,:,i)
                @avx @. colik = ((dxψbasisi * dxkψk) *ˡ coeff) * dgψ
                colki .= colik
            end
        end
        #############################################################################
        # Compute ∂^2_k R^k(x1:k) = ψ1 ⊗ … ⊗ ψk″(xk) c g′(ψ1 ⊗ … ⊗ ψk′(xk) c)

        d2xkψk = grad_xk_basis(R.f, X, 2, Nx, Nx, R.f.idx)
        colkk = view(out, :, length(dim), length(dim))
        if Nx>1
            @avx @. colkk = ((ψoff * d2xkψk) *ˡ coeff) * dgψ
        else
            dxkψk = repeated_grad_xk_basis(R.f,  xlast)
            grad_x!(dgψ, R.g, (dxkψk)* coeff)
            @avx @. colkk = (d2xkψk *ˡ coeff) * dgψ
        end
    end
    return out
end

"""
$(TYPEDSIGNATURES)

Evaluates the hessian (with respect to the state) the `IntegratedFunction` `R` for the ensemble matrix `X`, i.e.
evaluates ∂²ₓ f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
Note that the hessian is only taken with respect to the active dimensions in the feature expansion.
"""
reduced_hess_x(R::IntegratedFunction, X) = reduced_hess_x!(zeros(size(X,2), length(active_dim(R.f)), length(active_dim(R.f))), R, X)

## Compute ∂_c int_0^{x_k} g(∂ₖf(x_{1:k-1}, t))dt
"""
$(TYPEDSIGNATURES)

Evaluates the following gradient with respect to the coefficients of `R`:  ∂_c ∫_0^{x_k} g(∂ₖf(x_{1:k-1}, t))dt for the different columns of the ensemble matrix `X`.
"""
function grad_coeff_integrate_xd(R::IntegratedFunction, X::Array{Float64,2})
    ψoff = evaluate_offdiagbasis(R.f, X)
    Nx = R.Nx
    Nψ = R.Nψ
    NxX, Ne = size(X)
    @assert NxX == Nx "Wrong dimension of the sample X"
    xk = deepcopy(X[Nx, :])
    dcdψ = zeros(Ne, Nψ)

    cache = zeros(Ne, Nψ)

    function integrand!(v::Matrix{Float64}, t::Float64)
        dcdψ .= repeated_grad_xk_basis(R.f,  t*xk) .* ψoff

        v .= vgrad_x(R.g, dcdψ*R.f.coeff) .* dcdψ
    end

    return xk .* quadgk!(integrand!, cache, 0.0, 1.0)[1]
end

# Compute ∂²_c int_0^{x_k} g(c^T F(t))dt = int_0^{x_k} ∂_c F(t) g′(c^T F(t)) + F(t)F(t)*g″(c^T F(t)) dt
"""
$(TYPEDSIGNATURES)

Evaluates the following hessian with respect to the coefficients of `R`:  ∂^2c ∫_0^{x_k} g(∂ₖf(x_{1:k-1}, t))dt for the different columns of the ensemble matrix `X`.
"""
function hess_coeff_integrate_xd(R::IntegratedFunction, X::Array{Float64,2})
     # We drop the first term for speed improvement since it is always equal to 0
    Nx = R.Nx
    Nψ = R.Nψ
    NxX, Ne = size(X)
    @assert NxX == Nx "Wrong dimension of the sample X"
    ψoff = evaluate_offdiagbasis(R.f, X)
    xk = deepcopy(X[Nx, :])
    dcdψ = zeros(Ne, Nψ)

    dcdψouter = zeros(Ne, Nψ, Nψ)

    cache = zeros(Ne*Nψ*Nψ)

    function integrand!(v::Vector{Float64}, t::Float64)
        dcdψ .= repeated_grad_xk_basis(R.f,  t*xk) .* ψoff
        @inbounds for i=1:Nψ
            for j=1:Nψ
                dcdψouter[:,i,j] = dcdψ[:,i] .* dcdψ[:, j]
            end
        end
        v .= reshape(vhess_x(R.g, (dcdψ ) * R.f.coeff) .* dcdψouter, (Ne*Nψ*Nψ))
    end

    return xk .* reshape(quadgk!(integrand!, cache, 0.0, 1.0)[1], (Ne, Nψ, Nψ))
end

"""
$(TYPEDSIGNATURES)

Evaluates the gradient (with respect to the coefficients) the `IntegratedFunction` `R` for the different columns of the ensemble matrix `X`, i.e.
evaluates ∂_c f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
"""
function grad_coeff(R::IntegratedFunction, X::Array{Float64,2})
    Nx = R.Nx
    Nψ = R.Nψ
    NxX, Ne = size(X)
    @assert NxX == Nx "Wrong dimension of the sample X"

    ψoff = evaluate_offdiagbasis(R.f, X)
    ψdiag = repeated_evaluate_basis(R.f, zeros(Ne))

    xk = deepcopy(X[Nx, :])
    dcdψ = zeros(Ne, Nψ)

    cache = zeros(Ne, Nψ)

    function integrand!(v::Matrix{Float64}, t::Float64)
        dcdψ .= repeated_grad_xk_basis(R.f,  t*xk) .* ψoff
        v .= vgrad_x(R.g, dcdψ*R.f.coeff) .* dcdψ
    end

    return ψoff .* ψdiag + xk .* quadgk!(integrand!, cache, 0.0, 1.0)[1]
end

"""
$(TYPEDSIGNATURES)

Evaluates the hessian (with respect to the coefficients) the `IntegratedFunction` `R` for the different columns of the ensemble matrix `X`, i.e.
evaluates ∂²_c f(x_{1:k-1}, 0) + ∫₀ xₖ g(∂ₖ f(x_{1:k-1}, t)) dt for the different columns of `X`.
"""
hess_coeff(R::IntegratedFunction, X::Array{Float64,2}) = hess_coeff_integrate_xd(R, X)


text:
using PSO
using Base.Test
using FactCheck

facts("Example minimization of 4th-order banana function") do
    function myfunc(x)
        x1, x2 = x[1], x[2]
        return x1^4 - 2*x2*x1^2 + x2^2 + x1^2 - 2*x1 + 5
    end

    lb = [-3.0, -1.0]
    ub = [2.0, 6.0]

    context("no constraints") do
        xopt1, fopt1 = pso(myfunc, lb, ub)
        x, y = xopt1
        @fact x --> roughly(1.0; atol=1e-3)
        @fact y --> roughly(1.0; atol=1e-3)
        @fact fopt1 --> roughly(4.0, atol=1e-8)

        println("The optimum is at:")
        println("    $(xopt1)")
        println("Optimal function value:")
        println("    myfunc: $(fopt1)")
    end

    function mycon(x)
        x1, x2 = x[1], x[2]
        return [-(x1 + 0.25)^2 + 0.75*x2]
    end

    context("with constraint") do
        # xopt2, fopt2 = pso(myfunc, initial, bounds, num_particles, maxiter)
        xopt2, fopt2 = pso(myfunc, lb, ub, constraints=mycon, swarmsize=1000)
        x, y = xopt2
        @fact x --> roughly(0.5; atol=1e-2)
        @fact y --> roughly(0.75; atol=1e-1)
        @fact fopt2 --> roughly(4.5, atol=1e-3)

        println("The optimum is at:")
        println("    $(xopt2)")
        println("Optimal function value:")
        println("    myfunc: $(fopt2)")
        println("    mycon : $(mycon(xopt2))")
    end
end

facts("Engineering example: minimization of twobar truss weight, subject to") do
    println("  Yield Stress <= 100 kpsi")
    println("  Yield Stress <= Buckling Stress")
    println("  Deflection   <= 0.25 inches")

    function weight(x, B, ρ, E, P)
        H, d, t = x  # all in inches
        return ρ*2*π*d*t*√((B/2)^2 + H^2)
    end

    function stress(x, B, ρ, E, P)
        H, d, t = x  # all in inches
        return (P * √((B/2)^2 + H^2)) / (2*t*π*d*H)
    end

    function buckling_stress(x, B, ρ, E, P)
        H, d, t = x  # all in inches
        return (π^2 * E * (d^2 + t^2)) / (8 * ((B/2)^2 + H^2))
    end

    function deflection(x, B, ρ, E, P)
        H, d, t = x  # all in inches
        return (P * (√((B/2)^2 + H^2))^3) / (2*t*π*d*H^2 * E)
    end

    function mycons(x, B, ρ, E, P)
        strs = stress(x, B, ρ, E, P)
        buck = buckling_stress(x, B, ρ, E, P)
        defl = deflection(x, B, ρ, E, P)
        return [100 - strs, buck - strs, 0.25 - defl]
    end

    B = 60  # inches
    rho = 0.3  # lb/in^3
    E = 30000  # kpsi
    P = 66  # lb (force)
    args = (B, rho, E, P)
    lb = [10.0, 1.0, 0.01]
    ub = [30.0, 3.0, 0.25]
    context("") do
        xopt4, fopt4 = pso(weight, lb, ub, constraints=mycons, args=args, swarmsize=1000)
        x, y, z = xopt4
        @fact x --> roughly(29.0; atol=1.0)
        @fact y --> roughly(2.5; atol=1.0)
        @fact fopt4 --> roughly(11.88, atol=1e-2)

        println("The optimum is at:")
        println("    $(xopt4)")
        println("Optimal function values:")
        println("    weight         : $(fopt4)")
        println("Constraint functions:")
        println("    stress         : $(stress(xopt4, args...))")
        println("    buckling stress: $(buckling_stress(xopt4, args...))")
        println("    deflection     : $(deflection(xopt4, args...))")
    end
end

ext: jl


text:
## ---------------------------------------------------------------
function plot_dead_explosion_map(simid)

    params = lglob(Dyn, simid, :params, :infinite_cg)
    @show params
    @extract params: Ds ϵs cg

    Ds_idxs = Dict(D => i for (i, D) in enumerate(Ds))
    ϵs_idxs = Dict(ϵ => i for (i, ϵ) in enumerate(ϵs))

    colGRAD = cgrad([colorant"black", colorant"blue", colorant"yellow"]);
    mat = fill(3, length(Ds_idxs), length(ϵs_idxs))

    for (D, ϵ) in Iterators.product(Ds, ϵs)
        simparams = (;D, ϵ, cg)

        status = get_status(simid, simparams)

        if status == DEAD_SIM_STATUS
            mat[Ds_idxs[D], ϵs_idxs[ϵ]] = 1
        elseif status == STST_SIM_STATUS
            mat[Ds_idxs[D], ϵs_idxs[ϵ]] = 2
        elseif status == EXPLODED_SIM_STATUS
            mat[Ds_idxs[D], ϵs_idxs[ϵ]] = 3
        else
            @warn("Unexpected status: ", simparams, status)
        end
    end

    @info("At", cg)
    # title = _textbf("c_{g}=Inf"),
    p = plot(;
        title = _textbf(""),
        xlabel = _textbf("D"),
        ylabel = _textbf("\\epsilon")
    )
    # labels
    for (label, color, islast) = [
            (_textbf("dead"), :black, false),
            (_textbf("steady state"), :blue, false),
            (_textbf("explosion"), :yellow, true),
        ]
        scatter!(p, [Ds[end ÷ 2]], [ϵs[end ÷ 2]];
            markerstrokewidth=0.3,
            label, color, m = :square
        )
        # fake label for spacing
        # !islast && scatter!(p, [Ds[end ÷ 2]], [ϵs[end ÷ 2]]; label = _textbf(" "), alpha = 0.0)
    end

    # map
    heatmap!(p, Ds, ϵs, mat';
        color = colGRAD,
        colorbar = false,
        label = "",
    )

    # sfig(ChE, p,
    #     simid, "D_ϵ_heatmap", (;cg), ".png"
    # )

    return p
end


text:
@everywhere using DataStructures
using Base.Test
using ArrayViews

type InverseHessian
  syrhos::Deque{(Vector{Float64}, Vector{Float64}, Float64)}
  tau::Int
end

function InverseHessian(tau::Int)
  return InverseHessian(Deque{(Vector{Float64}, Vector{Float64}, Float64)}(), tau)
end

function add!(self::InverseHessian, s::Vector{Float64}, y::Vector{Float64})
  rho = 1./dot(s, y)
  push!(self.syrhos, (s, y, rho))
  if length(self.syrhos) > self.tau
    shift!(self.syrhos)
  end
end

function test_mvp(self::InverseHessian, g::Vector{Float64}, delta::Float64)
  syrhos = collect(self.syrhos)
  (s, y, rho) = syrhos[end]
  p = length(s)
  gamma = dot(s, y)/dot(y, y)
  H = gamma / (1 + delta*gamma) * eye(p)
  for (s, y, rho) in syrhos
    V = eye(p) - rho * y * s'
    H = V' * H * V + rho * s * s'
  end
  return H*g
end

function mvp(self::InverseHessian, g::Vector{Float64}, delta::Float64)
  @assert length(self.syrhos) > 0
  q = copy(g)
  alphas = zeros(length(self.syrhos))
  for (i, (s, y, rho)) in enumerate(reverse(collect(self.syrhos))) # TODO: Implement reverse iterator for Deques
    alpha = rho * dot(s, q)
     alphas[length(self.syrhos) - i + 1] = alpha # TODO: implement reverse(enumerate)
     q -= alpha * y
  end
  (s, y, rho) = back(self.syrhos) # current iterate
  gamma = dot(s, y)/dot(y, y)
  r = gamma/(1 + delta*gamma)*q
  for (i, (s, y, rho)) in enumerate(collect(self.syrhos)) # TODO: Remove collect overhead
    beta = rho * dot(y, r)
    r += s * (alphas[i] - beta)
  end
  return r
end

function test_inverse_hessian()
  self = InverseHessian(3)
  for i = 1:4
    s = rand(2)
    y = rand(2)
    add!(self, s, y)
  end
  g = rand(2)
  first = mvp(self, g, 0.0)
  second = test_mvp(self, g, 0.0)
  @test norm(first - second) <= 1e-10
end

test_inverse_hessian()

# S is the number of epochs
function slbfgs(loss::Function, grad!::Function, hvp!::Function, data::Matrix, theta::Vector, config::Dict, diagnostics::HDF5Group)
  tau = config["memory_size"]
  delta = config["delta"]
  eta = config["stepsize"]
  S = config["epochs"]
  verbose = true
  println("===== running stochastic lbfgs =====\n")
  println("iter | lbfgs error")
  n = size(data, 1)
  m = 2*n
  p = length(theta)
  H = InverseHessian(tau)
  thetaold = copy(theta)
  dtheta = zeros(p)
  mu = zeros(p)
  z = zeros(p)
  zold = zeros(p)
  P = zeros(p)

  t = -1 # number of correction pairs currently computed
  stheta = zeros(p)
  sthetaold = zeros(p)

  hdf_loss = g_create(diagnostics, "loss")
  hdf_theta = g_create(diagnostics, "theta")

  if config["use_gradient_step"]
    println("warning: using gradient steps")
  end

  for s in 1:S
    if verbose
      l = loss(data, theta, config["lambda"])
      @printf "%4d | %9.4E \n" s l
      hdf_loss[@sprintf "%04d" s] = l
      hdf_theta[@sprintf "%04d" s] = theta
    end

    ctheta = copy(theta) # theta checkpoint
    # \tilde \mu = 1/n \sum \nable f_i(\tilde theta)
    fill!(mu, 0.0)
    for i = 1:n
      grad!(vec(data[:,i]), ctheta, config["lambda"], mu)
    end
    mu /= n
    if s == 1
      # Do a gradient step to get LBFGS started
      copy!(P, -mu/norm(mu))
    end
    copy!(theta, ctheta)
    for k in 1:m
      index = rand(1:n)
      copy!(zold, z)
      fill!(z, 0.0)
      # z_k = \nabla f_{i_k}(\theta) - \nabla f_{i_k}(\tilde\theta) + \tilde \mu
      fill!(dtheta, 0.0)
      grad!(vec(data[:,index]), theta, config["lambda"], dtheta)
      z += dtheta
      fill!(dtheta, 0.0)
      grad!(vec(data[:,index]), ctheta, config["lambda"], dtheta)
      if config["reduce_grad_var"]
        z -= dtheta
        z += mu
      end

      stheta += theta

      if k > 1 # neccessary?
        if mod(k, config["hess_period"]) == 0
          stheta /= config["hess_period"]
          t += 1
        end
        if mod(k, config["hess_period"]) == 0 && t >= 1
          svec = stheta - sthetaold
          r = zeros(p)
          if config["use_hvp"]
            for i = 1:config["hess_samples"]
              index = rand(1:n)
              hvp!(vec(data[:,index]), stheta, svec, config["lambda"], r)
            end
            r /= config["hess_samples"]
          else
                # delta a regularization parameter
                r = z - zold + delta*svec
            end
            add!(H, svec, r)
            copy!(sthetaold, stheta)
        end
        if t > 1
          P = -eta * mvp(H, z, delta)
        else
          P = -eta*z
        end
      end
      # take step
      copy!(thetaold, theta)
      if config["use_gradient_step"]
        theta -= eta*z
      else
        theta += P
      end
    end
  end
  return theta
end

max_stars_repo_name: pcmoritz/slbfgs

max_issues_repo_name: pcmoritz/slbfgs

max_forks_repo_name: pcmoritz/slbfgs


text:
using LaTeXStrings
using Juqbox
using Plots

# Function to loop over ϵ values and create a plot of the objective function for pertrubed Hamiltonians
function ep_plot(pcof::Vector{Float64}, params:: Juqbox.objparams, wa::Working_Arrays, ep_vals::AbstractArray)
    results = zeros(length(ep_vals),4)

    for i =1:length(ep_vals)
        epsh = ep_vals[i]

        # Additive noise
        for j = 2:size(params.Hconst,2)
            params.Hconst[j,j] += 0.01*epsh*(10.0^(j-2))
        end

        obj, _, _, secondaryobjf, traceinfid = Juqbox.traceobjgrad(pcof,params,wa,false, true)
        results[i,1] = obj
        results[i,2] = secondaryobjf
        results[i,3] = traceinfid
        results[i,4] = traceinfid+secondaryobjf

        # Reset
        for j = 2:size(params.Hconst,2)
            params.Hconst[j,j] -= 0.01*epsh*(10.0^(j-2))
        end

    end

    pl1 = Plots.plot(ep_vals,results[:,1],yaxis=:log,xlabel = L"\epsilon",ylabel = "Objective Function")
    return results,pl1
end

function evalctrl_no_carrier(params::objparams, pcof0:: Array{Float64, 1}, jFunc:: Int64, D1::Int64)

    # Evaluate the ctrl functions on this grid in time
    nplot = round(Int64, params.T*32)
    # is this resolution sufficient for the lab frame ctrl functions so we can get meaningful FFTs?
    td = collect(range(0, stop = params.T, length = nplot+1))

    nfreq = length(params.Cfreq)
    offset = (jFunc-1)*2*D1
    nCoeff = 2*D1
    pcof = copy(pcof0[offset+1:offset+nCoeff])

    if (params.use_bcarrier)
        # B-splines with carrier waves (zero out the frequency to plot just the splines)
        splinepar = Juqbox.bcparams(params.T, D1, params.Ncoupled, params.Nunc, zeros(1,1), pcof)
    else
        # regular B-splines
        splinepar = Juqbox.splineparams(params.T, D1, 2*(params.Ncoupled + params.Nunc), pcof)
    end

    # define inline function to enable vectorization over t
    controlplot(t, splinefunc) = Juqbox.controlfunc(t, splinepar, splinefunc)

    fact = 1.0/(2*pi) # conversion factor to GHz
    fact = 1.0 # conversion factor to rad/ns
    pj = fact.*controlplot.(td, 0)
    qj = fact.*controlplot.(td, 1)
    return pj, qj

end

# Maximum shift in Hamiltonian (in rad*GHz)
ep_max = 2*pi*2e-2

# For plotting purposes
max_ep_sweep = 2*pi*3e-2
len = 1001
ep_vals = range(-max_ep_sweep,stop=max_ep_sweep,length=len)
freshOptim = true

# Usual optimization
nquad = 1
include("swap-02-risk-neutral.jl")
if(freshOptim)
    pcof = Juqbox.run_optimizer(prob, pcof0)
    results,pl1 = ep_plot(pcof, params, wa, ep_vals)
    data = zeros(size(results,1),4)
    data[:,1] = ep_vals
    for j = 1:3
        data[:,j+1] = results[:,j]
    end
    writedlm("usual_control_GLQ2.dat", pcof)
    writedlm("usual_optim_OF_sweep_GLQ2.dat", data)
else
    pcof = vec(readdlm("usual_control_GLQ2.dat"))
    results = readdlm("usual_optim_OF_sweep_GLQ2.dat")
end

pcof_old = copy(pcof)

# Plot control functions
scalefactor = 1000/(2*pi)
unitStr = "MHz"

nplot = round(Int64, params.T*32)
td = collect(range(0, stop = params.T, length = nplot+1))
p_NF_1,q_NF_1 = evalctrl_no_carrier(params, pcof_old, 1, D1)
p_NF_2,q_NF_2 = evalctrl_no_carrier(params, pcof_old, 2, D1)

pfunc1 = scalefactor .* p_NF_1
qfunc1 = scalefactor .* q_NF_1
pmax1 = maximum(abs.(pfunc1))
qmax1 = maximum(abs.(qfunc1))
pfunc2 = scalefactor .* p_NF_2
qfunc2 = scalefactor .* q_NF_2
pmax2 = maximum(abs.(pfunc2))
qmax2 = maximum(abs.(qfunc2))
pmax = maximum([pmax1,pmax2])
qmax = maximum([qmax1,qmax2])

fnt = Plots.font("Helvetica", 12)
lfnt = Plots.font("Helvetica", 12)
Plots.default(titlefont=fnt, guidefont=fnt, tickfont=fnt, legendfont=lfnt, linewidth=1.5, size=(650, 350))

titlestr = "Rotating frame NF ctrl " * " Max-p=" *@sprintf("%.3e", pmax) * " Max-q=" *@sprintf("%.3e", qmax) * " " * unitStr

pl_ctrl_NF = Plots.plot(td, pfunc1, lab=L"p_{1,1}(t)", title = titlestr, xlabel="Time [ns]",
                                  ylabel=unitStr, legend= :outerright, linewidth=1.5, legendfontsize=14)
# add in the control function for the anti-symmetric Hamiltonian
Plots.plot!(pl_ctrl_NF,td, qfunc1, lab=L"q_{1,1}(t)", linewidth=1.5, legendfontsize=14)
Plots.plot!(pl_ctrl_NF,td, pfunc2, lab=L"p_{1,2}(t)", linewidth=1.5, legendfontsize=14)
Plots.plot!(pl_ctrl_NF,td, qfunc2, lab=L"q_{1,2}(t)", linewidth=1.5, legendfontsize=14)

# save plots of control functions in rotating frame without carrier waves
Plots.savefig(pl_ctrl_NF,  "robust_comparison_"* @sprintf("%1.1e", ep_max/(2*pi)) * "_"*@sprintf("%3.1f",T)*"_T_"*@sprintf("%3d",maxIter)*"_iters_"*@sprintf("%d",nquad)*"_N_"*@sprintf("%d",D1)*"_D1_NF_ctrl.png")

# Risk-neutral optimization
nquad = 9
include("swap-02-risk-neutral.jl")
if(freshOptim)
    pcof = Juqbox.run_optimizer(prob, pcof0)
    results2,pl2 = ep_plot(pcof, params, wa, ep_vals)
    data2 = zeros(size(results2,1),4)
    data2[:,1] = ep_vals
    for j = 1:3
        data2[:,j+1] = results2[:,j]
    end
    writedlm("robust_control_GLQ2.dat", pcof)
    writedlm("robust_optim_OF_sweep_GLQ2.dat", data2)
else
    pcof = vec(readdlm("robust_control_GLQ2.dat"))
    results2 = readdlm("robust_optim_OF_sweep_GLQ2.dat")
end

# Plot objective functions for various Hamiltonian perturbations
plc = Plots.plot(scalefactor.*ep_vals,results[:,4],yaxis=:log,xlabel = "Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="NF Infidelity")
Plots.plot!(plc,scalefactor.*ep_vals,results2[:,4],yaxis=:log,xlabel = "Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="RN Infidelity")
Plots.plot!(plc,scalefactor.*ep_vals,results[:,3],yaxis=:log,xlabel = "Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="NF Guard Level Pop.",linestyle=:dash)
Plots.plot!(plc,scalefactor.*ep_vals,results2[:,3],yaxis=:log,xlabel = "Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="RN Guard Level Pop.",linestyle=:dash,legend= :outerright)

plc_short = Plots.plot(scalefactor.*ep_vals,results[:,4],yaxis=:log,xlabel = "Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="NF Infidelity")
Plots.plot!(plc_short,scalefactor.*ep_vals,results2[:,4],yaxis=:log,xlabel = "Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="RN Infidelity")
Plots.plot!(plc_short,scalefactor.*ep_vals,results[:,3],yaxis=:log,xlabel = "Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="NF Guard Level Pop.",linestyle=:dash)
Plots.plot!(plc_short,scalefactor.*ep_vals,results2[:,3],xlims=((-20,20)),yaxis=:log,xlabel = "Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="RN Guard Level Pop.",linestyle=:dash,legend= :outerright)

plc2 = Plots.plot(scalefactor.*ep_vals,results[:,2].+results[:,3],yaxis=:log,xlabel ="Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="Total NF Objective")
Plots.plot!(plc2,scalefactor.*ep_vals,results2[:,2].+results2[:,3],yaxis=:log,xlabel ="Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="Total RN Objective",legend= :outerright)

plc2_short = Plots.plot(scalefactor.*ep_vals,results[:,2].+results[:,3],yaxis=:log,xlabel ="Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="Total NF Objective")
Plots.plot!(plc2_short,scalefactor.*ep_vals,results2[:,2].+results2[:,3],xlims=((-20,20)),yaxis=:log,xlabel ="Hamiltonian Perturbation [MHz]",ylabel = "Objective Function", lab="Total RN Objective",legend= :outerright)

# Save plots to file
Plots.savefig(plc,  "robust_comparison_"* @sprintf("%1.1e", ep_max/(2*pi)) * "_"*@sprintf("%3.1f",T)*"_T_"*@sprintf("%3d",maxIter)*"_iters_"*@sprintf("%d",nquad)*"_N_"*@sprintf("%d",D1)*"_D1_separate.png")
Plots.savefig(plc_short,  "robust_comparison_"* @sprintf("%1.1e", ep_max/(2*pi)) * "_"*@sprintf("%3.1f",T)*"_T_"*@sprintf("%3d",maxIter)*"_iters_"*@sprintf("%d",nquad)*"_N_"*@sprintf("%d",D1)*"_D1_separate_short.png")
Plots.savefig(plc2,  "robust_comparison_"* @sprintf("%1.1e", ep_max/(2*pi)) * "_"*@sprintf("%3.1f",T)*"_T_"*@sprintf("%3d",maxIter)*"_iters_"*@sprintf("%d",nquad)*"_N_"*@sprintf("%d",D1)*"_D1_total.png")
Plots.savefig(plc2_short,  "robust_comparison_"* @sprintf("%1.1e", ep_max/(2*pi)) * "_"*@sprintf("%3.1f",T)*"_T_"*@sprintf("%3d",maxIter)*"_iters_"*@sprintf("%d",nquad)*"_N_"*@sprintf("%d",D1)*"_D1_total_short.png")

p_RN_1,q_RN_1 = evalctrl_no_carrier(params, pcof, 1, D1)
p_RN_2,q_RN_2 = evalctrl_no_carrier(params, pcof, 2, D1)
pfunc1 = scalefactor .* p_RN_1
qfunc1 = scalefactor .* q_RN_1
pmax1 = maximum(abs.(pfunc2))
qmax1 = maximum(abs.(qfunc2))
pfunc2 = scalefactor .* p_RN_2
qfunc2 = scalefactor .* q_RN_2
pmax1 = maximum(abs.(pfunc2))
qmax1 = maximum(abs.(qfunc2))
pmax = maximum([pmax1,pmax2])
qmax = maximum([qmax1,qmax2])

titlestr = "Rotating frame RN ctrl " * " Max-p=" *@sprintf("%.3e", pmax) * " Max-q=" *@sprintf("%.3e", qmax) * " " * unitStr
pl_ctrl_RN = Plots.plot(td, pfunc1, lab=L"p_{1,1}(t)", title = titlestr, xlabel="Time [ns]",
                                  ylabel=unitStr, legend= :topleft, linewidth=1.5, legendfontsize=12)
# add in the control function for the anti-symmetric Hamiltonian
Plots.plot!(pl_ctrl_RN,td, qfunc1, lab=L"q_{1,1}(t)", linewidth=1.5, legendfontsize=12)
Plots.plot!(pl_ctrl_RN,td, pfunc2, lab=L"p_{1,2}(t)", linewidth=1.5, legendfontsize=12)
Plots.plot!(pl_ctrl_RN,td, qfunc2, lab=L"q_{1,2}(t)", linewidth=1.5, legendfontsize=12)
Plots.savefig(pl_ctrl_RN,  "robust_comparison_"* @sprintf("%1.1e", ep_max/(2*pi)) * "_"*@sprintf("%3.1f",T)*"_T_"*@sprintf("%3d",maxIter)*"_iters_"*@sprintf("%d",nquad)*"_N_"*@sprintf("%d",D1)*"D1_RN_ctrl.png")

# Save coefficients
save_pcof("robust_comparison_"* @sprintf("%1.1e", ep_max/(2*pi)) * "_"*@sprintf("%3.1f",T)*"_T_"*@sprintf("%3d",maxIter)*"_iters_"*@sprintf("%d",nquad)*"_N_"*@sprintf("%d",D1)*"D1_NF_pcof.jld2",pcof_old)
save_pcof("robust_comparison_"* @sprintf("%1.1e", ep_max/(2*pi)) * "_"*@sprintf("%3.1f",T)*"_T_"*@sprintf("%3d",maxIter)*"_iters_"*@sprintf("%d",nquad)*"_N_"*@sprintf("%d",D1)*"D1_RN_pcof.jld2",pcof)


text:
const GR_SUPPORTED_TYPES = Union{
    MIME"image/svg", MIME"image/svg+xml", MIME"image/png", MIME"image/jpeg",
    MIME"image/tiff", MIME"image/bmp", MIME"application/pdf",
    MIME"application/postscript", MIME"application/x-tex"
}

backend_showable(::GRBackend, ::GR_SUPPORTED_TYPES, scene::SceneLike) = true

function gr_save(io, scene, filetype)
    fp = tempname() * "." * filetype

    touch(fp)

    GR.beginprint(fp)
    gr_draw(scene)
    GR.endprint()

    write(io, read(fp))

    rm(fp)
end

function backend_show(::GRBackend, io::IO, ::MIME"image/png", scene::Scene)
    AbstractPlotting.update!(scene)
    gr_save(io, scene, "png")
end

function backend_show(::GRBackend, io::IO, ::MIME"image/jpeg", scene::Scene)
    AbstractPlotting.update!(scene)
    gr_save(io, scene, "jpeg")
end

function backend_show(::GRBackend, io::IO, ::MIME"image/bmp", scene::Scene)
    AbstractPlotting.update!(scene)
    gr_save(io, scene, "bmp")
end

function backend_show(::GRBackend, io::IO, ::MIME"image/tiff", scene::Scene)
    AbstractPlotting.update!(scene)
    gr_save(io, scene, "tiff")
end

function backend_show(::GRBackend, io::IO, ::Union{MIME"image/svg", MIME"image/svg+xml"}, scene::Scene)
    AbstractPlotting.update!(scene)
    gr_save(io, scene, "svg")
end

function backend_show(::GRBackend, io::IO, ::MIME"application/pdf", scene::Scene)
    AbstractPlotting.update!(scene)
    gr_save(io, scene, "pdf")
end

function backend_show(::GRBackend, io::IO, ::MIME"application/postscript", scene::Scene)
    AbstractPlotting.update!(scene)
    gr_save(io, scene, "eps")
end

function backend_show(::GRBackend, io::IO, ::MIME"application/x-tex", scene::Scene)
    AbstractPlotting.update!(scene)
    fp = tempname() * ".tex"
    withenv("GKS_WSTYPE" => "pgf", "GKS_FILEPATH" => fp) do
        GR.clearws()
        gr_draw(scene)
        GR.updatews()
    end

    write(io, read(fp))
end

function gr_record(f::Function, filename::String, scene::Scene, iter)
    ext = uppercase(splitext(filename)[2][2:end])
    @assert ext in ("GIF", "MOV", "MP4", "WEBM", "OGG") """
    Extension of file is incorrect!  Expected one of (\"GIF\", \"MOV\", \"MP4\", \"WEBM\", \"OGG\").
    Found $ext.
    """
    withenv("GKS_WSTYPE" => uppercase(ext), "GKS_FILEPATH" => filename) do
        for i in iter
            GR.clearws()
            f(i)
            gr_draw(scene)
        end
    end
end


text:
using Mmap

rmin = -2.50
rmax = 2.50
imin = -1.5
imax = 1.5

height = 200_000
width = 400_000

tStart=time()
cset= [ComplexF64(r,i) for r in LinRange(rmin ,rmax , width) for i in LinRange(imin, imax, height)]
tStop=time()
print("done. took ", tStop-tStart, " seconds\n");

# Create a file for mmapping
s = open("zset-$(width)-$(height).bin", "w+")
write(s, size(cset,1))
# Now write the data
write(s, cset)
close(s)

# Test by reading it back in
s = open("zset-$(width)-$(height).bin")   # default is read-only
n = read(s, Int)
println("size of zset is $n")
zs = Mmap.mmap(s, Array{ComplexF64,1}, n)
close(s)

max_stars_repo_name: rabarar/juliaSnippets

max_issues_repo_name: rabarar/juliaSnippets

max_forks_repo_name: rabarar/juliaSnippets


text:
include("testsetup.jl")

import InteractiveUtils

InteractiveUtils.versioninfo(stdout; verbose = true)

const START_TIME = time()

@show LoopVectorization.register_count()

@show RUN_SLOW_TESTS

@time @testset "LoopVectorization.jl" begin

    # @time Aqua.test_all(LoopVectorization)
    # @test isempty(detect_unbound_args(LoopVectorization))

    @time include("printmethods.jl")

    @time include("can_avx.jl")

    @time include("fallback.jl")

    @time include("utils.jl")

    @time include("arraywrappers.jl")

    @time include("check_empty.jl")

    @time include("loopinductvars.jl")

    @time include("shuffleloadstores.jl")

    @time include("zygote.jl")

    @time include("offsetarrays.jl")

    @time include("tensors.jl")

    @time include("map.jl")

    @time include("filter.jl")

    @time include("mapreduce.jl")

    @time include("ifelsemasks.jl")

    @time include("dot.jl")

    @time include("special.jl")

    @time include("gemv.jl")

    @time include("miscellaneous.jl")

    @time include("copy.jl")

    @time include("broadcast.jl")

    @time include("gemm.jl")

    @time include("threading.jl")

    @time include("tullio.jl")
end

const ELAPSED_MINUTES = (time() - START_TIME)/60
# @test ELAPSED_MINUTES < 180
@test ELAPSED_MINUTES < 240


text:
export restart, step_forward, get_state

function restart(roadway::Roadway, render::Bool)
    # Create scene
    scene = Scene()
    push!(scene,Vehicle(VehicleState(VecSE2(5.0 + rand()*10.0,randn()*0.5,0.0), roadway, 20.0),
                    VehicleDef(1, AgentClass.CAR, 4.826, 1.81)))
    push!(scene,Vehicle(VehicleState(VecSE2(0.,randn()*0.5,0.0), roadway, 10.0 + rand()*20.0),
                        VehicleDef(2, AgentClass.CAR, 4.826, 1.81)))

    # Specify models
    context = IntegratedContinuous(1/10,3) #1/framrate
    models = Dict{Int, DriverModel}()

    if render
    	model1 = LatLonSeparableDriver(context,
	                                ProportionalLaneTracker(σ=0.0, kp=5.0, kd=0.5),
	                                IntelligentDriverModel(σ=0.0, v_des=21.0))
	    model2 = LatLonSeparableDriver(context,
	                                ProportionalLaneTracker(σ=0.0, kp=5.0, kd=0.5),
	                                StaticLongitudinalDriver())
    else
	    model1 = LatLonSeparableDriver(context,
	                                ProportionalLaneTracker(σ=0.01, kp=5.0, kd=0.5),
	                                IntelligentDriverModel(σ=0.2, v_des=21.0))
	    model2 = LatLonSeparableDriver(context,
	                                ProportionalLaneTracker(σ=0.01, kp=5.0, kd=0.5),
	                                StaticLongitudinalDriver())
	end

    # Placeholder for actions
    return scene, model1, model2
end
function step_forward(scene::Scene, roadway::Roadway, model1::DriverModel, model2::DriverModel, action::Float64)

    observe!(model1, scene, roadway, 1)
    context = action_context(model1)
    scene[1].state = propagate(scene[1], rand(model1), context, roadway)

    model2.mlon.a = action
    observe!(model2, scene, roadway, 2)
    context = action_context(model2)
    scene[2].state = propagate(scene[2], rand(model2), context, roadway)

    return get_state(scene)
end
function get_state(scene::Scene)
	# Extract state values
	d = scene[1].state.posF.s - scene[2].state.posF.s - scene[1].def.length/2. - scene[2].def.length/2.
	r = scene[2].state.v - scene[1].state.v
    s = scene[1].state.v
    return d, r, s
end

max_stars_repo_name: yoshihikokuwahara/gail-driver

max_issues_repo_name: yoshihikokuwahara/gail-driver

max_forks_repo_name: yoshihikokuwahara/gail-driver


text:
"""
Analog to CircularBuffer, but with better performance
Fixed-size - no preallocation
New items are pushed to the back of the list, overwriting values in a circular fashion.
"""
mutable struct RingBuffer{T} <: AbstractVector{T}
    buffer::Vector{T}
    len::Int
    first::Int
    last::Int
    contentLen::Int
    # counter::Int
    function RingBuffer{T}(maxlen::Int) where T
        new{T}(Vector{T}(undef, maxlen), maxlen, 1, 0, 0) #, 0)
    end
end

Base.@propagate_inbounds function _buffer_index_checked(buf::RingBuffer, i::Int)
    @boundscheck if i < firstindex(buf) || i > lastindex(buf)
        throw(BoundsError(buf, i))
    end
    _buffer_index(buf, i)
end

@inline function _buffer_index(buf::RingBuffer, i::Int)
    #idx = mod1(buf.first + i - 1, buf.len) # idx = (buf.first + i - 2) % buf.len + 1
    idx = buf.first + i - 1
    # idx > buf.len ? idx - buf.len : idx
    ifelse(idx > buf.len, idx - buf.len, idx)
end

@inline Base.@propagate_inbounds function Base.getindex(buf::RingBuffer, i::Int)
    buf.buffer[_buffer_index_checked(buf, i)]
end

@inline Base.@propagate_inbounds function Base.setindex!(buf::RingBuffer, data, i::Int)
    buf.buffer[_buffer_index_checked(buf, i)] = data
    buf
end

# pop_back
@inline function Base.pop!(buf::RingBuffer)
    i = buf.last
    if buf.contentLen > 0
        buf.contentLen -= 1
        buf.last = _buffer_index(buf, buf.contentLen)
    else
        throw(ArgumentError("array must be non-empty"))
    end
    buf.buffer[i]
end
@inline function Base.pop!(buf::RingBuffer, Npop::Int)
    if Npop > buf.contentLen
        Npop = buf.contentLen
    end
    if buf.contentLen > 0
        buf.contentLen -= Npop
        buf.last = _buffer_index(buf, buf.contentLen)
        # buf.counter -= Npop
    end
    buf
end

# pop_front
@inline function Base.popfirst!(buf::RingBuffer)
    i = buf.first
    if buf.contentLen > 0
        buf.contentLen -= 1
        buf.first = ifelse(buf.first == buf.len, 1, buf.first + 1)
    else
        throw(ArgumentError("array must be non-empty"))
    end
    buf.buffer[i]
end
@inline function Base.popfirst!(buf::RingBuffer, Npop::Int)
    if Npop > buf.contentLen
        Npop = buf.contentLen
    end
    if buf.contentLen > 0
        buf.first = _buffer_index(buf, 1 + Npop)
        buf.contentLen -= Npop
    end
    buf
end

# push_back
@inline function Base.push!(buf::RingBuffer, data)
    buf.last = _buffer_index(buf, buf.contentLen + 1)
    buf.buffer[buf.last] = data
    if buf.contentLen == buf.len
        buf.first = _buffer_index(buf, 2) # (buf.first == buf.len ? 1 : buf.first + 1)
    else
        buf.contentLen += 1
    end
    # buf.counter += 1
    buf
end
# push_front
function Base.pushfirst!(buf::RingBuffer, data)
    buf.first = _buffer_index(buf, 0) # (buf.first == 1 ? buf.len : buf.first - 1)
    buf.buffer[buf.first] = data
    if buf.contentLen == buf.len
        buf.last = _buffer_index(buf, 0)
    else
        buf.contentLen += 1
    end
    buf
end

# push_back vector
function Base.append!(buf::RingBuffer, datavec::AbstractVector)
    Nadd = length(datavec)
    for i = 1:Nadd
        buf[i + buf.last] = datavec[i]
    end
    if buf.contentLen == buf.len
        buf.first = _buffer_index(buf, 1 + Nadd)
    else
        if buf.len > buf.contentLen + Nadd
            buf.contentLen += Nadd;
        else
            buf.first = _buffer_index(buf, 1 + buf.contentLen + Nadd - buf.len);
            buf.contentLen = buf.len;
        end
    end
    buf.last = _buffer_index(buf, buf.contentLen)
    # buf.counter += Nadd
    buf
end

Base.first(buf::RingBuffer) = buf.buffer[buf.first]
Base.last(buf::RingBuffer) = buf.buffer[buf.last]

Base.length(buf::RingBuffer) = buf.contentLen
Base.size(buf::RingBuffer) = (buf.contentLen,)
#Base.convert{T}(::Type{Array}, buf::RingBuffer{T}) = T[x for x in buf]
Base.convert(::Type{Array}, buf::RingBuffer{T}) where T = T[x for x in buf]
Base.isempty(buf::RingBuffer) = isempty(buf.buffer)

capacity(buf::RingBuffer) = buf.len # full capacity
isfull(buf::RingBuffer) = buf.contentLen == buf.len
# bounds(buf::RingBuffer) = (buf.first, buf.last) # bounds indexes

function Base.empty!(buf::RingBuffer)
    buf.contentLen = buf.last = 0
    buf.first = 1
    buf
end
# clear and resize
function reset!(buf::RingBuffer, len::Int = buf.len) # clear and resize
    empty!(buf)
    if buf.len != len
        buf.len = len
        buf.buffer = Vector{T}(undef, len)
    end
end


text:
using JuliaChem.JCModules

using Base.Threads
using LinearAlgebra
using Printf

#=
"""
	 index(a::Int64,b::Int64)
Summary
======
Triangular indexing determination.

Arguments
======
a = row index

b = column index
"""
=#
@inline function triangular_index(a::Int,b::Int)
  index = (a*(a-1)) >> 1 #bitwise divide by 2
  index += b
  return index
end

@inline function triangular_index(a::Int)
  return (a*(a-1)) >> 1
end

@inline function decompose(input::Int)
  #return ceil(Int,(-1.0+√(1+8*input))/2.0)
  return Base.fptosi(Int, Base.ceil_llvm((-1.0 +
    Base.Math.sqrt_llvm(float(1+8*input)))/2.0))
    #return ccall((:decompose, "/export/home/david/projects/Julia/JuliaChem.jl/src/eri/libjeri.so"),
  #  Int64, (Int64,), input)
end

function axial_normalization_factor(oei, ash, bsh, ncoord)
  ama = ash.am
  amb = bsh.am

  na = ash.nbas
  nb = bsh.nbas

  for icoord in 1:ncoord
    ab = (icoord-1)*na*nb
    for asize::Int64 in 0:(na-1), bsize::Int64 in 0:(nb-1)
      ab += 1

      anorm = axial_norm_fact[asize+1,ama]
      bnorm = axial_norm_fact[bsize+1,amb]

      abnorm = anorm*bnorm
      oei[ab] *= abnorm
    end
  end
end

function compute_dipole(mol::Molecule,
  basis::Basis, P::Matrix{Float64}, jeri_prop_engine)

  #== define initial variables ==#
  natoms = length(mol.atoms)
  dipole = zeros(Float64,(3,))

  ncoord = length(dipole)

  #== compute nuclear contribution to dipole ==#
  nuc_dipole = zeros(Float64,size(dipole))
  for atom in mol
    nuc_dipole .+= atom.atom_id .* atom.atom_center
  end

  #== compute electronic contribution to dipole ==#
  elec_dipole = zeros(Float64,size(dipole))
  for ash in 1:length(basis), bsh in 1:length(basis)
    abas = basis[ash].nbas
    bbas = basis[bsh].nbas

    apos = basis[ash].pos
    bpos = basis[bsh].pos

    elec_dipole_block_JERI = zeros(Float64,(abas*bbas*ncoord))

    JERI.compute_dipole_block(jeri_prop_engine, elec_dipole_block_JERI, ash,
      bsh, abas*bbas)

    axial_normalization_factor(elec_dipole_block_JERI, basis[ash], basis[bsh],
      ncoord)

    #println("Julia Shells: $ash, $bsh")
    #println("Julia Atoms: $iatom, $jatom")

    for icoord in 1:ncoord
      idx = 1

      #println("$ishlset, $icoord => $op, $shlset_idx")

      for ibas in 0:abas-1, jbas in 0:bbas-1
        aorb = apos + ibas
        borb = bpos + jbas

        #scale = iorb == jorb ? 2.0 : 1.0
        scale = 1.0
        elec_dipole[icoord] += scale * P[aorb,borb] *
          elec_dipole_block_JERI[abas*bbas*(icoord-1) + idx]

        idx += 1
      end
    end
  end

  #== compute total dipole ==#
  dipole .= 2.54174623 .* (nuc_dipole .- elec_dipole)
  return dipole
end

function eri_quartet_batch_size(max_am)
  return am_to_nbas_cart(max_am)^4
end


text:
#==============================================================================#
# Shield.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/shield-2016-06-02.normal.json
#==============================================================================#

__precompile__()

module Shield

using AWSCore

"""
    using AWSSDK.Shield.associate_drtlog_bucket
    associate_drtlog_bucket([::AWSConfig], arguments::Dict)
    associate_drtlog_bucket([::AWSConfig]; LogBucket=)

    using AWSCore.Services.shield
    shield([::AWSConfig], "AssociateDRTLogBucket", arguments::Dict)
    shield([::AWSConfig], "AssociateDRTLogBucket", LogBucket=)

# AssociateDRTLogBucket Operation

Authorizes the DDoS Response team (DRT) to access the specified Amazon S3 bucket containing your flow logs. You can associate up to 10 Amazon S3 buckets with your subscription.

To use the services of the DRT and make an `AssociateDRTLogBucket` request, you must be subscribed to the [Business Support plan](https://aws.amazon.com/premiumsupport/business-support/) or the [Enterprise Support plan](https://aws.amazon.com/premiumsupport/enterprise-support/).

# Arguments

## `LogBucket = ::String` -- *Required*
The Amazon S3 bucket that contains your flow logs.

# Returns

`AssociateDRTLogBucketResponse`

# Exceptions

`InternalErrorException`, `InvalidOperationException`, `NoAssociatedRoleException`, `LimitsExceededException`, `InvalidParameterException`, `AccessDeniedForDependencyException`, `OptimisticLockException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/AssociateDRTLogBucket)
"""
@inline associate_drtlog_bucket(aws::AWSConfig=default_aws_config(); args...) = associate_drtlog_bucket(aws, args)

@inline associate_drtlog_bucket(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "AssociateDRTLogBucket", args)

@inline associate_drtlog_bucket(args) = associate_drtlog_bucket(default_aws_config(), args)

"""
    using AWSSDK.Shield.associate_drtrole
    associate_drtrole([::AWSConfig], arguments::Dict)
    associate_drtrole([::AWSConfig]; RoleArn=)

    using AWSCore.Services.shield
    shield([::AWSConfig], "AssociateDRTRole", arguments::Dict)
    shield([::AWSConfig], "AssociateDRTRole", RoleArn=)

# AssociateDRTRole Operation

Authorizes the DDoS Response team (DRT), using the specified role, to access your AWS account to assist with DDoS attack mitigation during potential attacks. This enables the DRT to inspect your AWS WAF configuration and create or update AWS WAF rules and web ACLs.

You can associate only one `RoleArn` with your subscription. If you submit an `AssociateDRTRole` request for an account that already has an associated role, the new `RoleArn` will replace the existing `RoleArn`.

Prior to making the `AssociateDRTRole` request, you must attach the [AWSShieldDRTAccessPolicy](https://console.aws.amazon.com/iam/home?#/policies/arn:aws:iam::aws:policy/service-role/AWSShieldDRTAccessPolicy) managed policy to the role you will specify in the request. For more information see [Attaching and Detaching IAM Policies]( https://docs.aws.amazon.com/IAM/latest/UserGuide/access_policies_manage-attach-detach.html). The role must also trust the service principal `drt.shield.amazonaws.com`. For more information, see [IAM JSON Policy Elements: Principal](https://docs.aws.amazon.com/IAM/latest/UserGuide/reference_policies_elements_principal.html).

The DRT will have access only to your AWS WAF and Shield resources. By submitting this request, you authorize the DRT to inspect your AWS WAF and Shield configuration and create and update AWS WAF rules and web ACLs on your behalf. The DRT takes these actions only if explicitly authorized by you.

You must have the `iam:PassRole` permission to make an `AssociateDRTRole` request. For more information, see [Granting a User Permissions to Pass a Role to an AWS Service](https://docs.aws.amazon.com/IAM/latest/UserGuide/id_roles_use_passrole.html).

To use the services of the DRT and make an `AssociateDRTRole` request, you must be subscribed to the [Business Support plan](https://aws.amazon.com/premiumsupport/business-support/) or the [Enterprise Support plan](https://aws.amazon.com/premiumsupport/enterprise-support/).

# Arguments

## `RoleArn = ::String` -- *Required*
The Amazon Resource Name (ARN) of the role the DRT will use to access your AWS account.

Prior to making the `AssociateDRTRole` request, you must attach the [AWSShieldDRTAccessPolicy](https://console.aws.amazon.com/iam/home?#/policies/arn:aws:iam::aws:policy/service-role/AWSShieldDRTAccessPolicy) managed policy to this role. For more information see [Attaching and Detaching IAM Policies]( https://docs.aws.amazon.com/IAM/latest/UserGuide/access_policies_manage-attach-detach.html).

# Returns

`AssociateDRTRoleResponse`

# Exceptions

`InternalErrorException`, `InvalidOperationException`, `InvalidParameterException`, `AccessDeniedForDependencyException`, `OptimisticLockException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/AssociateDRTRole)
"""
@inline associate_drtrole(aws::AWSConfig=default_aws_config(); args...) = associate_drtrole(aws, args)

@inline associate_drtrole(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "AssociateDRTRole", args)

@inline associate_drtrole(args) = associate_drtrole(default_aws_config(), args)

"""
    using AWSSDK.Shield.create_protection
    create_protection([::AWSConfig], arguments::Dict)
    create_protection([::AWSConfig]; Name=, ResourceArn=)

    using AWSCore.Services.shield
    shield([::AWSConfig], "CreateProtection", arguments::Dict)
    shield([::AWSConfig], "CreateProtection", Name=, ResourceArn=)

# CreateProtection Operation

Enables AWS Shield Advanced for a specific AWS resource. The resource can be an Amazon CloudFront distribution, Elastic Load Balancing load balancer, Elastic IP Address, or an Amazon Route 53 hosted zone.

You can add protection to only a single resource with each CreateProtection request. If you want to add protection to multiple resources at once, use the [AWS WAF console](https://console.aws.amazon.com/waf/). For more information see [Getting Started with AWS Shield Advanced](https://docs.aws.amazon.com/waf/latest/developerguide/getting-started-ddos.html) and [Add AWS Shield Advanced Protection to more AWS Resources](https://docs.aws.amazon.com/waf/latest/developerguide/configure-new-protection.html).

# Arguments

## `Name = ::String` -- *Required*
Friendly name for the `Protection` you are creating.

## `ResourceArn = ::String` -- *Required*
The ARN (Amazon Resource Name) of the resource to be protected.

The ARN should be in one of the following formats:

*   For an Application Load Balancer: `arn:aws:elasticloadbalancing:*region*:*account-id*:loadbalancer/app/*load-balancer-name*/*load-balancer-id*`

*   For an Elastic Load Balancer (Classic Load Balancer): `arn:aws:elasticloadbalancing:*region*:*account-id*:loadbalancer/*load-balancer-name*`

*   For AWS CloudFront distribution: `arn:aws:cloudfront::*account-id*:distribution/*distribution-id*`

*   For Amazon Route 53: `arn:aws:route53::*account-id*:hostedzone/*hosted-zone-id*`

*   For an Elastic IP address: `arn:aws:ec2:*region*:*account-id*:eip-allocation/*allocation-id*`

# Returns

`CreateProtectionResponse`

# Exceptions

`InternalErrorException`, `InvalidResourceException`, `InvalidOperationException`, `LimitsExceededException`, `ResourceAlreadyExistsException`, `OptimisticLockException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/CreateProtection)
"""
@inline create_protection(aws::AWSConfig=default_aws_config(); args...) = create_protection(aws, args)

@inline create_protection(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "CreateProtection", args)

@inline create_protection(args) = create_protection(default_aws_config(), args)

"""
    using AWSSDK.Shield.create_subscription
    create_subscription([::AWSConfig], arguments::Dict)
    create_subscription([::AWSConfig]; )

    using AWSCore.Services.shield
    shield([::AWSConfig], "CreateSubscription", arguments::Dict)
    shield([::AWSConfig], "CreateSubscription", )

# CreateSubscription Operation

Activates AWS Shield Advanced for an account.

As part of this request you can specify `EmergencySettings` that automaticaly grant the DDoS response team (DRT) needed permissions to assist you during a suspected DDoS attack. For more information see [Authorize the DDoS Response Team to Create Rules and Web ACLs on Your Behalf](https://docs.aws.amazon.com/waf/latest/developerguide/authorize-DRT.html).

When you initally create a subscription, your subscription is set to be automatically renewed at the end of the existing subscription period. You can change this by submitting an `UpdateSubscription` request.

# Arguments

# Returns

`CreateSubscriptionResponse`

# Exceptions

`InternalErrorException` or `ResourceAlreadyExistsException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/CreateSubscription)
"""
@inline create_subscription(aws::AWSConfig=default_aws_config(); args...) = create_subscription(aws, args)

@inline create_subscription(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "CreateSubscription", args)

@inline create_subscription(args) = create_subscription(default_aws_config(), args)

"""
    using AWSSDK.Shield.delete_protection
    delete_protection([::AWSConfig], arguments::Dict)
    delete_protection([::AWSConfig]; ProtectionId=)

    using AWSCore.Services.shield
    shield([::AWSConfig], "DeleteProtection", arguments::Dict)
    shield([::AWSConfig], "DeleteProtection", ProtectionId=)

# DeleteProtection Operation

Deletes an AWS Shield Advanced [Protection](@ref).

# Arguments

## `ProtectionId = ::String` -- *Required*
The unique identifier (ID) for the [Protection](@ref) object to be deleted.

# Returns

`DeleteProtectionResponse`

# Exceptions

`InternalErrorException`, `ResourceNotFoundException` or `OptimisticLockException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DeleteProtection)
"""
@inline delete_protection(aws::AWSConfig=default_aws_config(); args...) = delete_protection(aws, args)

@inline delete_protection(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DeleteProtection", args)

@inline delete_protection(args) = delete_protection(default_aws_config(), args)

"""
    using AWSSDK.Shield.delete_subscription
    delete_subscription([::AWSConfig], arguments::Dict)
    delete_subscription([::AWSConfig]; )

    using AWSCore.Services.shield
    shield([::AWSConfig], "DeleteSubscription", arguments::Dict)
    shield([::AWSConfig], "DeleteSubscription", )

# DeleteSubscription Operation

Removes AWS Shield Advanced from an account. AWS Shield Advanced requires a 1-year subscription commitment. You cannot delete a subscription prior to the completion of that commitment.

# Arguments

# Returns

`DeleteSubscriptionResponse`

# Exceptions

`InternalErrorException`, `LockedSubscriptionException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DeleteSubscription)
"""
@inline delete_subscription(aws::AWSConfig=default_aws_config(); args...) = delete_subscription(aws, args)

@inline delete_subscription(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DeleteSubscription", args)

@inline delete_subscription(args) = delete_subscription(default_aws_config(), args)

"""
    using AWSSDK.Shield.describe_attack
    describe_attack([::AWSConfig], arguments::Dict)
    describe_attack([::AWSConfig]; AttackId=)

    using AWSCore.Services.shield
    shield([::AWSConfig], "DescribeAttack", arguments::Dict)
    shield([::AWSConfig], "DescribeAttack", AttackId=)

# DescribeAttack Operation

Describes the details of a DDoS attack.

# Arguments

## `AttackId = ::String` -- *Required*
The unique identifier (ID) for the attack that to be described.

# Returns

`DescribeAttackResponse`

# Exceptions

`InternalErrorException` or `InvalidParameterException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DescribeAttack)
"""
@inline describe_attack(aws::AWSConfig=default_aws_config(); args...) = describe_attack(aws, args)

@inline describe_attack(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DescribeAttack", args)

@inline describe_attack(args) = describe_attack(default_aws_config(), args)

"""
    using AWSSDK.Shield.describe_drtaccess
    describe_drtaccess([::AWSConfig], arguments::Dict)
    describe_drtaccess([::AWSConfig]; )

    using AWSCore.Services.shield
    shield([::AWSConfig], "DescribeDRTAccess", arguments::Dict)
    shield([::AWSConfig], "DescribeDRTAccess", )

# DescribeDRTAccess Operation

Returns the current role and list of Amazon S3 log buckets used by the DDoS Response team (DRT) to access your AWS account while assisting with attack mitigation.

# Arguments

# Returns

`DescribeDRTAccessResponse`

# Exceptions

`InternalErrorException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DescribeDRTAccess)
"""
@inline describe_drtaccess(aws::AWSConfig=default_aws_config(); args...) = describe_drtaccess(aws, args)

@inline describe_drtaccess(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DescribeDRTAccess", args)

@inline describe_drtaccess(args) = describe_drtaccess(default_aws_config(), args)

"""
    using AWSSDK.Shield.describe_emergency_contact_settings
    describe_emergency_contact_settings([::AWSConfig], arguments::Dict)
    describe_emergency_contact_settings([::AWSConfig]; )

    using AWSCore.Services.shield
    shield([::AWSConfig], "DescribeEmergencyContactSettings", arguments::Dict)
    shield([::AWSConfig], "DescribeEmergencyContactSettings", )

# DescribeEmergencyContactSettings Operation

Lists the email addresses that the DRT can use to contact you during a suspected attack.

# Arguments

# Returns

`DescribeEmergencyContactSettingsResponse`

# Exceptions

`InternalErrorException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DescribeEmergencyContactSettings)
"""
@inline describe_emergency_contact_settings(aws::AWSConfig=default_aws_config(); args...) = describe_emergency_contact_settings(aws, args)

@inline describe_emergency_contact_settings(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DescribeEmergencyContactSettings", args)

@inline describe_emergency_contact_settings(args) = describe_emergency_contact_settings(default_aws_config(), args)

"""
    using AWSSDK.Shield.describe_protection
    describe_protection([::AWSConfig], arguments::Dict)
    describe_protection([::AWSConfig]; ProtectionId=)

    using AWSCore.Services.shield
    shield([::AWSConfig], "DescribeProtection", arguments::Dict)
    shield([::AWSConfig], "DescribeProtection", ProtectionId=)

# DescribeProtection Operation

Lists the details of a [Protection](@ref) object.

# Arguments

## `ProtectionId = ::String` -- *Required*
The unique identifier (ID) for the [Protection](@ref) object that is described.

# Returns

`DescribeProtectionResponse`

# Exceptions

`InternalErrorException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DescribeProtection)
"""
@inline describe_protection(aws::AWSConfig=default_aws_config(); args...) = describe_protection(aws, args)

@inline describe_protection(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DescribeProtection", args)

@inline describe_protection(args) = describe_protection(default_aws_config(), args)

"""
    using AWSSDK.Shield.describe_subscription
    describe_subscription([::AWSConfig], arguments::Dict)
    describe_subscription([::AWSConfig]; )

    using AWSCore.Services.shield
    shield([::AWSConfig], "DescribeSubscription", arguments::Dict)
    shield([::AWSConfig], "DescribeSubscription", )

# DescribeSubscription Operation

Provides details about the AWS Shield Advanced subscription for an account.

# Arguments

# Returns

`DescribeSubscriptionResponse`

# Exceptions

`InternalErrorException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DescribeSubscription)
"""
@inline describe_subscription(aws::AWSConfig=default_aws_config(); args...) = describe_subscription(aws, args)

@inline describe_subscription(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DescribeSubscription", args)

@inline describe_subscription(args) = describe_subscription(default_aws_config(), args)

"""
    using AWSSDK.Shield.disassociate_drtlog_bucket
    disassociate_drtlog_bucket([::AWSConfig], arguments::Dict)
    disassociate_drtlog_bucket([::AWSConfig]; LogBucket=)

    using AWSCore.Services.shield
    shield([::AWSConfig], "DisassociateDRTLogBucket", arguments::Dict)
    shield([::AWSConfig], "DisassociateDRTLogBucket", LogBucket=)

# DisassociateDRTLogBucket Operation

Removes the DDoS Response team's (DRT) access to the specified Amazon S3 bucket containing your flow logs.

To make a `DisassociateDRTLogBucket` request, you must be subscribed to the [Business Support plan](https://aws.amazon.com/premiumsupport/business-support/) or the [Enterprise Support plan](https://aws.amazon.com/premiumsupport/enterprise-support/). However, if you are not subscribed to one of these support plans, but had been previously and had granted the DRT access to your account, you can submit a `DisassociateDRTLogBucket` request to remove this access.

# Arguments

## `LogBucket = ::String` -- *Required*
The Amazon S3 bucket that contains your flow logs.

# Returns

`DisassociateDRTLogBucketResponse`

# Exceptions

`InternalErrorException`, `InvalidOperationException`, `NoAssociatedRoleException`, `AccessDeniedForDependencyException`, `OptimisticLockException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DisassociateDRTLogBucket)
"""
@inline disassociate_drtlog_bucket(aws::AWSConfig=default_aws_config(); args...) = disassociate_drtlog_bucket(aws, args)

@inline disassociate_drtlog_bucket(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DisassociateDRTLogBucket", args)

@inline disassociate_drtlog_bucket(args) = disassociate_drtlog_bucket(default_aws_config(), args)

"""
    using AWSSDK.Shield.disassociate_drtrole
    disassociate_drtrole([::AWSConfig], arguments::Dict)
    disassociate_drtrole([::AWSConfig]; )

    using AWSCore.Services.shield
    shield([::AWSConfig], "DisassociateDRTRole", arguments::Dict)
    shield([::AWSConfig], "DisassociateDRTRole", )

# DisassociateDRTRole Operation

Removes the DDoS Response team's (DRT) access to your AWS account.

To make a `DisassociateDRTRole` request, you must be subscribed to the [Business Support plan](https://aws.amazon.com/premiumsupport/business-support/) or the [Enterprise Support plan](https://aws.amazon.com/premiumsupport/enterprise-support/). However, if you are not subscribed to one of these support plans, but had been previously and had granted the DRT access to your account, you can submit a `DisassociateDRTRole` request to remove this access.

# Arguments

# Returns

`DisassociateDRTRoleResponse`

# Exceptions

`InternalErrorException`, `InvalidOperationException`, `OptimisticLockException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/DisassociateDRTRole)
"""
@inline disassociate_drtrole(aws::AWSConfig=default_aws_config(); args...) = disassociate_drtrole(aws, args)

@inline disassociate_drtrole(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "DisassociateDRTRole", args)

@inline disassociate_drtrole(args) = disassociate_drtrole(default_aws_config(), args)

"""
    using AWSSDK.Shield.get_subscription_state
    get_subscription_state([::AWSConfig], arguments::Dict)
    get_subscription_state([::AWSConfig]; )

    using AWSCore.Services.shield
    shield([::AWSConfig], "GetSubscriptionState", arguments::Dict)
    shield([::AWSConfig], "GetSubscriptionState", )

# GetSubscriptionState Operation

Returns the `SubscriptionState`, either `Active` or `Inactive`.

# Arguments

# Returns

`GetSubscriptionStateResponse`

# Exceptions

`InternalErrorException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/GetSubscriptionState)
"""
@inline get_subscription_state(aws::AWSConfig=default_aws_config(); args...) = get_subscription_state(aws, args)

@inline get_subscription_state(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "GetSubscriptionState", args)

@inline get_subscription_state(args) = get_subscription_state(default_aws_config(), args)

"""
    using AWSSDK.Shield.list_attacks
    list_attacks([::AWSConfig], arguments::Dict)
    list_attacks([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.shield
    shield([::AWSConfig], "ListAttacks", arguments::Dict)
    shield([::AWSConfig], "ListAttacks", <keyword arguments>)

# ListAttacks Operation

Returns all ongoing DDoS attacks or all DDoS attacks during a specified time period.

# Arguments

## `ResourceArns = [::String, ...]`
The ARN (Amazon Resource Name) of the resource that was attacked. If this is left blank, all applicable resources for this account will be included.

## `StartTime = [ ... ]`
The start of the time period for the attacks. This is a `timestamp` type. The sample request above indicates a `number` type because the default used by WAF is Unix time in seconds. However any valid [timestamp format](http://docs.aws.amazon.com/cli/latest/userguide/cli-using-param.html#parameter-types) is allowed.
```
 StartTime = [
        "FromInclusive" =>  timestamp,
        "ToExclusive" =>  timestamp
    ]
```

## `EndTime = [ ... ]`
The end of the time period for the attacks. This is a `timestamp` type. The sample request above indicates a `number` type because the default used by WAF is Unix time in seconds. However any valid [timestamp format](http://docs.aws.amazon.com/cli/latest/userguide/cli-using-param.html#parameter-types) is allowed.
```
 EndTime = [
        "FromInclusive" =>  timestamp,
        "ToExclusive" =>  timestamp
    ]
```

## `NextToken = ::String`
The `ListAttacksRequest.NextMarker` value from a previous call to `ListAttacksRequest`. Pass null if this is the first call.

## `MaxResults = ::Int`
The maximum number of [AttackSummary](@ref) objects to be returned. If this is left blank, the first 20 results will be returned.

# Returns

`ListAttacksResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException` or `InvalidOperationException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/ListAttacks)
"""
@inline list_attacks(aws::AWSConfig=default_aws_config(); args...) = list_attacks(aws, args)

@inline list_attacks(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "ListAttacks", args)

@inline list_attacks(args) = list_attacks(default_aws_config(), args)

"""
    using AWSSDK.Shield.list_protections
    list_protections([::AWSConfig], arguments::Dict)
    list_protections([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.shield
    shield([::AWSConfig], "ListProtections", arguments::Dict)
    shield([::AWSConfig], "ListProtections", <keyword arguments>)

# ListProtections Operation

Lists all [Protection](@ref) objects for the account.

# Arguments

## `NextToken = ::String`
The `ListProtectionsRequest.NextToken` value from a previous call to `ListProtections`. Pass null if this is the first call.

## `MaxResults = ::Int`
The maximum number of [Protection](@ref) objects to be returned. If this is left blank the first 20 results will be returned.

# Returns

`ListProtectionsResponse`

# Exceptions

`InternalErrorException`, `ResourceNotFoundException` or `InvalidPaginationTokenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/ListProtections)
"""
@inline list_protections(aws::AWSConfig=default_aws_config(); args...) = list_protections(aws, args)

@inline list_protections(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "ListProtections", args)

@inline list_protections(args) = list_protections(default_aws_config(), args)

"""
    using AWSSDK.Shield.update_emergency_contact_settings
    update_emergency_contact_settings([::AWSConfig], arguments::Dict)
    update_emergency_contact_settings([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.shield
    shield([::AWSConfig], "UpdateEmergencyContactSettings", arguments::Dict)
    shield([::AWSConfig], "UpdateEmergencyContactSettings", <keyword arguments>)

# UpdateEmergencyContactSettings Operation

Updates the details of the list of email addresses that the DRT can use to contact you during a suspected attack.

# Arguments

## `EmergencyContactList = [["EmailAddress" => <required> ::String], ...]`
A list of email addresses that the DRT can use to contact you during a suspected attack.

# Returns

`UpdateEmergencyContactSettingsResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `OptimisticLockException` or `ResourceNotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/UpdateEmergencyContactSettings)
"""
@inline update_emergency_contact_settings(aws::AWSConfig=default_aws_config(); args...) = update_emergency_contact_settings(aws, args)

@inline update_emergency_contact_settings(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "UpdateEmergencyContactSettings", args)

@inline update_emergency_contact_settings(args) = update_emergency_contact_settings(default_aws_config(), args)

"""
    using AWSSDK.Shield.update_subscription
    update_subscription([::AWSConfig], arguments::Dict)
    update_subscription([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.shield
    shield([::AWSConfig], "UpdateSubscription", arguments::Dict)
    shield([::AWSConfig], "UpdateSubscription", <keyword arguments>)

# UpdateSubscription Operation

Updates the details of an existing subscription. Only enter values for parameters you want to change. Empty parameters are not updated.

# Arguments

## `AutoRenew = "ENABLED" or "DISABLED"`
When you initally create a subscription, `AutoRenew` is set to `ENABLED`. If `ENABLED`, the subscription will be automatically renewed at the end of the existing subscription period. You can change this by submitting an `UpdateSubscription` request. If the `UpdateSubscription` request does not included a value for `AutoRenew`, the existing value for `AutoRenew` remains unchanged.

# Returns

`UpdateSubscriptionResponse`

# Exceptions

`InternalErrorException`, `LockedSubscriptionException`, `ResourceNotFoundException`, `InvalidParameterException` or `OptimisticLockException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/shield-2016-06-02/UpdateSubscription)
"""
@inline update_subscription(aws::AWSConfig=default_aws_config(); args...) = update_subscription(aws, args)

@inline update_subscription(aws::AWSConfig, args) = AWSCore.Services.shield(aws, "UpdateSubscription", args)

@inline update_subscription(args) = update_subscription(default_aws_config(), args)

end # module Shield

#==============================================================================#
# End of file
#==============================================================================#


text:
#TODO Perhaps we also need stack_trace and nested
struct ClickHouseServerException <: Exception
    code::Int
    name::String
    message::String
end


text:
@testset "IPNewton Unconstrained" begin
    method = IPNewton()
    run_optim_tests_constrained(method; show_name=debug_printing, show_res=false, show_itcalls=debug_printing)
end