text:
struct BackTracking2{O, TF, TI}
    c_1::TF
    ρ_hi::TF
    ρ_lo::TF
    iterations::TI
    maxstep::TF
end
function BackTracking2(::Val{O} = Val(3); c_1::TF = 1e-4, ρ_hi::TF = 0.5, ρ_lo::TF = 0.1, iterations::TI = 1_000, maxstep::TF = Inf) where {O,TF,TI}
    BackTracking2{O,TF,TI}(c_1, ρ_hi, ρ_lo, iterations, maxstep)
end
function BackTracking2{TF}(::Val{O} = Val(3); c_1 = TF(1e-4), ρ_hi = TF(0.5), ρ_lo = TF(0.1), iterations::TI = 1_000, maxstep = TF(Inf)) where {O,TF,TI}
    BackTracking2{O,TF,TI}(c_1, ρ_hi, ρ_lo, iterations, maxstep)
end

# struct StaticOptimizationResults{T, Tf}
#     # initial_x::Tx
#     # minimizer::Tx
#     minimum::Tf
#     iterations::Int
#     g_tol::T
#     f_calls::Int
#     g_calls::Int
#     g_converged::Bool
#     # h::Th
# end
# function Base.show(io::IO, r::StaticOptimizationResults)
#     @printf io "Results of Static Optimization Algorithm\n"
#     @printf io " * Minimizer: [%s]\n" join(r.minimizer, ",")
#     @printf io " * Minimum: [%s]\n" join(r.minimum, ",")
#     @printf io " * Hf(x): [%s]\n" join(r.h, ",")
#     @printf io " * Number of iterations: [%s]\n" join(r.iterations, ",")
#     @printf io " * Number of function calls: [%s]\n" join(r.f_calls, ",")
#     @printf io " * Number of gradient calls: [%s]\n" join(r.g_calls, ",")
#     @printf io " * Converged: [%s]\n" join(r.g_converged, ",")
#     nothing
# end

abstract type AbstractBFGSState{P,T,L,LT} end

mutable struct BFGSState{P,T,L,LT} <: AbstractBFGSState{P,T,L,LT}
    invH::ConstantFixedSizeMatrix{P,P,T,L,LT}
    x_old::ConstantFixedSizeVector{P,T,L,L}
    x_new::ConstantFixedSizeVector{P,T,L,L}
    ∇_old::ConstantFixedSizeVector{P,T,L,L}
    # ∇_new::SizedSIMDVector{P,T,L}
    δ∇::ConstantFixedSizeVector{P,T,L,L}
    u::ConstantFixedSizeVector{P,T,L,L}
    s::ConstantFixedSizeVector{P,T,L,L}
    function BFGSState{P,T,L,LT}(::UndefInitializer) where {P,T,L,LT}
        new{P,T,L,LT}()
    end
    @generated function BFGSState{P}(::UndefInitializer) where {P}
        L = PaddedMatrices.calc_padding(P, Float64)
        :(BFGSState{$P,Float64,$L,$(P*L)}(undef))
    end
    @generated function BFGSState{P,T}(::UndefInitializer) where {P,T}
        L = PaddedMatrices.calc_padding(P, T)
        :(BFGSState{$P,$T,$L,$(P*L)}(undef))
    end
end
BFGSState(::Val{P}, ::Type{T} = Float64) where {P,T} = BFGSState{P,T}(undef)
@inline Base.pointer(s::BFGSState{P,T})  where {P,T} = Base.unsafe_convert(Ptr{T}, pointer_from_objref(s))

"""
This type exists primarily to be the field of another mutable struct, so that you can get a pointer to this object.
"""
struct ConstantBFGSState{P,T,L,LT} <: AbstractBFGSState{P,T,L,LT}
    invH::ConstantFixedSizeMatrix{P,P,T,L,LT}
    x_old::ConstantFixedSizeVector{P,T,L,L}
    x_new::ConstantFixedSizeVector{P,T,L,L}
    ∇_old::ConstantFixedSizeVector{P,T,L,L}
    # ∇_new::SizedSIMDVector{P,T,L}
    δ∇::ConstantFixedSizeVector{P,T,L,L}
    u::ConstantFixedSizeVector{P,T,L,L}
    s::ConstantFixedSizeVector{P,T,L,L}
end
struct PtrBFGSState{P,T,L,LT} <: AbstractBFGSState{P,T,L,LT}
    ptr::Ptr{T}
end
@inline Base.pointer(state::PtrBFGSState) = state.ptr

@inline ref_invH(s::AbstractBFGSState{P,T,L,LT}) where {P,T,L,LT} = PtrMatrix{P,P,T,L,LT,false}(pointer(s))
@inline ref_x_old(s::AbstractBFGSState{P,T,L,LT}) where {P,T,L,LT} = PtrVector{P,T,L,L,false}(pointer(s) + LT*sizeof(T))
@inline ref_x_new(s::AbstractBFGSState{P,T,L,LT}) where {P,T,L,LT} = PtrVector{P,T,L,L,false}(pointer(s) + (LT+L)*sizeof(T))
@inline ref_∇_old(s::AbstractBFGSState{P,T,L,LT}) where {P,T,L,LT} = PtrVector{P,T,L,L,false}(pointer(s) + (LT+2L)*sizeof(T))
@inline ref_δ∇(s::AbstractBFGSState{P,T,L,LT}) where {P,T,L,LT} = PtrVector{P,T,L,L,false}(pointer(s) + (LT+3L)*sizeof(T))
@inline ref_u(s::AbstractBFGSState{P,T,L,LT}) where {P,T,L,LT} = PtrVector{P,T,L,L,false}(pointer(s) + (LT+4L)*sizeof(T))
@inline ref_s(s::AbstractBFGSState{P,T,L,LT}) where {P,T,L,LT} = PtrVector{P,T,L,L,false}(pointer(s) + (LT+5L)*sizeof(T))

function initial_invH!(state::AbstractBFGSState{P,T}) where {P,T}
    invH = ref_invH(state)
    fill!(invH, zero(T))
    @inbounds for p = 1:P
        invH[p,p] = one(T)
    end
end
@inline optimum(s::AbstractBFGSState{P,T,L,LT}) where {P,T,L,LT} = PtrVector{P,T,L,L,false}(pointer(s) + LT*sizeof(T))

@noinline linesearch_failure(iterations) = error("Linesearch failed to converge, reached maximum iterations $(iterations).")

"""
Optimum value is stored in state.x_old.

"""
function optimize_light!(state, obj, x::AbstractFixedSizeVector{P,T,L}, ls::BackTracking2{order}, tol = 1e-8) where {P,T,L,order}
    # res = DiffResults.GradientResult(x)
    # ls = BackTracking()
    # order = ordernum(bto)
    # copyto!(state.xinit, x)
    # copyto!(state.x_new, x)
    # xinit = copy(x)
    # x_new = copy(x)
    x_old = ref_x_old(state)
    ∇_old = ref_∇_old(state)
    invH = ref_invH(state)
    δ∇ = ref_δ∇(state)
    s = ref_s(state)
    u = ref_u(state)
    x_new = ref_x_new(state)
    copyto!(x_old, x)
    initial_invH!(state)
    # @show x_old
    # hx = SMatrix{P,P,T}(I)
    # if !(hguess isa Nothing)
    #     hx = hguess * hx
    # end
    # hold = copy(hx)
    # jold = copy(x); s = copy(x)
    # @unpack c_1, ρ_hi, ρ_lo, iterations = ls
    c_1, ρ_hi, ρ_lo, iterations = ls.c_1, ls.ρ_hi, ls.ρ_lo, ls.iterations
    iterfinitemax = round(Int,-log2(eps(eltype(x))))
    sqrttol = sqrt(eps(T))
    α_0 = one(T)
    N = 200
    f_calls = 0
    g_calls = 0
    for n = 1:N
        # res = ForwardDiff.gradient!(res, f, x); f_calls +=1; g_calls +=1; # Obtain gradient
        # ∇ = gradient!(d, x_old)
        # # ϕ_0 = DiffResults.value(res)
        # ϕ_0 = value(d)

        ϕ_0 = fdf(obj, x_old); f_calls +=1; g_calls +=1;
        ∇ = gradient(obj)
        # @show ∇, x_old
        # @show maximum(abs, ∇), tol
        isfinite(ϕ_0) || return T(NaN)
        maximum(abs, ∇) < tol && return ϕ_0
        if n > 1 # update hessian
            # y = jx - jold
            # hx = norm(y) < eps(eltype(x)) ? hx : hx + y*y' / (y'*s) - (hx*(s*s')*hx)/(s'*hx*s)

            # invH = invH + c1 * (s * s') - c2 * (u * s' + s * u')
            # @show L
            # @show length(δ∇), length(∇), length(∇_old)
            # @show full_length(δ∇), full_length(∇), full_length(∇_old)
            # @show typeof(δ∇), typeof(∇), typeof(∇_old)
            @inbounds @simd for i ∈ 1:L
                δ∇[i] =  ∇[i] - ∇_old[i]
            end
            # SIMDArrays.vsub!(δ∇, ∇, ∇_old)
            # Update the inverse Hessian approximation using Sherman-Morrison
            dx_dg = real(dot(s, δ∇))
            # dx_dg == 0.0 && return true # force stop
            mul!(u, invH, δ∇)
            c2 = one(T) / dx_dg
            c1 = fma(real(dot(δ∇, u)), c2*c2, c2)
            BFGS_update!(invH, s, u, c1, c2)
        end
        mul!(s, invH, ∇)
        # SIMDArrays.scale!(s, -one(eltype(s)))
        # SIMDArrays.reflect!(s)
        @inbounds @simd for i ∈ 1:L
            s[i] *= -1
        end
        # s = -hx\jx # Obtain direction
        dϕ_0 = dot(∇, s)

        if dϕ_0 >= zero(T) # If bad, reset search direction
            initial_invH!(state)
            # hx = hold
            # s = -jx
            # @show L, s, ∇
            # @show length(s)
            # @show length(∇)
            # @inbounds @simd for i ∈ 1:L
            #     s[i] = -∇[i]
            # end
            # SIMDArrays.reflect!(s, ∇)
            @inbounds @simd for i ∈ 1:L
                s[i] = - ∇[i]
            end
            dϕ_0 = dot(∇, s)
        end
        #### Perform line search

        # Count the total number of iterations
        iteration = 0
        ϕx_0, ϕx_1 = ϕ_0, ϕ_0
        α_1, α_2 = α_0, α_0
        @inbounds @simd for i ∈ 1:L
            x_new[i] = x_old[i] + α_1*s[i]
        end
        # SIMDArrays.vadd!(x_new, x_old, α_1, s)
        # ϕx_1 = f(x + α_1*s); f_calls += 1;
        ϕx_1 = f(obj, x_new); f_calls += 1;

        # Hard-coded backtrack until we find a finite function value
        iterfinite = 0
        while !isfinite(ϕx_1) && iterfinite < iterfinitemax
            iterfinite += 1
            α_1 = α_2
            α_2 = T(0.5)*α_1
            @inbounds @simd for i ∈ 1:L
                x_new[i] = x_old[i] + α_2*s[i]
            end
            # SIMDArrays.vadd!(x_new, x_old, α_2, s)
            # ϕx_1 = f(x + α_2*s); f_calls += 1;
            ϕx_1 = f(obj, x_new); f_calls += 1;
        end

        # Backtrack until we satisfy sufficient decrease condition
        while ϕx_1 > ϕ_0 + c_1 * α_2 * dϕ_0
            # Increment the number of steps we've had to perform
            iteration += 1

            # Ensure termination
            iteration > iterations && linesearch_failure(iterations)

            # Shrink proposed step-size:
            @fastmath if order == 2 || iteration == 1
                # backtracking via quadratic interpolation:
                # This interpolates the available data
                #    f(0), f'(0), f(α)
                # with a quadractic which is then minimised; this comes with a
                # guaranteed backtracking factor 0.5 * (1-c_1)^{-1} which is < 1
                # provided that c_1 < 1/2; the backtrack_condition at the beginning
                # of the function guarantees at least a backtracking factor ρ.
                α_tmp = - (dϕ_0 * α_2*α_2) / ( T(2) * (ϕx_1 - ϕ_0 - dϕ_0*α_2) )
            else
                div = one(T) / (α_1*α_1 * α_2*α_2 * (α_2 - α_1))
                a = (α_1*α_1*(ϕx_1 - ϕ_0 - dϕ_0*α_2) - α_2*α_2*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div
                b = (-α_1*α_1*α_1*(ϕx_1 - ϕ_0 - dϕ_0*α_2) + α_2*α_2*α_2*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div

                if norm(a) <= eps(Float64) + sqrttol*norm(a)
                    α_tmp = dϕ_0 / (T(2)*b)
                else
                    # discriminant
                    d = max(b*b - T(3)*a*dϕ_0, zero(T))
                    # quadratic equation root
                    α_tmp = (sqrt(d) - b) / (T(3)*a)
                end
            end
            α_1 = α_2

            α_tmp = NaNMath.min(α_tmp, α_2*ρ_hi) # avoid too small reductions
            α_2 = NaNMath.max(α_tmp, α_2*ρ_lo) # avoid too big reductions
            # α_tmp = min(α_tmp, α_2*ρ_hi) # avoid too small reductions
            # α_2 = max(α_tmp, α_2*ρ_lo) # avoid too big reductions

            # Evaluate f(x) at proposed position
            # ϕx_0, ϕx_1 = ϕx_1, f(x + α_2*s); f_calls += 1;
            @fastmath @inbounds @simd for i ∈ 1:L
                x_new[i] = x_old[i] + α_2*s[i]
            end
            # SIMDArrays.vadd!(x_new, x_old, α_2, s)
            ϕx_0, ϕx_1 = ϕx_1, f(obj, x_new); f_calls += 1;
        end
        alpha, fpropose = α_2, ϕx_1

        # s = alpha*s
        # x = x + s # Update x
        # @inbounds @simd for i ∈ 1:L
        #     s[i] *= alpha
        #     x_old[i] += s[i]
        # end
        update_state!(s, x_old, alpha)
        # jold = copy(jx)
        copyto!(∇_old, ∇)
    end
    # return StaticOptimizationResults(NaN, N, tol, f_calls, g_calls, false), x_old
    NaN
end

@generated function update_state!(C::AbstractFixedSizeArray{S,T,N,R,L}, B::AbstractFixedSizeArray{S,T,N,R,L}, α::T) where {S,T,N,R,L}
    T_size = sizeof(T)
    VL = min(VectorizationBase.REGISTER_SIZE ÷ T_size, L)
    VLT = VL * T_size
    V = Vec{VL,T}

    iter = L ÷ VL
    q = quote
        ptr_C = pointer(C)
        ptr_B = pointer(B)
        vα = vbroadcast($V, α)
    end

    rep, rem = divrem(iter, 4)
    if (rep == 1 && rem == 0) || (rep >= 1 && rem != 0)
        rep -= 1
        rem += 4
    end
    if rep > 0
        push!(q.args,
            quote
                for i ∈ 0:$(4VLT):$(4VLT*(rep-1))
                    vs_0 = vmul(vload($V, ptr_C + i), vα)
                    vstore!(ptr_C + i, vs_0)
                    vstore!(ptr_B + i, vadd(vload($V, ptr_B + i), vs_0))

                    vs_1 = vmul(vload($V, ptr_C + i + $VLT), vα)
                    vstore!(ptr_C + i + $VLT, vs_1)
                    vstore!(ptr_B + i + $VLT, vadd(vload($V, ptr_B + i + $VLT), vs_1))

                    vs_2 = vmul(vload($V, ptr_C + i + $(2VLT)), vα)
                    vstore!(ptr_C + i + $(2VLT), vs_2)
                    vstore!(ptr_B + i + $(2VLT), vadd(vload($V, ptr_B + i + $(2VLT)), vs_2))

                    vs_3 = vmul(vload($V, ptr_C + i + $(3VLT)), vα)
                    vstore!(ptr_C + i + $(3VLT), vs_3)
                    vstore!(ptr_B + i + $(3VLT), vadd(vload($V, ptr_B + i + $(3VLT)), vs_3))
                end
            end
        )
    end
    for i ∈ 0:(rem-1)
        offset = VLT*(i + 4rep)
        push!(q.args,
            quote
                $(Symbol(:vs_,i)) = vmul(vload($V, ptr_C + $offset), vα)
                vstore!(ptr_C + $offset, $(Symbol(:vs_,i)))
                vstore!(ptr_B + $offset, vadd(vload($V, ptr_B + $offset), $(Symbol(:vs_,i))))
            end
        )
    end

    push!(q.args, :(nothing))
    q
end

"""
Similar to optimize_light!, but it scales the function so that
norm(gradient(f, initial_x)) ≈ 10
"""
function optimize_scale!(state, obj, x::AbstractFixedSizeVector{P,T,L}, ls::BackTracking2{order}, scale_target=T(10), tol = T(1e-8)) where {P,T,L,order}
    # res = DiffResults.GradientResult(x)
    # ls = BackTracking()
    # order = ordernum(bto)
    # copyto!(state.xinit, x)
    # copyto!(state.x_new, x)
    # xinit = copy(x)
    # x_new = copy(x)
    x_old = ref_x_old(state)
    ∇_old = ref_∇_old(state)
    invH = ref_invH(state)
    δ∇ = ref_δ∇(state)
    s = ref_s(state)
    u = ref_u(state)
    initial_invH!(state)
    x_new = ref_x_new(state)
    copyto!(x_old, x)
    # hx = SMatrix{P,P,T}(I)
    # if !(hguess isa Nothing)
    #     hx = hguess * hx
    # end
    # hold = copy(hx)
    # jold = copy(x); s = copy(x)
    # @unpack c_1, ρ_hi, ρ_lo, iterations = ls
    c_1, ρ_hi, ρ_lo, iterations = ls.c_1, ls.ρ_hi, ls.ρ_lo, ls.iterations
    iterfinitemax = round(Int,-log2(eps(eltype(x))))
    sqrttol = sqrt(eps(T))
    α_0 = one(T)
    N = 200
#    f_calls = 0
#    g_calls = 0
    local scale::T
    local ϕ_0::T
    for n = 1:N
        # res = ForwardDiff.gradient!(res, f, x); f_calls +=1; g_calls +=1; # Obtain gradient
        # ∇ = gradient!(d, x_old)
        # # ϕ_0 = DiffResults.value(res)
        # ϕ_0 = value(d)

        if n == 1 # calculate scale
            ϕ_0, scale = scale_fdf(obj, x_old, scale_target);# f_calls +=1; g_calls +=1;
            # @show scale
            ∇ = gradient(obj)

            isfinite(ϕ_0) || return T(NaN), scale
            maximum(abs, ∇) < tol && return ϕ_0, scale
        else # update hessian
            ϕ_old = ϕ_0
            ϕ_0 = scaled_fdf(obj, x_old, scale);# f_calls +=1; g_calls +=1;
            ∇ = gradient(obj)

            norm(ϕ_0 - ϕ_old) <= tol*max(norm(ϕ_0),norm(ϕ_old)) && return ϕ_0, scale
            isfinite(ϕ_0) || return T(NaN), scale
            maximum(abs, ∇) < tol && return ϕ_0, scale
            # y = jx - jold
            # hx = norm(y) < eps(eltype(x)) ? hx : hx + y*y' / (y'*s) - (hx*(s*s')*hx)/(s'*hx*s)

            # invH = invH + c1 * (s * s') - c2 * (u * s' + s * u')
            # @show L
            # @show length(δ∇), length(∇), length(∇_old)
            # @show full_length(δ∇), full_length(∇), full_length(∇_old)
            # @show typeof(δ∇), typeof(∇), typeof(∇_old)
            @inbounds @simd for i ∈ 1:L
                δ∇[i] =  ∇[i] - ∇_old[i]
            end
            # SIMDArrays.vsub!(δ∇, ∇, ∇_old)
            # Update the inverse Hessian approximation using Sherman-Morrison
            dx_dg = real(dot(s, δ∇))
            # dx_dg == 0.0 && return true # force stop
            mul!(u, invH, δ∇)
            c2 = one(T) / dx_dg
            c1 = fma(real(dot(δ∇, u)), c2*c2, c2)
            BFGS_update!(invH, s, u, c1, c2)
        end
        # @show size(s), size(invH), size(∇)
        # @show typeof(s), typeof(invH), typeof(∇)
        mul!(s, invH, ∇)
        # SIMDArrays.scale!(s, -one(eltype(s)))
        @inbounds @simd for i ∈ 1:L
            s[i] *= -1
        end
        # SIMDArrays.reflect!(s)
        # s = -hx\jx # Obtain direction
        dϕ_0 = dot(∇, s)

        if dϕ_0 >= zero(T) # If bad, reset search direction
            initial_invH!(state)
            # hx = hold
            # s = -jx
            # @show L, s, ∇
            # @show length(s)
            # @show length(∇)
            @inbounds @simd for i ∈ 1:L
                s[i] = -∇[i]
            end
            # SIMDArrays.reflect!(s, ∇)
            dϕ_0 = dot(∇, s)
        end
        #### Perform line search

        # Count the total number of iterations
        iteration = 0
        ϕx_0, ϕx_1 = ϕ_0, ϕ_0
        α_1, α_2 = α_0, α_0
        @inbounds @simd for i ∈ 1:L
            x_new[i] = x_old[i] + α_1*s[i]
        end
        # SIMDArrays.vadd!(x_new, x_old, α_1, s)
        # ϕx_1 = f(x + α_1*s); f_calls += 1;
        ϕx_1_before = ϕx_1
        ϕx_1 = f(obj, x_new) * scale;# f_calls += 1;
        # @show typeof(ϕx_1_before), typeof(ϕx_1), typeof(f(obj, state.x_new)), typeof(scale)

        # Hard-coded backtrack until we find a finite function value
        iterfinite = 0
        while !isfinite(ϕx_1) && iterfinite < iterfinitemax
            iterfinite += 1
            α_1 = α_2
            α_2 = T(0.5)*α_1
            @inbounds @simd for i ∈ 1:L
                x_new[i] = x_old[i] + α_2*s[i]
            end
            # SIMDArrays.vadd!(x_new, x_old, α_2, s)
            # ϕx_1 = f(x + α_2*s); f_calls += 1;
            ϕx_1 = f(obj, x_new) * scale;# f_calls += 1;
        end

        # Backtrack until we satisfy sufficient decrease condition
        while ϕx_1 > ϕ_0 + c_1 * α_2 * dϕ_0
            # Increment the number of steps we've had to perform
            iteration += 1

            # Ensure termination
            iteration > iterations && linesearch_failure(iterations)

            # Shrink proposed step-size:
            @fastmath if order == 2 || iteration == 1
                # backtracking via quadratic interpolation:
                # This interpolates the available data
                #    f(0), f'(0), f(α)
                # with a quadractic which is then minimised; this comes with a
                # guaranteed backtracking factor 0.5 * (1-c_1)^{-1} which is < 1
                # provided that c_1 < 1/2; the backtrack_condition at the beginning
                # of the function guarantees at least a backtracking factor ρ.
                # @show typeof(dϕ_0), typeof(ϕx_1), typeof(ϕ_0)
                α_tmp = - (dϕ_0 * α_2*α_2) / ( T(2) * (ϕx_1 - ϕ_0 - dϕ_0*α_2) )
            else
                div = one(T) / (α_1*α_1 * α_2*α_2 * (α_2 - α_1))
                a = (α_1*α_1*(ϕx_1 - ϕ_0 - dϕ_0*α_2) - α_2*α_2*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div
                b = (-α_1*α_1*α_1*(ϕx_1 - ϕ_0 - dϕ_0*α_2) + α_2*α_2*α_2*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div
                # @show norm(a) <= eps(T) + sqrttol*norm(a)
                if norm(a) <= eps(T) + sqrttol*norm(a)
                    α_tmp = dϕ_0 / (T(2)*b)
                else
                    # discriminant
                    d = max(b*b - T(3)*a*dϕ_0, zero(T))
                    # quadratic equation root
                    α_tmp = (sqrt(d) - b) / (T(3)*a)
                end
            end
            α_1 = α_2

            # @show typeof(α_tmp)
            α_tmp = NaNMath.min(α_tmp, α_2*ρ_hi) # avoid too small reductions
            α_2 = NaNMath.max(α_tmp, α_2*ρ_lo) # avoid too big reductions
            # α_tmp = min(α_tmp, α_2*ρ_hi) # avoid too small reductions
            # α_2 = max(α_tmp, α_2*ρ_lo) # avoid too big reductions

            # Evaluate f(x) at proposed position
            # ϕx_0, ϕx_1 = ϕx_1, f(x + α_2*s); f_calls += 1;
            @fastmath @inbounds @simd for i ∈ 1:L
                x_new[i] = x_old[i] + α_2*s[i]
            end
            # SIMDArrays.vadd!(x_new, x_old, α_2::T, s)
            ϕx_0, ϕx_1 = ϕx_1, f(obj, x_new) * scale; #f_calls += 1;
        end
        alpha, fpropose = α_2, ϕx_1

        # s = alpha*s
        # x = x + s # Update x
        # @inbounds @simd for i ∈ 1:L
        #     s[i] *= alpha
        #     x_old[i] += s[i]
        # end
        update_state!(s, x_old, alpha)
        # jold = copy(jx)
        copyto!(∇_old, ∇)
    end
    # return StaticOptimizationResults(NaN, N, tol, f_calls, g_calls, false), x_old
    T(NaN), scale
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Provide the [`SharedArray`](@ref) type. It represents an array, which is shared across multiple processes, on a single machine.
"""
module SharedArrays

using Mmap, Distributed, Random

import Base: length, size, ndims, IndexStyle, reshape, convert, deepcopy_internal,
             show, getindex, setindex!, fill!, similar, reduce, map!, copyto!, unsafe_convert
import Random
using Serialization
using Serialization: serialize_cycle_header, serialize_type, writetag, UNDEFREF_TAG, serialize, deserialize
import Serialization: serialize, deserialize
import Distributed: RRID, procs, remotecall_fetch
import Base.Filesystem: JL_O_CREAT, JL_O_RDWR, S_IRUSR, S_IWUSR

export SharedArray, SharedVector, SharedMatrix, sdata, indexpids, localindices

mutable struct SharedArray{T,N} <: DenseArray{T,N}
    id::RRID
    dims::NTuple{N,Int}
    pids::Vector{Int}
    refs::Vector

    # The segname is currently used only in the test scripts to ensure that
    # the shmem segment has been unlinked.
    segname::String

    # Fields below are not to be serialized
    # Local shmem map.
    s::Array{T,N}

    # idx of current worker's pid in the pids vector, 0 if this shared array is not mapped locally.
    pidx::Int

    # the local partition into the array when viewed as a single dimensional array.
    # this can be removed when @distributed or its equivalent supports looping on
    # a subset of workers.
    loc_subarr_1d::SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}

    function SharedArray{T,N}(d,p,r,sn,s) where {T,N}
        S = new(RRID(),d,p,r,sn,s,0,view(Array{T}(undef, ntuple(d->0,N)), 1:0))
        sa_refs[S.id] = WeakRef(S)
        S
    end
end

const sa_refs = Dict{RRID, WeakRef}()

"""
    SharedArray{T}(dims::NTuple; init=false, pids=Int[])
    SharedArray{T,N}(...)

Construct a `SharedArray` of a bits type `T` and size `dims` across the
processes specified by `pids` - all of which have to be on the same
host.  If `N` is specified by calling `SharedArray{T,N}(dims)`, then
`N` must match the length of `dims`.

If `pids` is left unspecified, the shared array will be mapped across all processes on the
current host, including the master. But, `localindices` and `indexpids` will only refer to
worker processes. This facilitates work distribution code to use workers for actual
computation with the master process acting as a driver.

If an `init` function of the type `initfn(S::SharedArray)` is specified, it is called on all
the participating workers.

The shared array is valid as long as a reference to the `SharedArray` object exists on the node
which created the mapping.

    SharedArray{T}(filename::AbstractString, dims::NTuple, [offset=0]; mode=nothing, init=false, pids=Int[])
    SharedArray{T,N}(...)

Construct a `SharedArray` backed by the file `filename`, with element
type `T` (must be a bits type) and size `dims`, across the processes
specified by `pids` - all of which have to be on the same host. This
file is mmapped into the host memory, with the following consequences:

- The array data must be represented in binary format (e.g., an ASCII
  format like CSV cannot be supported)

- Any changes you make to the array values (e.g., `A[3] = 0`) will
  also change the values on disk

If `pids` is left unspecified, the shared array will be mapped across
all processes on the current host, including the master. But,
`localindices` and `indexpids` will only refer to worker
processes. This facilitates work distribution code to use workers for
actual computation with the master process acting as a driver.

`mode` must be one of `"r"`, `"r+"`, `"w+"`, or `"a+"`, and defaults
to `"r+"` if the file specified by `filename` already exists, or
`"w+"` if not. If an `init` function of the type
`initfn(S::SharedArray)` is specified, it is called on all the
participating workers. You cannot specify an `init` function if the
file is not writable.

`offset` allows you to skip the specified number of bytes at the
beginning of the file.
"""
SharedArray

function SharedArray{T,N}(dims::Dims{N}; init=false, pids=Int[]) where {T,N}
    isbitstype(T) || throw(ArgumentError("type of SharedArray elements must be bits types, got $(T)"))

    pids, onlocalhost = shared_pids(pids)

    local shm_seg_name = ""
    local s = Array{T}(undef, ntuple(d->0,N))
    local S
    local shmmem_create_pid
    try
        # On OSX, the shm_seg_name length must be <= 31 characters (including the terminating NULL character)
        shm_seg_name = "/jl$(lpad(string(getpid() % 10^6), 6, "0"))$(randstring(20))"
        if onlocalhost
            shmmem_create_pid = myid()
            s = shm_mmap_array(T, dims, shm_seg_name, JL_O_CREAT | JL_O_RDWR)
        else
            # The shared array is created on a remote machine
            shmmem_create_pid = pids[1]
            remotecall_fetch(pids[1]) do
                shm_mmap_array(T, dims, shm_seg_name, JL_O_CREAT | JL_O_RDWR)
                nothing
            end
        end

        func_mapshmem = () -> shm_mmap_array(T, dims, shm_seg_name, JL_O_RDWR)

        refs = Vector{Future}(undef, length(pids))
        for (i, p) in enumerate(pids)
            refs[i] = remotecall(func_mapshmem, p)
        end

        # Wait till all the workers have mapped the segment
        for ref in refs
            wait(ref)
        end

        # All good, immediately unlink the segment.
        if (prod(dims) > 0) && (sizeof(T) > 0)
            if onlocalhost
                rc = shm_unlink(shm_seg_name)
            else
                rc = remotecall_fetch(shm_unlink, shmmem_create_pid, shm_seg_name)
            end
            systemerror("Error unlinking shmem segment " * shm_seg_name, rc != 0)
        end
        S = SharedArray{T,N}(dims, pids, refs, shm_seg_name, s)
        initialize_shared_array(S, onlocalhost, init, pids)
        shm_seg_name = ""

    finally
        if !isempty(shm_seg_name)
            remotecall_fetch(shm_unlink, shmmem_create_pid, shm_seg_name)
        end
    end
    S
end

SharedArray{T,N}(I::Integer...; kwargs...) where {T,N} =
    SharedArray{T,N}(I; kwargs...)
SharedArray{T}(d::NTuple; kwargs...) where {T} =
    SharedArray{T,length(d)}(d; kwargs...)
SharedArray{T}(I::Integer...; kwargs...) where {T} =
    SharedArray{T,length(I)}(I; kwargs...)
SharedArray{T}(m::Integer; kwargs...) where {T} =
    SharedArray{T,1}(m; kwargs...)
SharedArray{T}(m::Integer, n::Integer; kwargs...) where {T} =
    SharedArray{T,2}(m, n; kwargs...)
SharedArray{T}(m::Integer, n::Integer, o::Integer; kwargs...) where {T} =
    SharedArray{T,3}(m, n, o; kwargs...)

function SharedArray{T,N}(filename::AbstractString, dims::NTuple{N,Int}, offset::Integer=0;
                          mode=nothing, init=false, pids::Vector{Int}=Int[]) where {T,N}
    if !isabspath(filename)
        throw(ArgumentError("$filename is not an absolute path; try abspath(filename)?"))
    end
    if !isbitstype(T)
        throw(ArgumentError("type of SharedArray elements must be bits types, got $(T)"))
    end

    pids, onlocalhost = shared_pids(pids)

    # If not supplied, determine the appropriate mode
    have_file = onlocalhost ? isfile(filename) : remotecall_fetch(isfile, pids[1], filename)
    if mode === nothing
        mode = have_file ? "r+" : "w+"
    end
    workermode = mode == "w+" ? "r+" : mode  # workers don't truncate!

    # Ensure the file will be readable
    if !(mode in ("r", "r+", "w+", "a+"))
        throw(ArgumentError("mode must be readable, but $mode is not"))
    end
    if init !== false
        typeassert(init, Function)
        if !(mode in ("r+", "w+", "a+"))
            throw(ArgumentError("cannot initialize unwritable array (mode = $mode)"))
        end
    end
    if mode == "r" && !isfile(filename)
        throw(ArgumentError("file $filename does not exist, but mode $mode cannot create it"))
    end

    # Create the file if it doesn't exist, map it if it does
    refs = Vector{Future}(undef, length(pids))
    func_mmap = mode -> open(filename, mode) do io
        Mmap.mmap(io, Array{T,N}, dims, offset; shared=true)
    end
    s = Array{T}(undef, ntuple(d->0,N))
    if onlocalhost
        s = func_mmap(mode)
        refs[1] = remotecall(pids[1]) do
            func_mmap(workermode)
        end
    else
        refs[1] = remotecall_wait(pids[1]) do
            func_mmap(mode)
        end
    end

    # Populate the rest of the workers
    for i = 2:length(pids)
        refs[i] = remotecall(pids[i]) do
            func_mmap(workermode)
        end
    end

    # Wait till all the workers have mapped the segment
    for ref in refs
        wait(ref)
    end

    S = SharedArray{T,N}(dims, pids, refs, filename, s)
    initialize_shared_array(S, onlocalhost, init, pids)
    S
end

SharedArray{T}(filename::AbstractString, dims::NTuple{N,Int}, offset::Integer=0;
               mode=nothing, init=false, pids::Vector{Int}=Int[]) where {T,N} =
    SharedArray{T,N}(filename, dims, offset; mode=mode, init=init, pids=pids)

function initialize_shared_array(S, onlocalhost, init, pids)
    if onlocalhost
        init_loc_flds(S)
    else
        S.pidx = 0
    end

    # if present, init function is called on each of the parts
    if isa(init, Function)
        @sync begin
            for p in pids
                @async remotecall_wait(init, p, S)
            end
        end
    end

    finalizer(finalize_refs, S)
    S
end

function finalize_refs(S::SharedArray{T,N}) where T where N
    if length(S.pids) > 0
        for r in S.refs
            finalize(r)
        end
        empty!(S.pids)
        empty!(S.refs)
        init_loc_flds(S)
        S.s = Array{T}(undef, ntuple(d->0,N))
        delete!(sa_refs, S.id)
    end
    S
end

const SharedVector{T} = SharedArray{T,1}
const SharedMatrix{T} = SharedArray{T,2}

SharedVector(A::Vector) = SharedArray(A)
SharedMatrix(A::Matrix) = SharedArray(A)

length(S::SharedArray) = prod(S.dims)
size(S::SharedArray) = S.dims
ndims(S::SharedArray) = length(S.dims)
IndexStyle(::Type{<:SharedArray}) = IndexLinear()

function reshape(a::SharedArray{T}, dims::NTuple{N,Int}) where {T,N}
    if length(a) != prod(dims)
        throw(DimensionMismatch("dimensions must be consistent with array size"))
    end
    refs = Vector{Future}(undef, length(a.pids))
    for (i, p) in enumerate(a.pids)
        refs[i] = remotecall(p, a.refs[i], dims) do r,d
            reshape(fetch(r),d)
        end
    end

    A = SharedArray{T,N}(dims, a.pids, refs, a.segname, reshape(a.s, dims))
    init_loc_flds(A)
    A
end

"""
    procs(S::SharedArray)

Get the vector of processes mapping the shared array.
"""
procs(S::SharedArray) = S.pids

"""
    indexpids(S::SharedArray)

Returns the current worker's index in the list of workers
mapping the `SharedArray` (i.e. in the same list returned by `procs(S)`), or
0 if the `SharedArray` is not mapped locally.
"""
indexpids(S::SharedArray) = S.pidx

"""
    sdata(S::SharedArray)

Returns the actual `Array` object backing `S`.
"""
sdata(S::SharedArray) = S.s
sdata(A::AbstractArray) = A

"""
    localindices(S::SharedArray)

Returns a range describing the "default" indices to be handled by the
current process.  This range should be interpreted in the sense of
linear indexing, i.e., as a sub-range of `1:length(S)`.  In
multi-process contexts, returns an empty range in the parent process
(or any process for which [`indexpids`](@ref) returns 0).

It's worth emphasizing that `localindices` exists purely as a
convenience, and you can partition work on the array among workers any
way you wish. For a `SharedArray`, all indices should be equally fast
for each worker process.
"""
localindices(S::SharedArray) = S.pidx > 0 ? range_1dim(S, S.pidx) : 1:0

unsafe_convert(::Type{Ptr{T}}, S::SharedArray{T}) where {T} = unsafe_convert(Ptr{T}, sdata(S))
unsafe_convert(::Type{Ptr{T}}, S::SharedArray   ) where {T} = unsafe_convert(Ptr{T}, sdata(S))

function SharedArray(A::Array)
    S = SharedArray{eltype(A),ndims(A)}(size(A))
    copyto!(S, A)
end
function SharedArray{T}(A::Array) where T
    S = SharedArray{T,ndims(A)}(size(A))
    copyto!(S, A)
end
function SharedArray{TS,N}(A::Array{TA,N}) where {TS,TA,N}
    S = SharedArray{TS,ndims(A)}(size(A))
    copyto!(S, A)
end

convert(T::Type{<:SharedArray}, a::Array) = T(a)

function deepcopy_internal(S::SharedArray, stackdict::IdDict)
    haskey(stackdict, S) && return stackdict[S]
    R = SharedArray{eltype(S),ndims(S)}(size(S); pids = S.pids)
    copyto!(sdata(R), sdata(S))
    stackdict[S] = R
    return R
end

function shared_pids(pids)
    if isempty(pids)
        # only use workers on the current host
        pids = procs(myid())
        if length(pids) > 1
            pids = filter(x -> x != 1, pids)
        end

        onlocalhost = true
    else
        if !check_same_host(pids)
            throw(ArgumentError("SharedArray requires all requested processes to be on the same machine."))
        end

        onlocalhost = myid() in procs(pids[1])
    end
    pids, onlocalhost
end

function range_1dim(S::SharedArray, pidx)
    l = length(S)
    nw = length(S.pids)
    partlen = div(l, nw)

    if l < nw
        if pidx <= l
            return pidx:pidx
        else
            return 1:0
        end
    elseif pidx == nw
        return (((pidx-1) * partlen) + 1):l
    else
        return (((pidx-1) * partlen) + 1):(pidx*partlen)
    end
end

sub_1dim(S::SharedArray, pidx) = view(S.s, range_1dim(S, pidx))

function init_loc_flds(S::SharedArray{T,N}, empty_local=false) where T where N
    if myid() in S.pids
        S.pidx = findfirst(isequal(myid()), S.pids)
        if isa(S.refs[1], Future)
            refid = remoteref_id(S.refs[S.pidx])
        else
            refid = S.refs[S.pidx]
        end
        c = channel_from_id(refid)
        S.s = fetch(c)
        S.loc_subarr_1d = sub_1dim(S, S.pidx)
    else
        S.pidx = 0
        if empty_local
            S.s = Array{T}(undef, ntuple(d->0,N))
        end
        S.loc_subarr_1d = view(Array{T}(undef, ntuple(d->0,N)), 1:0)
    end
end

# Don't serialize s (it is the complete array) and
# pidx, which is relevant to the current process only
function serialize(s::AbstractSerializer, S::SharedArray)
    serialize_cycle_header(s, S) && return

    destpid = worker_id_from_socket(s.io)
    if S.id.whence == destpid
        # The shared array was created from destpid, hence a reference to it
        # must be available at destpid.
        serialize(s, true)
        serialize(s, S.id.whence)
        serialize(s, S.id.id)
        return
    end
    serialize(s, false)
    for n in fieldnames(SharedArray)
        if n in [:s, :pidx, :loc_subarr_1d]
            writetag(s.io, UNDEFREF_TAG)
        elseif n == :refs
            v = getfield(S, n)
            if isa(v[1], Future)
                # convert to ids to avoid distributed GC overhead
                ids = [remoteref_id(x) for x in v]
                serialize(s, ids)
            else
                serialize(s, v)
            end
        else
            serialize(s, getfield(S, n))
        end
    end
end

function deserialize(s::AbstractSerializer, t::Type{<:SharedArray})
    ref_exists = deserialize(s)
    if ref_exists
        sref = sa_refs[RRID(deserialize(s), deserialize(s))]
        if sref.value !== nothing
            return sref.value
        end
        error("Expected reference to shared array instance not found")
    end

    S = invoke(deserialize, Tuple{AbstractSerializer,DataType}, s, t)
    init_loc_flds(S, true)
    return S
end

function show(io::IO, S::SharedArray)
    if length(S.s) > 0
        invoke(show, Tuple{IO,DenseArray}, io, S)
    else
        show(io, remotecall_fetch(sharr->sharr.s, S.pids[1], S))
    end
end

function show(io::IO, mime::MIME"text/plain", S::SharedArray)
    if length(S.s) > 0
        invoke(show, Tuple{IO,MIME"text/plain",DenseArray}, io, MIME"text/plain"(), S)
    else
        # retrieve from the first worker mapping the array.
        summary(io, S); println(io, ":")
        Base.print_array(io, remotecall_fetch(sharr->sharr.s, S.pids[1], S))
    end
end

Array(S::SharedArray) = S.s

# pass through getindex and setindex! - unlike DArrays, these always work on the complete array
getindex(S::SharedArray, i::Real) = getindex(S.s, i)

setindex!(S::SharedArray, x, i::Real) = setindex!(S.s, x, i)

function fill!(S::SharedArray, v)
    vT = convert(eltype(S), v)
    f = S->fill!(S.loc_subarr_1d, vT)
    @sync for p in procs(S)
        @async remotecall_wait(f, p, S)
    end
    return S
end

function Random.rand!(S::SharedArray{T}) where T
    f = S->map!(x -> rand(T), S.loc_subarr_1d, S.loc_subarr_1d)
    @sync for p in procs(S)
        @async remotecall_wait(f, p, S)
    end
    return S
end

function Random.randn!(S::SharedArray)
    f = S->map!(x -> randn(), S.loc_subarr_1d, S.loc_subarr_1d)
    @sync for p in procs(S)
        @async remotecall_wait(f, p, S)
    end
    return S
end

# convenience constructors
function shmem_fill(v, dims; kwargs...)
    SharedArray{typeof(v),length(dims)}(dims; init = S->fill!(S.loc_subarr_1d, v), kwargs...)
end
shmem_fill(v, I::Int...; kwargs...) = shmem_fill(v, I; kwargs...)

# rand variant with range
function shmem_rand(TR::Union{DataType, UnitRange}, dims; kwargs...)
    if isa(TR, UnitRange)
        SharedArray{Int,length(dims)}(dims; init = S -> map!(x -> rand(TR), S.loc_subarr_1d, S.loc_subarr_1d), kwargs...)
    else
        SharedArray{TR,length(dims)}(dims; init = S -> map!(x -> rand(TR), S.loc_subarr_1d, S.loc_subarr_1d), kwargs...)
    end
end
shmem_rand(TR::Union{DataType, UnitRange}, i::Int; kwargs...) = shmem_rand(TR, (i,); kwargs...)
shmem_rand(TR::Union{DataType, UnitRange}, I::Int...; kwargs...) = shmem_rand(TR, I; kwargs...)

shmem_rand(dims; kwargs...) = shmem_rand(Float64, dims; kwargs...)
shmem_rand(I::Int...; kwargs...) = shmem_rand(I; kwargs...)

function shmem_randn(dims; kwargs...)
    SharedArray{Float64,length(dims)}(dims; init = S-> map!(x -> randn(), S.loc_subarr_1d, S.loc_subarr_1d), kwargs...)
end
shmem_randn(I::Int...; kwargs...) = shmem_randn(I; kwargs...)

similar(S::SharedArray, T::Type, dims::Dims) = similar(S.s, T, dims)
similar(S::SharedArray, T::Type) = similar(S.s, T, size(S))
similar(S::SharedArray, dims::Dims) = similar(S.s, eltype(S), dims)
similar(S::SharedArray) = similar(S.s, eltype(S), size(S))

reduce(f, S::SharedArray) =
    mapreduce(fetch, f, Any[ @spawnat p reduce(f, S.loc_subarr_1d) for p in procs(S) ])

reduce(::typeof(vcat), S::SharedVector) = invoke(reduce, Tuple{Any,SharedArray}, vcat, S)
reduce(::typeof(hcat), S::SharedVector) = invoke(reduce, Tuple{Any,SharedArray}, hcat, S)

function map!(f, S::SharedArray, Q::SharedArray)
    if (S !== Q) && (procs(S) != procs(Q) || localindices(S) != localindices(Q))
        throw(ArgumentError("incompatible source and destination arguments"))
    end
    @sync for p in procs(S)
        @spawnat p begin
            for idx in localindices(S)
                S.s[idx] = f(Q.s[idx])
            end
        end
    end
    return S
end

copyto!(S::SharedArray, A::Array) = (copyto!(S.s, A); S)

function copyto!(S::SharedArray, R::SharedArray)
    length(S) == length(R) || throw(BoundsError())
    ps = intersect(procs(S), procs(R))
    isempty(ps) && throw(ArgumentError("source and destination arrays don't share any process"))
    l = length(S)
    length(ps) > l && (ps = ps[1:l])
    nw = length(ps)
    partlen = div(l, nw)

    @sync for i = 1:nw
        p = ps[i]
        idx = i < nw ?  ((i-1)*partlen+1:i*partlen) : ((i-1)*partlen+1:l)
        @spawnat p begin
            S.s[idx] = R.s[idx]
        end
    end

    return S
end

function print_shmem_limits(slen)
    try
        if Sys.islinux()
            pfx = "kernel"
        elseif Sys.isapple()
            pfx = "kern.sysv"
        elseif Sys.KERNEL == :FreeBSD || Sys.KERNEL == :DragonFly
            pfx = "kern.ipc"
        elseif Sys.KERNEL == :OpenBSD
            pfx = "kern.shminfo"
        else
            # seems NetBSD does not have *.shmall
            return
        end

        shmmax_MB = div(parse(Int, split(read(`sysctl $(pfx).shmmax`, String))[end]), 1024*1024)
        page_size = parse(Int, split(read(`getconf PAGE_SIZE`, String))[end])
        shmall_MB = div(parse(Int, split(read(`sysctl $(pfx).shmall`, String))[end]) * page_size, 1024*1024)

        println("System max size of single shmem segment(MB) : ", shmmax_MB,
            "\nSystem max size of all shmem segments(MB) : ", shmall_MB,
            "\nRequested size(MB) : ", div(slen, 1024*1024),
            "\nPlease ensure requested size is within system limits.",
            "\nIf not, increase system limits and try again."
        )
    catch e
        nothing # Ignore any errors in this
    end
end

# utilities
function shm_mmap_array(T, dims, shm_seg_name, mode)
    local s = nothing
    local A = nothing

    if (prod(dims) == 0) || (sizeof(T) == 0)
        return Array{T}(undef, dims)
    end

    try
        A = _shm_mmap_array(T, dims, shm_seg_name, mode)
    catch
        print_shmem_limits(prod(dims)*sizeof(T))
        rethrow()

    finally
        if s !== nothing
            close(s)
        end
    end
    A
end

# platform-specific code

if Sys.iswindows()
function _shm_mmap_array(T, dims, shm_seg_name, mode)
    readonly = !((mode & JL_O_RDWR) == JL_O_RDWR)
    create = (mode & JL_O_CREAT) == JL_O_CREAT
    s = Mmap.Anonymous(shm_seg_name, readonly, create)
    Mmap.mmap(s, Array{T,length(dims)}, dims, zero(Int64))
end

# no-op in windows
shm_unlink(shm_seg_name) = 0

else # !windows
function _shm_mmap_array(T, dims, shm_seg_name, mode)
    fd_mem = shm_open(shm_seg_name, mode, S_IRUSR | S_IWUSR)
    systemerror("shm_open() failed for " * shm_seg_name, fd_mem < 0)

    s = fdio(fd_mem, true)

    # On OSX, ftruncate must to used to set size of segment, just lseek does not work.
    # and only at creation time
    if (mode & JL_O_CREAT) == JL_O_CREAT
        rc = ccall(:jl_ftruncate, Cint, (Cint, Int64), fd_mem, prod(dims)*sizeof(T))
        systemerror("ftruncate() failed for shm segment " * shm_seg_name, rc != 0)
    end

    Mmap.mmap(s, Array{T,length(dims)}, dims, zero(Int64); grow=false)
end

shm_unlink(shm_seg_name) = ccall(:shm_unlink, Cint, (Cstring,), shm_seg_name)
shm_open(shm_seg_name, oflags, permissions) = ccall(:shm_open, Cint,
    (Cstring, Cint, Base.Cmode_t), shm_seg_name, oflags, permissions)

end # os-test

end # module

ext: jl

max_stars_repo_name: ashleysommer/julia

max_issues_repo_name: ashleysommer/julia

max_forks_repo_name: ashleysommer/julia


text:
using Documenter, BroadcastableStructs

makedocs(;
    modules=[BroadcastableStructs],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tkf/BroadcastableStructs.jl/blob/{commit}{path}#L{line}",
    sitename="BroadcastableStructs.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
    strict=v"1.2" <= VERSION < v"1.3",
)

deploydocs(;
    repo="github.com/tkf/BroadcastableStructs.jl",
)


text:
"""
    create(; kwargs...)
    create(config::Config, args...; kwargs...)
    create(env::Environment; kwargs...)

Creates a new playground `Environment` including initializing its package directory and installing
any package in the REQUIRE file passed in.

# Optional Arguments
You can optionally pass in an `Environment` instance of a `Config` and args to build one.

# Keywords Arguments
* `julia::AbstractString` - a julia binary to use in this playground environment.
* `reqs_file::AbstractPath` - path to a REQUIRE file of packages to install in this environment.
* `registry::AbstractString` - url to the package registry to be cloned.
* `branch::AbstractString` - registry branch to be checked out.
"""
create(; kwargs...) = create(Environment(); kwargs...)
create(config::Config, args...; kwargs...) = create(Environment(config, args...); kwargs...)

function create(env::Environment; kwargs...)
    debug(logger, "Environment Config: $(env.config)")
    init(env)
    opts = Dict(kwargs)

    julia_exec = if haskey(opts, :julia) && !isempty(opts[:julia])
        join(env.config.bin, opts[:julia])
    else
        join(Path(Base.JULIA_HOME), p"julia")
    end

    debug(logger, "$(julia(env)) -> $julia_exec")
    symlink(julia_exec, julia(env), exist_ok=true, overwrite=true)

    withenv(env) do
        registry = if isempty(get(opts, :registry, ""))
            env.config.default_registry
        else
            opts[:registry]
        end
        branch = if isempty(get(opts, :branch, ""))
            env.config.default_branch
        else
            opts[:branch]
        end
        init_cmd = "Pkg.init(\"$registry\", \"$branch\")"
        Playground.log_output(`$(julia(env)) -e $init_cmd`)

        reqs_file = if haskey(opts, :reqs_file) && !isempty(opts[:reqs_file])
            opts[:reqs_file]
        else
            join(env.config.root, "REQUIRE")
        end

        if exists(reqs_file)
            info(logger, "Installing packages from REQUIRE file $reqs_file...")

            for v in readdir(pkg(env))
                copy(reqs_file, join(pkg(env), v, "REQUIRE"); exist_ok=true, overwrite=true)

                try
                    Playground.log_output(`$(julia(env)) -e 'Pkg.resolve()'`)
                catch
                    warn(logger, string(
                        "Failed to resolve requirements. ",
                        "Perhaps there is something wrong with your REQUIRE file."
                    ))
                end
            end
        end
    end
end


text:
show(import_leadfield(verbose = true))

Xcords, Ycords, Zcords = readSRF(verbose = true)
# scatter(Xcords, -Ycords)


text:
using LinearAlgebra, ForwardDiff, StaticArrays
using IterativeSolvers
using Plots

"""
    particle dynamics
    - 3D particle subject to contact forces

    - configuration: q = (x, y, z) ∈ R³
    - impact force (magnitude): n ∈ R₊
    - friction force: β ∈ R⁴₊
        - contact force: λ = (β, n) ∈ R⁴₊ × R₊
        - friction coefficient: μ ∈ R₊

    Discrete Mechanics and Variational Integrators
        pg. 363
"""
struct Particle
      m # mass
      g # gravity
      μ # friction coefficient

      nq # configuration dimension
end

# mass matrix
function mass_matrix(model)
    m = model.m

    Diagonal(@SVector [m, m, m])
end

# gravity
function gravity(model, q)
    m = model.m
    g = model.g

    @SVector [0.0, 0.0, m * g]
end

# signed distance function
function signed_distance(model, q)
    q[3]
end

# contact force Jacobian
function jacobian(model, q)
    Diagonal(@SVector ones(model.nq))
end

# dynamics
function dynamics(model, q1, q2, q3, λ, h)
    nq = model.nq
    SVector{nq}(mass_matrix(model) * (2.0 * q2 - q1 - q3) / h
        - h * gravity(model, q2)
        + h * jacobian(model, q3)' * λ)
end

"""
    step particle system
        solves 1-step feasibility problem
"""
function _step(q1, q2, h;
    tol = 1.0e-8, z0_scale = 0.001, step_type = :gmres, max_iter = 100)
    # 1-step optimization problem:
    #     find z
    #     s.t. r(z) = 0
    #
    # z = (q, n, b, ψ, η, s)
    #     ψ, η are dual variables for friction subproblem
    #     s are slack variables for convenience

    # initialize
    z = 1.0e-1 * ones(15)
    z[1:3] = copy(q2)

    # if z[3] < 0.1
    #     z[3] = 0.1
    # end

    ρ = 1.0 # barrier parameter
    flag = false

    for k = 1:5
        function r(z)
            # system variables
            q3 = view(z, 1:3)
            n = z[4]
            b = view(z, 5:8)

            λ = [(b[1] - b[3]); (b[2] - b[4]); n] # contact forces
            ϕ = signed_distance(model, q3)        # signed-distance function
            vT = (view(q3, 1:2) - view(q2, 1:2)) ./ h

            # friction subproblem variables
            ψ = z[9]
            η = view(z, 10:13)

            # slack
            s1 = z[14]
            s2 = z[15]

            # action optimality conditions
            [dynamics(model, q1, q2, q3, λ, h);
             s1 - ϕ;
             n * s1 - ρ;

             # maximum dissipation optimality conditions
             [vT; -vT] + ψ * ones(4) - η;
             s2 - (model.μ * n - sum(b));
             ψ * s2 - ρ;
             b .* η .- ρ]
        end

        # Jacobian
        function R(z)
            # differentiate r
            _R = ForwardDiff.jacobian(r, z)

            return _R
        end

        function check_variables(z)
            # q3 = view(z, 1:3)
            n = z[4]
            b = view(z, 5:8)

            # friction subproblem variables
            ψ = z[9]
            η = view(z, 10:13)
            s1 = z[14]
            s2 = z[15]

            if n <= 0.0
                return true
            end

            if any(b .<= 0.0)
                return true
            end

            if any(η .<= 0.0)
                return true
            end

            if ψ <= 0.0
                return true
            end

            if s1 <= 0.0
                return true
            end

            if s2 <= 0.0
                return true
            end

            return false
        end

        extra_iters = 0

        for i = 1:max_iter
            # compute residual, residual Jacobian
            res = r(z)

            if norm(res) < tol
                println("   iter ($i) - norm: $(norm(res))")
                # return z, true
                flag = true
                continue
            end

            jac = R(z)

            # compute step
            # if step_type == :gmres
            #     Δ = gmres(jac, res, abstol = 1.0e-12, maxiter = i + extra_iters)
            # else
            # Δ = (jac' * jac + 0.0e-6 * I) \ (jac' * res) # damped least-squares direction
            Δ = jac \ res
            # end

            # line search the step direction
            α = 1.0

            iter = 0
            while check_variables(z - α * Δ) # backtrack inequalities
                α = 0.5 * α
                # println("   α = $α")
                iter += 1
                if iter > 50
                    @error "backtracking line search fail"
                    flag = false
                    @show n = (z - α * Δ)[4]
                    @show b = view(z - α * Δ, 5:8)
                    @show ψ = (z - α * Δ)[9]
                    @show η = view(z - α * Δ, 10:13)
                    @show s = (z - α * Δ)[14:15]
                    return z, false
                end
            end

            while norm(r(z - α * Δ))^2.0 >= (1.0 - 0.001 * α) * norm(res)^2.0
                α = 0.5 * α
                # println("   α = $α")
                iter += 1
                if iter > 50
                    @error "line search fail"
                    flag = false
                    return z, false
                end
            end

            # update
            z .-= α * Δ
        end

        ρ = 0.1 * ρ
        println("ρ: $ρ")
    end

    return z, flag
end

"""
    simulate
    - solves 1-step feasibility problem for T time steps
    - initial configurations: q1, q2 (note this can encode initial velocity)
    - time step: h
"""
function simulate(q1, q2, T, h;
    z0_scale = 0.0, step_type = :gmres)
    println("simulation")

    # initialize histories
    q = [q1, q2]
    n = [0.0]
    b = [zeros(4)]

    # step
    for t = 1:T
        println("   t = $t")
        z_sol, status = _step(q[end-1], q[end], h,
            z0_scale = z0_scale, step_type = step_type)

        if !status
            @error "failed step (t = $t)"
            return q, n, b
        else
            push!(q, view(z_sol, 1:3))
            push!(n, z_sol[4])
            push!(b, view(z_sol, 5:8))
        end
    end

    return q, n, b
end

# simulation setup
model = Particle(1.0, 9.81, 1.0, 3)
h = 0.01

# initial conditions
# v1 = [1.0; 1.0; 0.0]
# q1 = [0.0; 0.0; 1.0]

v1 = [1.0; 10.0; 0.0]
q1 = [0.0; 0.0; 1.0]

v2 = v1 - gravity(model, q1) * h
q2 = q1 + 0.5 * (v1 + v2) * h

q_sol, y_sol, b_sol = simulate(q1, q2, 1000, h)

plot(hcat(q_sol...)[3:3, :]', xlabel = "", label = "z")
plot!(h .* hcat(y_sol...)', xlabel = "", label = "n", linetype = :steppost)

plot(hcat(q_sol...)[1:2, :]', xlabel = "", label = ["x" "y"])
plot!(h * hcat(b_sol...)', label = ["b1" "b2" "b3" "b4"], linetype = :steppost)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

function bat_findmedian_impl(samples::DensitySampleVector)
    median_params = median(samples)
    (result = median_params,)
end


text:
using JuMP, Ipopt, GLPK
using DRMSMIP
using DualDecomposition
using Random

const DD = DualDecomposition

const rng = Random.MersenneTwister(1234)
"""
a: interest rate
π: unit stock price
ρ: unit dividend price

K = 3 number of stages
L = 2 number of investment vehicles
2^L scenarios in each stage
2^L^(K-1)=16 scenarios in total
ρ = 0.05 * π
bank: interest rate 0.01
asset1: 1.03 or 0.97
asset2: 1.06 or 0.94

In each node, we have Np=10 samples from a log-normal distribution
"""

function create_tree(K::Int, L::Int, Np::Int)::DRMSMIP.DR_Tree
    π = ones(L)
    π_samp = generate_sample(L, π, Np)
    set = DRMSMIP.WassersteinSet(π_samp, 1.0, DRMSMIP.norm_L1)
    cost = zeros(L+1)
    tree = DRMSMIP.DR_Tree(π, set, cost)
    add_nodes!(K, L, tree, 1, 1, Np)
    return tree
end

function generate_sample(L::Int, π::Array{Float64}, Np::Int)::Array{DRMSMIP.Sample}
    # generates random samples following a lognormal distribution
    ret = Array{DRMSMIP.Sample}(undef, Np)
    for ii in 1:Np
        ξ = Array{Float64}(undef, L)
        for l in 1:L
            sig = sqrt( log( 0.5+sqrt( ( 0.03*l )^2+0.25 ) ) )
            rnd = sig * randn(rng) .+ log(π[l])
            ξ[l] = exp(rnd)
        end
        ret[ii] = DRMSMIP.Sample(ξ, 1/Np)
    end
    return ret
end

function add_nodes!(K::Int, L::Int, tree::DRMSMIP.DR_Tree, id::Int, k::Int, Np::Int)
    if k < K-1
        ls = iterlist(L,tree.nodes[id].ξ)
        for π in ls
            π_samp = generate_sample(L, π, Np)
            set = DRMSMIP.WassersteinSet(π_samp, 1.0, DRMSMIP.norm_L1)
            cost = zeros(L+1)
            DD.addchild!(tree, id, π, set, cost)
            childid = length(tree.nodes)
            add_nodes!(K, L, tree, childid, k+1, Np)
        end
    elseif k == K-1
        ls = iterlist(L,tree.nodes[id].ξ)
        for π in ls
            cost = vcat(-π, [-1])
            DD.addchild!(tree, id, π, nothing, cost)
        end
    end
end

function iterlist(L::Int, π::Array{Float64})::Array{Array{Float64}}
    # generates all combinations of up and down scenarios
    ls = [Float64[] for _ in 1:2^L]
    ii = 1

    function foo(L::Int, l::Int, arr::Vector{Float64})
        up = (1.0 + 0.03 * l) * π[l]
        dn = (1.0 - 0.03 * l) * π[l]

        if l < L
            arr1 = copy(arr)
            arr1[l] = up
            foo(L, l+1, arr1)

            arr2 = copy(arr)
            arr2[l] = dn
            foo(L, l+1, arr2)
        else
            arr1 = copy(arr)
            arr1[l] = up
            ls[ii] = arr1
            ii+=1

            arr2 = copy(arr)
            arr2[l] = dn
            ls[ii] = arr2
            ii+=1
        end
    end

    foo(L, 1, Array{Float64}(undef, L))
    return ls
end

"""
    max     B_K+∑_{l=1}^{L}π_{K,l}y_{K,l}

    s.t.    B_1+∑_{l=1}^{L}π_{1,l}x_{1,l} = b_1

            b_k+(1+a)B_{k-1}+∑_{l=1}^{L}ρ_{k,l}y_{k-1,l} = B_k+∑_{l=1}^{L}π_{k,l}x_{k,l}, ∀ k=2,…,K

            y_{1,l} = x_{1,l}, ∀ l=1,…,L

            y_{k-1,l}+x_{k,l} = y_{k,l}, ∀ k=2,…,K, l=1,…,L

            x_{k,l} ∈ ℤ , ∀ k=1,…,K, l=1,…,L

            y_{k,l} ≥ 0, ∀ k=1,…,K, l=1,…,L

            B_k ≥ 0, ∀ k=1,…,K.
"""
const K = 3
const L = 2
const a = 0.01
const b1 = 100
const b2 = 30
const Np = 10

function create_scenario_model(K::Int, L::Int, tree::DRMSMIP.DR_Tree, id::Int)
    hist = DD.get_history(tree, id)
    m = Model(GLPK.Optimizer)
    @variable(m, x[1:K,1:L], integer=true)
    #@variable(m, x[1:K,1:L])
    @variable(m, y[1:K,1:L]>=0)
    @variable(m, B[1:K]>=0)

    π = tree.nodes[1].ξ

    con = @constraint(m, B[1] + sum( π[l] * x[1,l] for l in 1:L) == b1)
    set_name(con, "con[1]")

    for l in 1:L
        bal = @constraint(m, y[1,l]-x[1,l]==0)
        set_name(bal, "bal[1,$(l)]")
    end

    for k = 2:K
        π = tree.nodes[hist[k]].ξ
        ρ = tree.nodes[hist[k-1]].ξ * 0.05

        con = @constraint(m, B[k] + sum( π[l] * x[k,l] - ρ[l] * y[k-1,l] for l in 1:L)
            - (1+a) * B[k-1] == b2)
        set_name(con, "con[$(k)]")
        for l in 1:L
            bal = @constraint(m, y[k,l]-x[k,l]-y[k-1,l]==0)
            set_name(bal, "bal[$(k),$(l)]")
        end
    end

    for l in 1:L
        @constraint(m, x[K,l]==0)
    end
    #π = tree.nodes[id].ξ
    #@objective(m, Max, B[K] + sum( π[l] * x[K,l] for l in 1:L )
    @objective(m, Min, 0 )
    return m
end

function main_comp()
    tree = create_tree(K,L,Np)

    #DEmodel = det_eq(L, tree)
    #det_eq_results(tree, DEmodel)

    dual_decomp(L, tree)
end

function dual_decomp(L::Int, tree::DRMSMIP.DR_Tree)

    # Initialize MPI
    DD.parallel.init()

    # Create DualDecomposition instance.
    #algo = DRMSMIP.DRMS_LagrangeDual(tree, BM.TrustRegionMethod)
    algo = DRMSMIP.DRMS_LagrangeDual(tree)

    # Add Lagrange dual problem for each scenario s.
    nodelist = DD.get_stage_id(tree)
    leafdict = DD.leaf2block(nodelist[K])

    # partition scenarios into processes
    DD.parallel.partition(length(nodelist[K]))

    models = Dict{Int,JuMP.Model}(nodelist[K][s] => create_scenario_model(K,L,tree,nodelist[K][s]) for s in DD.parallel.getpartition())
    for s in DD.parallel.getpartition()
        id = nodelist[K][s]
        DD.add_block_model!(algo, s, models[id])
    end

    coupling_variables = Vector{DD.CouplingVariableRef}()

    for s in DD.parallel.getpartition()
        id = nodelist[K][s]
        model = models[id]
        hist = DD.get_history(tree, id)
        yref = model[:y]
        Bref = model[:B]
        for k in 1:K-1
            root = hist[k]
            for l in 1:L
                push!(coupling_variables, DD.CouplingVariableRef(leafdict[id], [root, l], yref[k, l]))
            end
            push!(coupling_variables, DD.CouplingVariableRef(leafdict[id], [root, L+1], Bref[k]))
        end
        # dummy coupling variables
        for l in 1:L
            push!(coupling_variables, DD.CouplingVariableRef(s, [id, l], yref[K, l]))
        end
        Bref = model[:B]
        push!(coupling_variables, DD.CouplingVariableRef(s, [id, L+1], Bref[K]))
    end

    # Set nonanticipativity variables as an array of symbols.
    DD.set_coupling_variables!(algo, coupling_variables)

    bundle_init = DRMSMIP.initialize_bundle(tree, algo)

    # Solve the problem with the solver; this solver is for the underlying bundle method.
    DD.run!(algo, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0), bundle_init)
    return algo
end

function det_eq(L::Int, tree::DRMSMIP.DR_Tree)
    m = Model(GLPK.Optimizer)
    lenN = length(tree.nodes)
    node = tree.nodes[1]
    @variable(m, x[1:lenN,1:L], integer=true)
    #@variable(m, x[1:lenN,1:L])
    @variable(m, y[1:lenN,1:L]>=0)
    @variable(m, B[1:lenN]>=0)
    @variable(m, α[1:lenN]>=0)
    @variable(m, β[1:lenN,1:node.set.N])

    @objective(m, Min, sum( node.cost[l]*y[1,l] for l in 1:L) + node.cost[L+1]*B[1]
         + node.set.ϵ*α[1] + sum( node.set.samples[ss].p*β[1,ss] for ss in 1:node.set.N) )

    π = tree.nodes[1].ξ
    @constraint(m, B[1] + sum( π[l] * x[1,l] for l in 1:L) == b1)
    for l in 1:L
        @constraint(m, y[1,l]-x[1,l]==0)
    end

    function iterate_children(id::Int)
        pnode = tree.nodes[id]
        children = DD.get_children(tree, id)

        for child in children
            cnode = tree.nodes[child]

            π = cnode.ξ
            ρ = pnode.ξ * 0.05

            @constraint(m, B[child] + sum( π[l] * x[child,l] - ρ[l] * y[id,l] for l in 1:L)
                - (1+a) * B[id] == b2)
            for l in 1:L
                @constraint(m, y[child,l]-x[child,l]-y[id,l]==0)
            end

            if length(DD.get_children(tree, child)) > 0
                for s in 1:pnode.set.N
                    foo = pnode.set.norm_func(cnode.ξ, pnode.set.samples[s].ξ)
                    @constraint(m, foo*α[id] + β[id,s]
                        - sum( cnode.cost[l]*y[child,l] for l in 1:L) - cnode.cost[L+1]*B[child] - cnode.set.ϵ*α[child]
                        - sum( cnode.set.samples[ss].p*β[child,ss] for ss in 1:cnode.set.N) >= 0)
                end
            else
                for s in 1:pnode.set.N
                    foo = pnode.set.norm_func(cnode.ξ, pnode.set.samples[s].ξ)
                    @constraint(m, foo*α[id] + β[id,s]
                        - sum( cnode.cost[l]*y[child,l] for l in 1:L) - cnode.cost[L+1]*B[child] >= 0)
                end
            end
            iterate_children(child)
        end

    end
    iterate_children(1)
    nodelist = DD.get_stage_id(tree)
    for id in nodelist[K]
        for l in 1:L
            @constraint(m, x[id,l]==0)
        end
    end
    JuMP.optimize!(m)
    return m
end

function det_eq_results(tree::DRMSMIP.DR_Tree, model::Model)
    open("examples/investment_results/det_eq.lp", "w") do f
        print(f, model)
    end
    nodelist = DD.get_stage_id(tree)

    xref = value.(model[:x])
    Bref = value.(model[:B])
    yref = value.(model[:y])
    αref = value.(model[:α])
    βref = value.(model[:β])

    open("examples/investment_results/det_eq_results.csv", "w") do io

        for leaf in nodelist[K]

            hist = DD.get_history(tree, leaf)
            write(io, "stage, ") + sum(write(io, "$(k), ") for k in 1:K) + write(io, "\n")
            write(io, "B, ") + sum(write(io, string(Bref[id])*", ") for id in hist) + write(io, "\n")
            for l in 1:L
                write(io, "π[$(l)], ") + sum(write(io, string(tree.nodes[id].ξ[l])*", ") for id in hist) + write(io, "\n")
                write(io, "x[$(l)], ") + sum(write(io, string(xref[id,l])*", ") for id in hist) + write(io, "\n")
                write(io, "y[$(l)], ") + sum(write(io, string(yref[id,l])*", ") for id in hist) + write(io, "\n")
            end

            tot = 0
            for id in hist
                tot += sum( tree.nodes[id].cost[l]*yref[id,l] for l in 1:L) + tree.nodes[id].cost[L+1]*Bref[id]
            end
            write(io, "total, " * string(-tot) * "\n")
            write(io, "\n")
        end
        write(io, "objective, "*string(-objective_value(model)))
    end

end

main_comp()


text:
export
point_contained,
point_contained!

"""
point_contained(simplex::ST, point::Vector{VT}) where {T <: AbstractSimplex, VT}

Check if `point` is contained in the interior of `simplex`.

The `point` is contained within `simplex` if it can be expressed as a
convex combination of the vertices of `simplex`. To check this condition,
we iteratively check the signed volume of a modified simplex. Specifically,
at the i-th iteration step, substitute the i-th vertex of `simplex` with
`point` and take the signed volume when when one vertex is substituted.
In the next step, replace the substitated vertex by the original, then
replace the (i + 1)th vertex by `point` and repeat the process.
If any two consecutive signs are different, the point is not contained.

When the point is very close to the edge of a simplex, numerical accuracies
will play a role, so the function sometimes (very rarely) puts `point` inside
`simplex`.
"""
function point_contained(simplex::T, point::Vector{VT}) where {T <: AbstractSimplex, VT}
    dim = dimension(simplex)

    # The signs of the linear combination coefficients
    coefficient_signs = zeros(Float64, dim + 1)

    # A temporary simplex into which we substitute `point`
    temp_simplex = MutableSimplex([zeros(Float64, dim) for i = 1:(dim + 1)])

    # Using the temporary simplex, replace the the first vertex by `point`
    # and fill the remaining vertices with the original `simplex` vertices.
    temp_simplex[1] = point
    [temp_simplex[i] = simplex[i] for i = 2:(dim + 1)]

    # The signed volume of the simplex
    coefficient_signs[1] = sign(orientation(temp_simplex))

    # Iteratively replace the other vertices and see if the signed
    # volume changes. It it does, the point can not be expressed as a
    # convex combination of the vertices of `simplex`.
    for κ = 2:dim

        # Replace the i-th vertex of the temporary simplex with `point`
        temp_simplex[κ] = point

        # Replace the remaining vertices of the temporary simplex
        # with the original vertices
        remaining_vertex_idxs = setdiff(1:(dim + 1), κ)
        [temp_simplex[i] = simplex[i] for i in remaining_vertex_idxs]

        # Check the signed volume and stop if it has changed.
        coefficient_signs[κ] = sign(orientation(temp_simplex))

        if !(coefficient_signs[κ-1] == coefficient_signs[κ])
            return false
        end
    end

    # Check when replacing the last point.
    temp_simplex[end] = point
    [temp_simplex[i] = simplex[i] for i = 1:dim]

    coefficient_signs[end] = sign(orientation(temp_simplex))

    if !(coefficient_signs[end-1] == coefficient_signs[end])
        return false
    else
        return true
    end
end

"""
    point_contained!(temp_simplex::MutableSimplex, simplex::ST, point::Vector{VT}) where {T <: AbstractSimplex, VT}

Check if `point` is contained in the interior of `simplex`.

The `point` is contained within `simplex` if it can be expressed as a
convex combination of the vertices of `simplex`. To check this condition,
we iteratively check the signed volume of a modified simplex. Specifically,
at the i-th iteration step, substitute the i-th vertex of `simplex` with
`point` and take the signed volume when when one vertex is substituted.
In the next step, replace the substitated vertex by the original, then
replace the (i + 1)th vertex by `point` and repeat the process.
If any two consecutive signs are different, the point is not contained.

When the point is very close to the edge of a simplex, numerical accuracies
will play a role, so the function sometimes (very rarely) puts `point` inside
`simplex`.
"""
function point_contained!(temp_simplex::MutableSimplex, simplex::T, point::Vector{VT}) where {T <: AbstractSimplex, VT}
    dim = dimension(simplex)

    # The signs of the linear combination coefficients
    coefficient_signs = zeros(Float64, dim + 1)

    # Using the temporary simplex, replace the the first vertex by `point`
    # and fill the remaining vertices with the original `simplex` vertices.
    temp_simplex[1] = point
    [temp_simplex[i] = simplex[i] for i = 2:(dim + 1)]

    # The signed volume of the simplex
    coefficient_signs[1] = sign(orientation(temp_simplex))

    # Iteratively replace the other vertices and see if the signed
    # volume changes. It it does, the point can not be expressed as a
    # convex combination of the vertices of `simplex`.
    for κ = 2:dim

        # Replace the i-th vertex of the temporary simplex with `point`
        temp_simplex[κ] = point

        # Replace the remaining vertices of the temporary simplex
        # with the original vertices
        remaining_vertex_idxs = setdiff(1:(dim + 1), κ)
        [temp_simplex[i] = simplex[i] for i in remaining_vertex_idxs]

        # Check the signed volume and stop if it has changed.
        coefficient_signs[κ] = sign(orientation(temp_simplex))

        if !(coefficient_signs[κ-1] == coefficient_signs[κ])
            return false
        end
    end

    # Check when replacing the last point.
    temp_simplex[end] = point
    [temp_simplex[i] = simplex[i] for i = 1:dim]

    coefficient_signs[end] = sign(orientation(temp_simplex))

    if !(coefficient_signs[end-1] == coefficient_signs[end])
        return false
    else
        return true
    end
end


text:
@userplot PlotToroidalMomentumError
@recipe function f(p::PlotToroidalMomentumError; nplot=1, latex=true)
    if length(p.args) == 2 && typeof(p.args[1]) <: TimeSeries && typeof(p.args[2]) <: DataSeries
        t  = p.args[1]
        pφ = p.args[2]
        @assert t.n == pφ.nt
    else
        error("Toroidal momentum error plot should be given a timeseries and a data series. Got: $(typeof(p.args))")
    end

    legend := :none
    size := (800,400)

    @series begin
        if latex
            xlabel := L"t"
            ylabel := L"[P_\varphi (t) - P_\varphi (0)] / P_\varphi (0)"
        else
            xlabel := "t"
            ylabel := "[P(t) - P(0)] / P(0)"
        end
        xlims  := (t[0], Inf)
        yformatter := :scientific
        guidefont := font(18)
        tickfont := font(12)
        right_margin := 10mm
        t[0:nplot:end], pφ[0:nplot:end]
    end
end


text:
#### Subsidence model
export AbstractSubsidence,
       NoSubsidence,
       ConstantSubsidence,
       subsidence_velocity

abstract type AbstractSubsidence{FT<:AbstractFloat} end

struct NoSubsidence{FT} <: AbstractSubsidence{FT} end

struct ConstantSubsidence{FT} <: AbstractSubsidence{FT}
  D::FT
end

subsidence_velocity(::NoSubsidence{FT}, z::FT) where {FT} = FT(0)
subsidence_velocity(subsidence::ConstantSubsidence{FT}, z::FT) where {FT} = subsidence.D*z

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
# tests from chapter 13 of the draft of "Speech and Language
# Processing" 3rd edition by Jurafsky & Martin

@testset "J&M 3ed Ch 13" begin

    sent_f13_7 = ["book", "me", "the", "morning", "flight"]

    # leftovers:
    sent = sent_f13_7

    word2id = Dict(word => id for (id, word) in enumerate(sent))
    word2id["root"] = 0

    words(ids) = [(id == 0 ? "root" : sent[id]) for id in ids]

    @testset "Figure 13.7" begin
        @test projective_only(ArcStandard())

        # figure 13.7 in jurafsky & martin SLP 3rd ed., aug 2018 draft
        table = [
            # stack                                    buffer                                  # transition
            (["ROOT"],                                 ["book","me","the","morning","flight"], Shift()),
            (["ROOT","book"],                          ["me","the","morning","flight"],        Shift()),
            (["ROOT","book","me"],                     ["the","morning","flight"],             RightArc("indobj")),
            (["ROOT","book"],                          ["the","morning","flight"],             Shift()),
            (["ROOT","book","the"],                    ["morning","flight"],                   Shift()),
            (["ROOT","book","the","morning"],          ["flight"],                             Shift()),
            (["ROOT","book","the","morning","flight"], [],                                     LeftArc("adv")),
            (["ROOT","book","the","flight"],           [],                                     LeftArc("dt")),
            (["ROOT","book","flight"],                 [],                                     RightArc("dobj")),
            (["ROOT","book",],                         [],                                     RightArc("pred")),
            (["ROOT",],                                [],                                     nothing)]
        sent = copy(sent_f13_7)

        gold_tree = test_sentence("morningflight.conll")

        for arc in (untyped, typed)
            o = cfg -> static_oracle(cfg, gold_tree, arc)
            cfg = initconfig(ArcStandard(), gold_tree)
            for (stk, buf, t) in table
                @test form.(tokens(cfg, stack(cfg))) == stk
                @test form.(tokens(cfg, buffer(cfg))) == buf
                t == nothing && break
                if arc == untyped
                    t isa LeftArc && (t = LeftArc())
                    t isa RightArc && (t = RightArc())
                end
                @test o(cfg) == t
                cfg = t(cfg)
            end
        end
    end

    @testset "Figure 13.8" begin
        @test projective_only(ArcStandard())

        # figure 13.8 in jurafsky & martin SLP 3rd ed., aug 2018 draft
        table = [
            # stack                                    buffer                                  # transition
            (["ROOT"],                                 ["book","the","flight","through","houston"], Shift()),
            (["ROOT","book"],                          ["the","flight","through","houston"],   Shift()),
            (["ROOT","book","the"],                    ["flight","through","houston"],         Shift()),
            (["ROOT","book","the","flight"],           ["through","houston"],                  LeftArc("det")),
            (["ROOT","book","flight"],                 ["through","houston"],                  Shift()),
            (["ROOT","book","flight","through"],       ["houston"],                            Shift()),
            (["ROOT","book","flight","through","houston"], [],                                 LeftArc("case")),
            (["ROOT","book","flight","houston"],       [],                                     RightArc("nmod")),
            (["ROOT","book","flight"],                 [],                                     RightArc("dobj")),
            (["ROOT","book",],                         [],                                     RightArc("root")),
            (["ROOT",],                                [],                                     nothing)]
        sent = copy(sent_f13_7)

        gold_tree = test_sentence("flightthroughhouston.conll")

        for arc in (untyped, typed)
            o = cfg -> static_oracle(cfg, gold_tree, arc)
            cfg = initconfig(ArcStandard(), gold_tree)
            for (stk, buf, t) in table
                @test form.(tokens(cfg, stack(cfg))) == stk
                @test form.(tokens(cfg, buffer(cfg))) == buf
                t == nothing && break
                if arc == untyped
                    t isa LeftArc && (t = LeftArc())
                    t isa RightArc && (t = RightArc())
                    @test endswith(showstr(t), "()")
                end
                @test o(cfg) == t
                cfg = t(cfg)
            end
        end
    end

    @testset "Figure 13.10" begin

        @test DependencyTrees.projective_only(ArcEager())

        table = [
            # stack                                    buffer                                  # transition
            (["ROOT"],                                 ["book","the","flight","through","houston"], RightArc("root")),
            (["ROOT","book"],                          ["the","flight","through","houston"],   Shift()),
            (["ROOT","book","the"],                    ["flight","through","houston"],         LeftArc("det")),
            (["ROOT","book"],                          ["flight","through","houston"],         RightArc("dobj")),
            (["ROOT","book","flight"],                 ["through","houston"],                  Shift()),
            (["ROOT","book","flight","through"],       ["houston"],                            LeftArc("case")),
            (["ROOT","book","flight"] ,                ["houston"],                            RightArc("nmod")),
            (["ROOT","book","flight","houston"],       [],                                     Reduce()),
            (["ROOT","book","flight"],                 [],                                     Reduce()),
            (["ROOT","book"],                          [],                                     Reduce()),
            (["ROOT"],                                 [],                                     nothing)]

        gold_tree = test_sentence("flightthroughhouston.conll")

        for arc in (untyped, typed)
            o = cfg -> static_oracle(cfg, gold_tree, arc)
            cfg = initconfig(ArcEager(), gold_tree)
            for (stk, buf, t) in table
                @test form.(tokens(cfg, stack(cfg))) == stk
                @test form.(tokens(cfg, buffer(cfg))) == buf
                t == nothing && break
                if arc == untyped
                    t isa LeftArc && (t = LeftArc())
                    t isa RightArc && (t = RightArc())
                    @test endswith(showstr(t), "()")
                end
                @test o(cfg) == t
                if o(cfg) != t
                    @show stk buf t cfg o(cfg)
                end
                cfg = t(cfg)
            end
            @test isfinal(cfg)
        end
    end
end


text:
"""
    ```prim = primitive_cell(xtal)```

Returns the minimal representation (primitive unit cell) of a crystal structure.
"""
function primitive_cell(xtal::Crystal)
    # check for pymatgen.io.cif dependency
    if isnothing(rc[:pymatgen])
        error("Python dependency pymatgen not loaded.")
    else
        pymatgen = rc[:pymatgen]
    end
    tempfile = joinpath(pwd(), ".temp_$(uuid1()).cif")
    # copy out xtal and convert it to primitive cell
    write_cif(xtal, tempfile)
    pymatgen.CifParser(tempfile).get_structures()[1].get_primitive_structure().to(filename=tempfile)
    # load the result
    primitive = Crystal(tempfile)
    # clean up
    rm(tempfile)
    name = split(xtal.name, ".cif")[1] * "_primitive_cell.cif"
    return Crystal(name, primitive.box, primitive.atoms, primitive.charges)
end

max_stars_repo_name: SimonEnsemble/XtalsPyTools

max_issues_repo_name: SimonEnsemble/XtalsPyTools


text:
module Univariates

#パッケージの読み込み
using DataFrames,StatsBase#,Plots,StatsPlots,GR
include("univariate.jl")
include("groupbycolunivariate.jl")

export numericsummary,summarymerge,univariate
export groupbycolnumericsummary,groupbycolnonnumericsummary,groupbycolunivariate
end # module


text:
struct PTDF{Ax,L<:NTuple{2,Dict}} <: PowerNetworkMatrix{Float64}
    data::Array{Float64,2}
    axes::Ax
    lookup::L
end

function _buildptdf(branches, nodes, dist_slack::Array{Float64}=[0.1])

    buscount = length(nodes)
    linecount = length(branches)
    num_bus = Dict{Int32,Int32}()

    for (ix,b) in enumerate(nodes)
        num_bus[get_number(b)] = ix
    end

    A = zeros(Float64,buscount,linecount);
    inv_X = zeros(Float64,linecount,linecount);

   #build incidence matrix
   #incidence_matrix = A

    for (ix,b) in enumerate(branches)

        if isa(b,DCBranch)
            @warn("PTDF construction ignores DC-Lines")
            continue
        end

        A[num_bus[get_arc(b) |> get_from |> get_number], ix] =  1;

        A[num_bus[get_arc(b) |> get_to |> get_number], ix] = -1;

        if isa(b,Transformer2W)
            inv_X[ix,ix] = 1/get_x(b);

        elseif isa(b,TapTransformer)
            inv_X[ix,ix] = 1/(get_x(b)*get_tap(b));

        elseif isa(b, Line)
            inv_X[ix,ix] = 1/get_x(b);

        elseif isa(b,PhaseShiftingTransformer)
            y = 1 / (get_r(b) + get_x(b) * 1im)
            y_a = y / (get_tap(b) * exp(get_α(b) * 1im * (π / 180)))
            inv_X[ix,ix] = 1/imag(y_a)
        end

    end

    slacks = [num_bus[get_number(n)] for n in nodes if get_bustype(n) == REF::BusType]
    slack_position = slacks[1]
    B = gemm('N','T', gemm('N','N',A[setdiff(1:end, slack_position),1:end] ,inv_X), A[setdiff(1:end, slack_position),1:end])
    if dist_slack[1] == 0.1 && length(dist_slack) ==1
        (B, bipiv, binfo) = getrf!(B)
        S_ = gemm('N','N', gemm('N','T', inv_X, A[setdiff(1:end, slack_position), :]), getri!(B, bipiv) )
        S = hcat(S_[:,1:slack_position-1],zeros(linecount,),S_[:,slack_position:end])

    elseif dist_slack[1] != 0.1 && length(dist_slack)  == buscount
        @info "Distributed bus"
        (B, bipiv, binfo) = getrf!(B)
        S_ = gemm('N','N', gemm('N','T', inv_X, A[setdiff(1:end, slack_position), :]), getri!(B, bipiv) )
        S = hcat(S_[:,1:slack_position-1],zeros(linecount,),S_[:,slack_position:end])
        slack_array =dist_slack/sum(dist_slack)
        slack_array = reshape(slack_array,buscount,1)
        S = S - gemm('N','N',gemm('N','N',S,slack_array),ones(1,buscount))

    elseif length(slack_position) == 0
        @warn("Slack bus not identified in the Bus/Nodes list, can't build PTDF")
        S = Array{Float64,2}(undef,linecount,buscount)
    end

    return  S, A

end

function PTDF(branches, nodes, dist_slack::Array{Float64}=[0.1])

    #Get axis names
    line_ax = [get_name(branch) for branch in branches]
    bus_ax = [get_number(bus) for bus in nodes]
    S, A = _buildptdf(branches, nodes, dist_slack)

    axes = (line_ax, bus_ax)
    look_up = (_make_ax_ref(line_ax),_make_ax_ref(bus_ax))

    return PTDF(S, axes, look_up)

end

function PTDF(sys::System, dist_slack::Array{Float64}=[0.1])
    branches = get_components(ACBranch, sys)
    nodes = get_components(Bus, sys)

    return PTDF(branches, nodes, dist_slack)

end


text:
#Read d00_intro.jl

#Define bit manipulation functions
include("f01_init_file.jl")

"""
So now we have the relevant functions to work with bits. Now onto the format.

From 3.2.3
1. Bit 1 - BFINAL - is this the last block?
2. Bits 2,3 - BTYPE - specifies how the data are compressed, as follows:

    00 - no compression
    01 - compressed with fixed Huffman codes
    10 - compressed with dynamic Huffman codes
    11 - reserved (error)

So, let's read our file.
"""

BFINAL = pop_bits(1) #BFINAL = 1 - this is the final block
BTYPE  = pop_bits(2) #BTYPE  = 3 - Dynamic Huffman Coding

"""
Since BTYPE = 3 we will be using dynamic Huffman coding

we need to look at 3.2.7 of the RFC which states:
    5 Bits: HLIT, # of Literal/Length codes - 257 (257 - 286)
    5 Bits: HDIST, # of Distance codes - 1        (1 - 32)
    4 Bits: HCLEN, # of Code Length codes - 4     (4 - 19)
"""

HLIT = pop_bits(5)
HDIST = pop_bits(5)
HCLEN = pop_bits(4)

n_literal_codes = HLIT + 257
n_distance_codes = HDIST + 1
n_code_length_codes = HCLEN + 4

MAX_NUMBER_CODE_LENGTH_CODES = 19

code_lengths = zeros(UInt8, MAX_NUMBER_CODE_LENGTH_CODES)

code_lengths_order = UInt8[16,17,18,0,8,7,9,6,10,5,11,4,12,3,13,2,14,1,15]
#Convert from 1 based to zero based
code_lengths_order += 1

length_counts = zeros(UInt8, 8)

for i = 1:n_code_length_codes
    index = code_lengths_order[i]
    length_value = pop_bits(3)
    code_lengths[index] = length_value

    #We'll use this later
    length_counts[length_value+1] += 1
end

for i = 1:MAX_NUMBER_CODE_LENGTH_CODES
    @printf("%d: %d\n",i-1,code_lengths[i])
end

#d03_code_length_results.jl
#-------------------------------------

#This is not the best approach but it works

code_values = zeros(UInt8, MAX_NUMBER_CODE_LENGTH_CODES)

last_value = -1
last_length = -1
#Must loop over lengths on outer, not inner
for current_length = 1:7
    for i = 1:MAX_NUMBER_CODE_LENGTH_CODES
    #We skip 0, it doesn't take any value
    #7 is the max because we used 3 bytes (encoding 0-7)
        if code_lengths[i] == current_length
            #Then we will assign the code some value
            if last_value == -1
                code_values[i] = 0
                last_value = 0
            elseif last_length != current_length
                last_length = current_length
                code_values[i] = 2*(last_value+1)
                last_value = code_values[i]
            else
                code_values[i] = last_value+1
                last_value += 1
            end
        end
    end
end

#summary
for i = 1:MAX_NUMBER_CODE_LENGTH_CODES
    if code_lengths[i] != 0
        len = Int64(code_lengths[i])
        bin_value = bin(code_values[i],len)
        pad_string = lpad(bin_value,7)
        #Subtract 1 due to IDs typically being 0 based
        index = lpad(i-1,2)
        #@printf("bin %s\n",bin_value)
        #@printf("len %d\n",7-len)
        #@printf("lpad %s\n",lpad(bin_value,7-len))
        @printf("%s, %s\n",index,pad_string)
    end
end

#=
This yields:
 0,     010
 4,    1100
 5,    1101
 6,     011
 7,      00
 8,     100
 9,     101
10,    1110
11,   11110
16,  111110
17, 1111110
18, 1111111

Note that I've right aligned the bit strings.

Normally when lookings at bits, we consider the first bit
to be the bit on the right. Then we move to the left.

Now consider the bit string "101"
If we see that, which ID does it represent?

You might say #9, but what about #5?
=#

#The same code as above, but now right padded
for i = 1:MAX_NUMBER_CODE_LENGTH_CODES
    if code_lengths[i] != 0
        len = Int64(code_lengths[i])
        bin_value = bin(code_values[i],len)
        pad_string = rpad(bin_value,7)
        #Subtract 1 due to IDs typically being 0 based
        index = lpad(i-1,2)
        #@printf("bin %s\n",bin_value)
        #@printf("len %d\n",7-len)
        #@printf("lpad %s\n",lpad(bin_value,7-len))
        @printf("%s, %s\n",index,pad_string)
    end
end

#d04_I1_lengths_part1.jl

#Sort index by length
indices_to_process = zeros(UInt8,n_code_length_codes)

MAX_LENGTH = 8
current_start_indices = zeros(UInt8, MAX_LENGTH)
current_start_indices[1] = 1
for i = 1:MAX_LENGTH-1
    current_start_indices[i+1] = current_start_indices[i] + length_counts[i]
end

for i = 1:MAX_NUMBER_CODE_LENGTH_CODES
    if code_lengths[i] != 0
        len_b1 = code_lengths[i]+1
        current_index = current_start_indices[len_b1]
        @printf("%d of len %d to %d\n",i,code_lengths[i],current_index)
        indices_to_process[current_index] = i
        current_start_indices[len_b1] += 1
    end
end

#indices_to_process
#
#Now we know the order in which to process the values

code_values2 = zeros(UInt8,MAX_NUMBER_CODE_LENGTH_CODES)

current_value = UInt16(0)
last_length = -1
for i = 2:n_code_length_codes
    current_index = indices_to_process[i]
    current_length = code_lengths[current_index]

    #x2 - does nothing since we are going in reverse
    #add 1
    #   - find first zero
    #   - everything beyond needs to be zerod
    #   - add 1

    #001 => 1 << 0
    #010 => 1 << 1 (length = 2)
    #100 => 1 << 2 (length = 3)

    if current_length == last_length
        #Add 1 at current length
        bit = UInt16(1) << (current_length-1)
    else
        #add 1 at old length then multiply by 2
        #
        # when reversed increasing multiplying
        # by 2 has no effect
        #
        # =>  01
        # => 001 #Now "3" bits but still same value
        bit = UInt16(1) << (current_length-2)
    end

    #Adding 1
    #---------------------------------
    while (current_value & bit != 0)
        bit >>= 1
    end
    #Keep all lower bits
    current_value &= bit-1
    #Add in current bit
    current_value |= bit

    last_length = current_length
    code_values2[current_index] = current_value
end

#Building lookup table
i1_table = zeros(UInt16,256)
for i = 1:MAX_NUMBER_CODE_LENGTH_CODES
#code_values2
    current_value = UInt16(code_values2[i])
    current_length = code_lengths[i]
    if current_length != 0
        #5
        @printf("%d\n",i)
        increment = UInt16(1) << current_length
        while current_value < 256
            #TODO: problem if bit string is 0
            #Unfortunately we'll have to test if zero
            #and if not, then do our indexing
            if current_value != 0
            i1_table[current_value] = current_value
            end
            current_value += increment
        end
    end
end


text:
(pwd() != @__DIR__) && cd(@__DIR__) # allow starting app from bin/ dir

using BitempWebApp
push!(Base.modules_warned_for, Base.PkgId(BitempWebApp))
BitempWebApp.main()

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
"""
hausleitner.jl

Bond order tight-binding approach by Hausleitner and Hafner.

References:

* C. Hausleitner and J. Hafner: Phys. Rev. B, 45 (1992) 115-127.
"""

immutable HHTB <: TBPerturbation
  N::Int # number of components
  Z::Int # total coordination number on the Bethe lattice
  σ::Float64 # ordering parameter
  x::Vector{Float64} # composition
  p::Matrix{Float64} # probability
  Nd::Vector{Float64} # number of d-electrons
  Ed::Vector{Float64} # d-state energy
  Wd::Vector{Float64} # band width
  r₀::Vector{Float64} # atomic radius
end

function HHTB(Nd::Real, Ed::Real, Wd::Real, r₀::Real; Z::Int = 12)
  HHTB(1, Z, 0, ones(1), ones(1,1), [Nd], [Ed], [Wd], [r₀])
end

function HHTB(x::Vector{Float64}, Nd::Vector{Float64}, Ed::Vector{Float64}, Wd::Vector{Float64}, r₀::Vector{Float64}; Z::Int = 12, σ::Real = 0)
  N = length(x)
  p = [ x[j] for i in 1:N, j in 1:N ]
  HHTB(N, Z, σ, x, p, Nd, Ed, Wd, r₀)
end

ncomp(hhtb::HHTB) = hhtb.N

function coreradius(hhtb)::Matrix{Float64}
  @attach(hhtb, N, r₀)
  [ 0.5 * (r₀[i] + r₀[j]) / 2 for i in 1:N, j in 1:N ]
end

function cutoffradius(hhtb)::Matrix{Float64}
  @attach(hhtb, N, r₀)
  [ 10 * (r₀[i] + r₀[j]) / 2 for i in 1:N, j in 1:N ]
end

function cutoffenergy(hhtb)::Tuple{Float64,Float64}
  @attach(hhtb, Ed, Wd)
  Emin = minimum(Ed) - maximum(Wd)
  Emax = maximum(Ed) + maximum(Wd)
  (Emin, Emax)
end

function nntransferintegral(hhtb::HHTB)::Matrix{Float64}
  @attach(hhtb, Z, N, Wd)
  hd = Matrix{Float64}(N,N)
  for i in 1:N
    hd[i,i] = Wd[i] / 12
  end
  for i in 1:N, j in 1:N
    i == j && continue
    hd[i,j] = √( hd[i,i] * hd[j,j] )
  end
  hd
end

function transferintegral(hhtb::HHTB, a::Real, n::Int)::Matrix{Function}
  @attach(hhtb, N, Wd, r₀)
  ret = Matrix{Function}(N,N)
  η_ddσ = -6 * 2/5
  η_ddπ =  4 * 2/5
  η_ddδ = -1 * 2/5
  for i in 1:N
    h(r) = a * √( 1/5 * (η_ddσ^2 + 2η_ddπ^2 + 2η_ddδ^2) ) * Wd[i] * r₀[i]^5 / r^n
    ret[i,i] = h
  end
  for i in 1:N, j in 1:N
    i == j && continue
    h(r) = √( ret[i,i](r) * ret[j,j](r) )
    ret[i,j] = h
  end
  ret
end

function transferintegral(hhtb::HHTB)::Matrix{Function}
  transferintegral(hhtb, 1, 5)
end

function transferintegral_derivative(hhtb::HHTB)::Matrix{Function}
  transferintegral(hhtb, -5, 6)
end

function transfermatrix(hhtb::HHTB)::Matrix{Function}
  @attach(hhtb, Z, N, x, Ed, p)
  h::Matrix{Float64} = nntransferintegral(hhtb)
  ret = Matrix{Function}(N,N)
  if N == 1 || prod(x) == 0
    α = findfirst(x -> x ≠ 0, x)
    function S(E::Real)
      A₂ = Z*h[α,α]
      A₁ = -(E - Ed[α])
      A₀ = h[α,α]
      roots = Polynomials.roots( Poly([A₀, A₁, A₂]) )
      roots[1]
    end
    ret[α,α] = S
    if N == 2
      β = findfirst(x -> x == 0, x)
      ret[α,β] = E -> h[α,β]/h[α,α]*ret[α,α](E)
      ret[β,α] = E -> 0.
      ret[β,β] = E -> 0.
    end
    return ret
  end
  function S₁₁(E::Real)
    A₄ = (Z-1)^2 * ( p[1,1]^2*p[2,2] - p[1,1]*p[1,2]*p[2,1] )
    A₃ = (Z-1) * ( (p[1,2]*p[2,1] - 2p[1,1]*p[2,2]) * (E - Ed[1]) / h[1,1] + p[1,1]*p[1,2] * (E - Ed[2]) / h[2,2] )
    A₂ = p[2,2] * ( (E - Ed[1]) / h[1,1] )^2 - p[1,2] * (E - Ed[1])*(E - Ed[2]) / ( h[1,1] * h[2,2] ) + (Z-1) * (2*p[1,1]*p[2,2] + p[1,2]^2 - p[1,2]*p[2,1])
    A₁ = p[1,2]*(E - Ed[2])/h[2,2] - 2*p[2,2]*(E - Ed[1])/h[1,1]
    A₀ = p[2,2]
    roots::Vector{Number} = cardano(A₃, A₂, A₁, A₀)
    roots[2]
  end
  S₁₂(E) = h[1,2]/h[1,1]*S₁₁(E)
  S₂₁(E) = (-(Z-1)*p[1,1]*h[1,1]*S₁₁(E)^2 + (E-Ed[1])*S₁₁(E) - h[1,1]) / ((Z-1)*p[1,2]*h[1,2]*S₁₁(E))
  S₂₂(E) = h[2,2]/h[2,1]*S₂₁(E)
  [ S₁₁ S₁₂; S₂₁ S₂₂ ]
end

function selfenergy(hhtb::HHTB)::Vector{Function}
  @attach(hhtb, Z, N, x, p)
  h::Matrix{Float64} = nntransferintegral(hhtb)
  S::Matrix{Function} = transfermatrix(hhtb)
  [ Δ(E) = Z*sum( p[i,j]*h[i,j]*S[j,i](E) for j in 1:N ) for i in 1:N ]
end

function diagonalgreenfunction(hhtb::HHTB)::Vector{Function}
  @attach(hhtb, Z, N, Ed)
  Δ::Vector{Function} = selfenergy(hhtb)
  [ Gᵢᵢ(E) = 1 / (E - Ed[i] - Δ[i](E)) for i in 1:N ]
end

function offdiagonalgreenfunction(hhtb::HHTB)::Matrix{Function}
  @attach(hhtb, Z, N, Ed)
  h::Matrix{Float64} = nntransferintegral(hhtb)
  Δ::Vector{Function} = selfenergy(hhtb)
  S::Matrix{Function} = transfermatrix(hhtb)
  [ Gᵢⱼ(E) = h[i,j] / ( E - Ed[i] - (Z-1)/Z * Δ[i](E) ) / ( E - Ed[j] - Δ[j](E) ) for i in 1:N, j in 1:N ]
end

function partialdensityofstate(hhtb::HHTB)::Vector{Function}
  @attach(hhtb, N, Z, x, Ed)
  G::Vector{Function} = diagonalgreenfunction(hhtb)
  [ D(E) = 10*x[i]*1/π*imag(G[i](E)) for i in 1:N ]
end

function totaldensityofstate(hhtb::HHTB)::Function
  @attach(hhtb, N)
  Dᵢ = partialdensityofstate(hhtb)
  D(E) = sum(Dᵢ[i](E) for i in 1:N)
end

densityofstate(hhtb) = totaldensityofstate(hhtb)

function fermienergy(hhtb::HHTB)::Float64
  @attach(hhtb, N, x, Nd, Ed, Wd)
  D::Function = totaldensityofstate(hhtb)
  Emin, Emax = cutoffenergy(hhtb)
  fopt(E) = abs(sum(x .* Nd) - ∫(D, Emin, E))
  opt = Optim.optimize(fopt, Emin, Emax)
  Ef = Optim.minimizer(opt)
end

function bondorder(hhtb::HHTB)::Matrix{Float64}
  @attach(hhtb, N, Ed, Wd)
  G::Matrix{Function} = offdiagonalgreenfunction(hhtb)
  Ef::Float64 = fermienergy(hhtb)
  Emin, Emax = cutoffenergy(hhtb)
  Θ = Matrix{Float64}(N,N)
  for i in 1:N
    Θ[i,i] = 10 * 2/π * ∫(E -> imag(G[i,i](E)), Emin, Ef)
  end
  for i in 1:N, j in 1:N
    i == j && continue
    Θ[i,j] = 10 * 2/π * 1/2 * ∫(E -> imag(G[i,j](E) + G[j,i](E)), Emin, Ef)
  end
  Θ
end

function pairpotential_bond(hhtb::HHTB)::Matrix{Function}
  @attach(hhtb, N)
  h::Matrix{Function} = transferintegral(hhtb)
  Θ::Matrix{Float64} = bondorder(hhtb)
  [ u_bond(r) = h[i,j](r) * Θ[i,j] for i in 1:N, j in 1:N ]
end

function pairpotential_rep(hhtb::HHTB, a::Real, n::Int)::Matrix{Function}
  @attach(hhtb, N, Nd, Wd, r₀)
  [ u_rep(r) = a * 16/25*√(Nd[i]*Nd[j])*Wd[i]*Wd[j]*r₀[i]^5*r₀[j]^5 / r^n
    for i in 1:N, j in 1:N ]
end

function pairpotential_rep(hhtb::HHTB)::Matrix{Function}
  pairpotential_rep(hhtb, 1, 8)
end

function pairpotential_rep_derivative(hhtb::HHTB)::Matrix{Function}
  pairpotential_rep(hhtb, -8, 9)
end

function pairpotential(hhtb::HHTB)::Matrix{Function}
  @attach(hhtb, N)
  u_bond = pairpotential_bond(hhtb)
  u_rep = pairpotential_rep(hhtb)
  [ u(r) = u_rep[i,j](r) + u_bond[i,j](r) for i in 1:N, j in 1:N ]
end

function pairpotential_derivative(hhtb::HHTB)::Matrix{Function}
  @attach(hhtb, N)
  u′_bond = pairpotential_bond_derivative(hhtb)
  u′_rep = pairpotential_rep_derivative(hhtb)
  [ u′(r) = u′_bond[i,j](r) + u′_rep[i,j](r) for i in 1:N, j in 1:N ]
end

function pairpotential_bond_derivative(hhtb::HHTB)::Matrix{Function}
  @attach(hhtb, N)
  h′::Matrix{Function} = transferintegral_derivative(hhtb)
  Θ::Matrix{Float64} = bondorder(hhtb)
  [ u′_bond(r) = h′[i,j](r) * Θ[i,j] for i in 1:N, j in 1:N ]
end

function chargetransfer(hhtb::HHTB)::Vector{Float64}
  @attach(hhtb, N, x, Nd)
  Emin, Emax = cutoffenergy(hhtb)
  Ef::Float64 = fermienergy(hhtb)
  D::Vector{Function} = partialdensityofstate(hhtb)
  [ ∫(D[i], Emin, Ef) / x[i] - Nd[i] for i in 1:N ]
end

function bandenergy(hhtb::HHTB)::Float64
  D::Function = totaldensityofstate(hhtb)
  Ef::Float64 = fermienergy(hhtb)
  Emin, Emax = cutoffenergy(hhtb)
  E_band = ∫(E -> D(E)*E, Emin, Ef)
end

function bondenergy(hhtb::HHTB)::Float64
  @attach(hhtb, N, Ed)
  D::Vector{Function} = partialdensityofstate(hhtb)
  Ef::Float64 = fermienergy(hhtb)
  Emin, Emax = cutoffenergy(hhtb)
  E_bond = sum( ∫(E -> D[i](E)*(E - Ed[i]), Emin, Ef) for i in 1:N )
end

function onsiteenergy(hhtb::HHTB)::Float64
  @attach(hhtb, N, x, Nd, Ed)
  ΔNd = chargetransfer(hhtb)
  E_site = sum( x[i]*(Nd[i] + ΔNd[i])*Ed[i] for i in 1:N )
end

function internal(hhtb::HHTB, ref::ReferenceSystem)::Float64
  u_rep = pairpotential_rep(hhtb)
  U_band = bandenergy(hhtb)
  U_rep = internal(ref, hhtb, u_rep)
  U = U_band + U_rep
end


text:
module Extents

export Extent, extent, bounds

"""
    Extent

    Extent(; kw...)
    Extent(values::NamedTuple)

A wrapper for a `NamedTuple` of tuples holding
the lower and upper bounds for each dimension of the object.

`keys(extent)` will return the dimension name Symbols,
in the order the dimensions are used in the object.

`values` will return a tuple of tuples: `(lowerbound, upperbound)` for each dimension.
"""
struct Extent{T<:NamedTuple}
    bounds::T
end
Extent(; kw...) = Extent(values(kw))

bounds(ext::Extent) = getfield(ext, :bounds)

function Base.getproperty(ext::Extent, key::Symbol)
    haskey(bounds(ext), key) || throw(ErrorException("Extent has no field $key"))
    getproperty(bounds(ext), key)
end
function Base.getindex(ext::Extent, key::Symbol)
    haskey(bounds(ext), key) || throw(ErrorException("Extent has no field $key"))
    getindex(bounds(ext), key)
end
function Base.getindex(ext::Extent, i::Int)
    haskey(bounds(ext), i) || throw(ErrorException("Extent has no field $i"))
    getindex(bounds(ext), i)
end
Base.keys(ext::Extent{<:NamedTuple}) = keys(bounds(ext))
Base.values(ext::Extent) = values(bounds(ext))

function Base.:(==)(a::Extent{<:NamedTuple{K1}}, b::Extent{<:NamedTuple{K2}}) where {K1, K2}
    length(K1) == length(K2) || return false
    all(map(k -> k in K1, K2)) || return false
    return all(map((k -> a[k] == b[k]), K1))
end
Base.:(==)(a::Extent, b::Extent) = false

function Base.show(io::IO, ::MIME"text/plain", extent::Extent)
    print(io, "Extent")
    show(io, bounds(extent))
end

"""
    extent(x)

Returns an [`Extent`](@ref), holding the bounds for each dimension of the object.
"""
function extent end

extent(extent) = nothing
extent(extent::Extent) = extent

end


text:
export RandomPatch

"""
    RandomPatch()

A random patch. For this we first draw entries of a vector `v` independently from a
complex normal distribution (`randn(ComplexF64)`). And then normalize `v` with respect
to the 2-norm.
"""
struct RandomPatch <: AbstractAffinePatch end

struct RandomPatchState{T,N} <: AbstractAffinePatchState{N}
    v::PVector{T,N}
end

is_global_patch(::RandomPatch) = true

function state(::RandomPatch, x::PVector)
    v = similar(x, ComplexF64)
    Random.randn!(v)
    LinearAlgebra.normalize!(v)
    RandomPatchState(v)
end

onpatch!(x::AbstractVector, state::RandomPatchState) = onpatch!(x, state.v)
evaluate!(u, state::RandomPatchState, x::PVector) = evaluate_patch!(u, state.v, x)
jacobian!(U, state::RandomPatchState, x::PVector) = jacobian_patch!(U, state.v, x)


text:
using GeometricIntegrators.Integrators
using GeometricIntegrators.Tableaus
using GeometricIntegrators.Utils
using Test

include("coefficients_tests.jl")
include("rk_tableaus_tests.jl")
include("splitting_tableaus_tests.jl")


text:
# check the generlised Euler (Stieltjes) numbers
#    this is mostly just a check that code is correct
#    most errors are arising because of errors in the functions being used
using Polylogarithms
using BenchmarkTools
using DataFrames, CSV
using PyPlot

S = Symbol("S_n")
data1 = CSV.read("../data/stieltjes_test_data_1.csv")
m = size(data1,1)
S1 = zeros(Float64, m)
S1_ref = zeros(Float64, m)
error1 = zeros(Float64, m)
rel_error1 = zeros(Float64, m)
for i=1:m
    S1[i] = stieltjes( data1[i,:n] )
    S1_ref[i] = data1[i,S]
    error1[i] = Float64( S1[i] - S1_ref[i]  )
    rel_error1[i] = error1[i]/Float64( S1_ref[i] )
end
println("max abs. error      = $(maximum( abs.(error1) ))")
println("max rel. abs. error = $(maximum( abs.(rel_error1) ))")