text:
const winningpositions = [[1, 2, 3], [4, 5, 6], [7, 8, 9], [1, 4, 7],
    [2, 5, 8], [3, 6, 9],[1, 5, 9], [7, 5, 3]]

function haswon(brd, xoro)
    marked = findall(x -> x == xoro, brd)
    for pos in winningpositions
        if length(pos) <= length(marked) && pos == sort(marked)[1:3]
            return true
        end
    end
    false
end

function readcharwithprompt(prompt, expected)
    ret = '*'
    while !(ret in expected)
        print("\n", prompt, " ->  ")
        ret = lowercase(chomp(readline()))[1]
    end
    ret
end

availablemoves(brd) = findall(x -> x == ' ', brd)
cornersopen(brd) = [x for x in [1, 3, 7, 9] if brd[x] == ' ']
int2char(x) = Char(x + UInt8('0'))
char2int(x) = UInt8(x) - UInt8('0')
getyn(query) = readcharwithprompt(query, ['y', 'n'])
gettheirmove(brd) = char2int(readcharwithprompt("Your move(1-9)", int2char.(availablemoves(brd))))

function findwin(brd, xoro)
    tmpbrd = deepcopy(brd)
    for mv in availablemoves(tmpbrd)
        tmpbrd[mv] = xoro
        if haswon(tmpbrd, xoro)
            return mv
        end
        tmpbrd[mv] = ' '
    end
    return nothing
end

function choosemove(brd, mychar, theirchar)
    if all(x -> x == ' ', brd)
        brd[rand(cornersopen(brd))] = mychar # corner trap if starting game
    elseif availablemoves(brd) == [] # no more moves
        println("Game is over. It was a draw.")
        exit(0)
    elseif (x = findwin(brd, mychar)) != nothing || (x = findwin(brd, theirchar)) != nothing
        brd[x] = mychar # win if possible, block their win otherwise if their win is possible
    elseif brd[5] == ' '
        brd[5] = mychar # take center if open and not doing corner trap
    elseif (corners = cornersopen(brd)) != []
        brd[rand(corners)] = mychar # choose a corner over a side middle move
    else
        brd[rand(availablemoves(brd))] = mychar # random otherwise
    end
end

function display(brd)
    println("+-----------+")
    println("| ", brd[1], " | ", brd[2], " | ", brd[3], " |")
    println("| ", brd[4], " | ", brd[5], " | ", brd[6], " |")
    println("| ", brd[7], " | ", brd[8], " | ", brd[9], " |")
    println("+-----------+")
end

function tictactoe()
    board = fill(' ', 9)
    println("Board move grid:\n 1 2 3\n 4 5 6\n 7 8 9")
    yn = getyn("Would you like to move first (y/n)?")
    if yn == 'y'
        mychar = 'O'
        theirchar = 'X'
        board[gettheirmove(board)] = theirchar
    else
        mychar = 'X'
        theirchar = 'O'
    end
    while true
        choosemove(board, mychar, theirchar)
        println("Computer has moved.")
        display(board)
        if haswon(board, mychar)
            println("Game over. Computer wins!")
            exit(0)
        elseif availablemoves(board) == []
            break
        end
        board[gettheirmove(board)] = theirchar
        println("Player has moved.")
        display(board)
        if haswon(board, theirchar)
            println("Game over. Player wins!")
            exit(0)
        elseif availablemoves(board) == []
            break
        end
    end
    println("Game over. It was a draw.")
end

tictactoe()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using DESincEig, SincFun, Base.Test

# Testing the main function SincEigen.

println("Testing Example 1")
f(x) = (35/4)./x.^2  .- 2 .+ x.^2 /16
(RESULTS, All_Abs_Error_Approx , hoptimal , n, MatrixSizes) = SincEigen( f , ones , SemiInfinite1{Float64}() , [1.5;0.03125] , [1.0;2.0] , pi/4, tol=1e-10 )
@test norm(RESULTS[1:5,3]-collect(0.0:6.0)) < sqrt(eps())

println("Testing quartic anharmonic oscillators")
# Exact values for the ground states from E. J. Weniger, Ann. Phys. 246:133--165, 1996.
V(x,β) = x.^2 .+ β.*x.^4
β = [0.2;1.0;4.0]
exact_values = [1.1182926543696632;1.392351641530291855;1.9031369454587908]

for i=1:length(β)
B = sqrt(β[i])*(1/2)^3 / 3.0
    (RESULTS, All_Abs_Error_Approx , hoptimal , n, MatrixSizes) = SincEigen( x->V(x,β[i]) , ones ,Infinite2{Float64}(), [B;B] , [3.0;3.0] , pi/6 )
    @test norm(RESULTS[1,3] - exact_values[i]) < 1e-8
end

println("Testing sextic anharmonic oscillators")

# Exact values for the ground state from Chaudhuri1991a and Tater1993, and Flessas1979.
Vsextic(x,c) = c[1]*x.^2 .+ c[2]*x.^4 .+ c[3]*x.^6
c = [1.0 -4.0 1.0 ; 4.0 -6.0 1.0 ; -7.0 0.0 1.0]
exact_values = [-2.0 ; -9.0 ; -2*sqrt(2) ]
enum = [0 , 1 , 0]
ga = 4.0
B = (1/2)^ga / ga
for i=1:size(c)[2]
    (RESULTS, All_Abs_Error_Approx , hoptimal , n, MatrixSizes) = SincEigen( x->Vsextic(x,c[i,:]) , ones ,Infinite2{Float64}(), [B;B] , [ga,ga] , pi/2ga )
    @test norm(RESULTS[enum[i]+1,3] - exact_values[i]) < 1e-8
end

#TODO add examples of SincEigenStop with converged eigenfunctions.

ext: jl


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using AbstractAlgebra
AbstractAlgebra.PrettyPrinting.set_html_as_latex(true)

using AbstractAlgebra
AbstractAlgebra.PrettyPrinting.set_html_as_latex(true)

eqstr(x) = "\\ \\text{is " * string(x) * "}"
eqstr(x::SetElem) = "=" * sprint(x) do io, x; show(io, "text/latex", x) end

function dispeq(lhs, rhs)
    latexstr = "\$\\displaystyle $lhs " * eqstr(rhs) * "\$"
    display("text/html", latexstr)
end

function euclidean(f, g, stop = 1)
    while true
        r = mod(f, g)
        iszero(r) && return g
        degree(r) ≤ stop && return r
        f, g = g, r
    end
end

function rootdeg1(f)
    C = coefficients(f) |> collect
    -first(C)//last(C)
end

function commonroot(f, g)
    m = euclidean(f, g, 1)
    rootdeg1(m)
end

function revpoly(z, x, f)
    m = degree(f)
    C = coefficients(f)
    sum(a*z^(i-1)*x^(m-i+1) for (i, a) in enumerate(C))
end

function val(h::FracElem, θ, L)
    n = L(numerator(h)(θ))
    d = L(denominator(h)(θ))
    n//d
end

function testval(h::FracElem, θ, β, L)
    n = L(numerator(h)(θ))
    d = L(denominator(h)(θ))
    n == β*d
end

function calcallresults(z, x, f, g)
    f_plus = (-1)^degree(f)*f(z - x)
    F_plus = resultant(f_plus, g)
    β_plus = commonroot(g, f_plus)

    f_mult = (-1)^degree(f)*revpoly(z, x, f)
    F_mult = resultant(f_mult, g)
    β_mult = commonroot(g, f_mult)

    K = FractionField(base_ring(z))
    B1, α = K["α"]
    L1 = ResidueField(B1, numerator(f(z))(α))
    B2, β = L1["β"]
    L = ResidueField(B2, numerator(g(z))(β))

    val_F_plus = val(F_plus, α + β, L)
    test_β_plus = testval(β_plus, α + β, β, L)
    val_F_mult = val(F_mult, α * β, L)
    test_β_mult = testval(β_mult, α * β, β, L)

    S_plus = subresultant((-1)^degree(f)*f(z-x), g, 1)
    rootS_plus = rootdeg1(S_plus)

    F_plus, β_plus, F_mult, β_mult, val_F_plus, test_β_plus, val_F_mult, test_β_mult, S_plus, rootS_plus
end

function dispallresults(z, x, f, g)
    F_plus, β_plus, F_mult, β_mult, val_F_plus, test_β_plus, val_F_mult, test_β_mult, S_plus, rootS_plus =
        @time calcallresults(z, x, f, g)
    flush(stdout)
    dispeq("F_\\alpha(x)", f)
    dispeq("F_\\beta(x)", g)
    dispeq("R_{\\alpha + \\beta}(z)", F_plus)
    dispeq("R_{\\alpha + \\beta}(\\alpha + \\beta)", val_F_plus)
    dispeq("\\beta_\\mathrm{plus}(z)", β_plus)
    dispeq("\\beta_\\mathrm{plus}(\\alpha + \\beta) = \\beta", test_β_plus)
    dispeq("R_{\\alpha\\beta}(z)", F_mult)
    dispeq("R_{\\alpha\\beta}(\\alpha\\beta)", val_F_mult)
    dispeq("\\beta_\\mathrm{mult}(z)", β_mult)
    dispeq("\\beta_\\mathrm{mult}(\\alpha\\beta) = \\beta", test_β_mult)
    dispeq("\\text{1-subresultant}", S_plus)
    dispeq("\\text{root of 1-subresultant}", rootS_plus)
end

safecoeff(f, k) = 0 ≤ k ≤ degree(f) ? coeff(f, k) : zero(base_ring(f))

function subresultant_matrix(p::PolyElem{T}, q::PolyElem{T}, k) where T <: RingElement
    check_parent(p, q)
    R = parent(p)
    m = degree(p)
    n = degree(q)
    if length(p) == 0 || length(q) == k || m + n < 2k
       return zero_matrix(R, 0, 0)
    end
    M = zero_matrix(R, m + n - 2k, m + n - 2k)
    x = gen(R)
    for i in 1:n-k
        for j in 1:m+n-2k-1
           M[i, j] = safecoeff(p, m + (i-1) - (j-1))
        end
        M[i, end] = x^(n-k-i)*p
    end
    for i in 1:m-k
        for j in 1:m+n-2k-1
           M[n-k+i, j] = safecoeff(q, n + (i-1) - (j-1))
        end
        M[n-k+i, end] = x^(m-k-i)*q
    end
    return M
end

subresultant(p, q, k) = det(subresultant_matrix(p, q, k))

# %%
R, (a, b, c, p, q, r, s) = ZZ["a", "b", "c", "p", "q", "r", "s"]
K = FractionField(R)
Rz, z = R["z"]
Kz = FractionField(Rz)
Rx, x = Kz["x"]

# %%
f = x^3 + a*x^2 + b*x + c
g = x^4 + p*x^3 + q*x^2 + r*x + s
subresultant_matrix(f, g, 0)

# %%
S0 = sylvester_matrix(f, g)

# %%
subresultant(f, g, 0) == resultant(f, g) == det(S0)

# %% tags=[]
subresultant_matrix(f, g, 1)

# %%
S1 = Rx[
    1 a b c 0
    0 1 a b c*x
    0 0 1 a b*x+c
    1 p q r s*x
    0 1 p q r*x+s
]

# %%
subresultant(f, g, 1)

# %% tags=[]
subresultant(f, g, 1) == det(S1)

# %%
S2 = Rx[
    1 a b*x^2+c*x
    0 1 a*x^2+b*x+c
    1 p q*x^2+r*x+s
]

# %%
subresultant(f, g, 2)

# %%
subresultant(f, g, 2) == det(S2)

# %%
f = x^3 - a
g = x^4 - p
dispallresults(z, x, f, g)

# %%
subresultant_matrix(-f(z-x), g, 1)

# %%
s1 = Rx[
    1 -3z 3z^2 -z^3+a 0
    0 1   -3z  3z^2   (-z^3+a)*x
    0 0   1    -3z    3z^2*x-z^3+a
    1 0   0    0      -p*x
    0 1   0    0      -p
]

# %%
S_plus = subresultant(-f(z-x), g, 1)

# %%
subresultant(-f(z-x), g, 1) == det(s1)

# %%
rootS_plus = rootdeg1(S_plus)
dispeq("\\text{root of 1-subresultant}", rootS_plus)

# %% tags=[]
f = x^3 - a
g = x^3 + p*x + q
dispallresults(z, x, f, g)

# %% tags=[]
f = x^3 - a
g = x^3 + p*x + q
dispallresults(z, x, g, f)

# %%
h1 = g(z-x) + f

# %%
(3z)^2*f - 3z*x*h1

# %%
h2 = (3z)^2*f - (3z*x + 3z^2 + p)*h1

# %%
subresultant(f, g(z-x), 1)

# %%
f = x^2 - a
g = x^2 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^3 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^4 - p
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^5 - p
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^3 - p
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^4 - p
dispallresults(z, x, f, g)

# %%
f = x^3 - a
g = x^5 - p
dispallresults(z, x, f, g)

# %%
f = x^3 + a*x + b
g = x^3 + p*x + q
dispallresults(z, x, f, g)

# %%
f = x^2 - a
g = x^3 + p*x^2 + q*x + r
dispallresults(z, x, f, g)

# %%
f = x^3 + a*x + b
g = x^4 + p*x^2 + q*x + r
dispallresults(z, x, f, g)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using Documenter, Gnuplot

makedocs(sitename="Gnuplot.jl",
         authors = "Giorgio Calderone",
         #format = Documenter.HTML(prettyurls = false),  # uncomment for local use, comment for deployment
         modules=[Gnuplot],
         pages = [
             "Home" => "index.md",
             "Installation" => "install.md",
             "Basic usage" => "basic.md",
             "Advanced usage" => "advanced.md",
             "Package options" => "options.md",
             "Style guide" => "style.md",
             "Gnuplot terminals" => "terminals.md",
             "Plot recipes" => "recipes.md",
             "Examples" => "examples.md",
             "API" => "api.md"
         ])


text:
@doc raw"""
    GradientProblem <: Problem
specify a problem for gradient based algorithms.

# Fields
* `M`            – a manifold $\mathcal M$
* `cost` – a function $F\colon\mathcal M\to\mathbb R$ to minimize
* `gradient`     – the gradient $\nabla F\colon\mathcal M
  \to \mathcal T\mathcal M$ of the cost function $F$

# See also
[`gradient_descent`](@ref)
[`GradientDescentOptions`](@ref)

# """
struct GradientProblem{mT<:Manifold,TCost,TGradient} <: Problem
    M::mT
    cost::TCost
    gradient::TGradient
end
"""
    get_gradient(p,x)

evaluate the gradient of a [`GradientProblem`](@ref)`p` at the point `x`.
"""
function get_gradient(p::GradientProblem, x)
    return p.gradient(x)
end

"""
    DirectionUpdateRule

A general functor, that handles direction update rules. It's field(s) is usually
only a [`StoreOptionsAction`](@ref) by default initialized to the fields required
for the specific coefficient, but can also be replaced by a (common, global)
individual one that provides these values.
"""
abstract type DirectionUpdateRule end

"""
    IdentityUpdateRule <: DirectionUpdateRule

The default gradient direction update is the identity, i.e. it just evaluates the gradient.
"""
struct IdentityUpdateRule <: DirectionUpdateRule end

"""
    AbstractGradientDescentOptions <: Options

A generic [`Options`](@ref) type for gradient based options data.
"""
abstract type AbstractGradientDescentOptions <: Options end

"""
    GradientDescentOptions{P,T} <: AbstractGradientDescentOptions

Describes a Gradient based descent algorithm, with

# Fields
a default value is given in brackets if a parameter can be left out in initialization.

* `x0` – an a point (of type `P`) on a manifold as starting point
* `stopping_criterion` – ([`StopAfterIteration`](@ref)`(100)`) a [`StoppingCriterion`](@ref)
* `stepsize` – ([`ConstantStepsize`](@ref)`(1.)`)a [`Stepsize`](@ref)
* `direction` - ([`IdentityUpdateRule`](@ref)) a processor to compute the gradient
* `retraction_method` – (`ExponentialRetraction()`) the rectraction to use, defaults to
  the exponential map

# Constructor

    GradientDescentOptions(x, stop, s [, retr=ExponentialRetraction()])

construct a Gradient Descent Option with the fields and defaults as above

# See also
[`gradient_descent`](@ref), [`GradientProblem`](@ref)
"""
mutable struct GradientDescentOptions{
    P,TStop<:StoppingCriterion,TStepsize<:Stepsize,TRTM<:AbstractRetractionMethod
} <: AbstractGradientDescentOptions
    x::P
    direction::DirectionUpdateRule
    stop::TStop
    stepsize::TStepsize
    ∇::P
    retraction_method::TRTM
    function GradientDescentOptions{P}(
        initialX::P,
        s::StoppingCriterion=StopAfterIteration(100),
        stepsize::Stepsize=ConstantStepsize(1.0),
        retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
        direction::DirectionUpdateRule=IdentityUpdateRule(),
    ) where {P}
        o = new{P,typeof(s),typeof(stepsize),typeof(retraction_method)}()
        o.x = initialX
        o.stop = s
        o.retraction_method = retraction_method
        o.stepsize = stepsize
        o.direction = direction
        return o
    end
end
function GradientDescentOptions(
    x::P;
    stopping_criterion::StoppingCriterion=StopAfterIteration(100),
    stepsize::Stepsize=ConstantStepsize(1.0),
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    direction::DirectionUpdateRule=IdentityUpdateRule(),
) where {P}
    return GradientDescentOptions{P}(
        x, stopping_criterion, stepsize, retraction_method, direction
    )
end
#
# Processors
#
function (s::IdentityUpdateRule)(p::GradientProblem, o::GradientDescentOptions, i)
    return get_stepsize(p, o, i), get_gradient(p, o.x)
end

"""
    MomentumGradient <: DirectionUpdateRule

Append a momentum to a gradient processor, where the last direction and last iterate are
stored and the new is composed as ``\nabla_i = m*\nabla_{i-1}' - s d_i``,
where ``sd_i`` is the current (inner) direction and ``\nabla_{i-1}'`` is the vector transported
last direction multiplied by momentum ``m``.

# Fields
* `∇` – (`zero_tangent_vector(M,x0)`) the last gradient/direction update added as momentum
* `last_iterate` - remember the last iterate for parallel transporting the last direction
* `momentum` – (`0.2`) factor for momentum
* `direction` – internal [`DirectionUpdateRule`](@ref) to determine directions to
  add the momentum to.
* `vector_transport_method` vector transport method to use

# Constructors
    MomentumGradient(
        p::GradientProlem,
        x0,
        s::DirectionUpdateRule=Gradient();
        ∇=zero_tangent_vector(p.M, o.x), momentum=0.2
       vector_transport_method=ParallelTransport(),
    )

Add momentum to a gradient problem, where by default just a gradient evaluation is used
Equivalently you can also use a `Manifold` `M` instead of the [`GradientProblem`](@ref) `p`.

    MomentumGradient(
        p::StochasticGradientProblem
        x0
        s::DirectionUpdateRule=StochasticGradient();
        ∇=zero_tangent_vector(p.M, x0), momentum=0.2
       vector_transport_method=ParallelTransport(),
    )

Add momentum to a stochastic gradient problem, where by default just a stochastic gradient evaluation is used
"""
mutable struct MomentumGradient{P,T,R<:Real,VTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    ∇::T
    last_iterate::P
    momentum::R
    direction::DirectionUpdateRule
    vector_transport_method::VTM
end
function MomentumGradient(
    p::GradientProblem,
    x0::P,
    s::DirectionUpdateRule=IdentityUpdateRule();
    last_iterate=x0,
    vector_transport_method::VTM=ParallelTransport(),
    ∇=zero_tangent_vector(p.M, x0),
    momentum=0.2,
) where {P,VTM<:AbstractVectorTransportMethod}
    return MomentumGradient{P,typeof(∇),typeof(momentum),VTM}(
        deepcopy(x0), ∇, momentum, s, vector_transport_method
    )
end
function MomentumGradient(
    M::Manifold,
    x0::P,
    s::DirectionUpdateRule=IdentityUpdateRule();
    ∇=zero_tangent_vector(M, x0),
    last_iterate=x0,
    momentum=0.2,
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM<:AbstractVectorTransportMethod}
    return MomentumGradient{P,typeof(∇),typeof(momentum),VTM}(
        deepcopy(x0), ∇, momentum, s, vector_transport_method
    )
end
function (m::MomentumGradient)(p::Problem, o::AbstractGradientDescentOptions, i)
    s, d = m.direction(p, o, i) #get inner direction and step size
    old_d =
        m.momentum *
        vector_transport_to(p.M, m.last_iterate, m.∇, o.x, m.vector_transport_method)
    m.∇ = old_d - s .* d
    m.last_iterate = deepcopy(o.x)
    return s, -m.∇
end

"""
    AverageGradient <: DirectionUpdateRule

Add an average of gradients to a gradient processor. A set of previous directions (from the
inner processor) and the last iterate are stored, average is taken after vector transporting
them to the current iterates tangent space.

# Fields
* `gradients` – (fill(`zero_tangent_vector(M,x0),n)`) the last `n` gradient/direction updates
* `last_iterate` – last iterate (needed to transport the gradients)
* `direction` – internal [`DirectionUpdateRule`](@ref) to determine directions to
  apply the averaging to
* `vector_transport_method` - vector transport method to use

# Constructors
    AverageGradient(
        p::GradientProlem,
        x0,
        n::Int=10
        s::DirectionUpdateRule=IdentityUpdateRule();
        gradients = fill(zero_tangent_vector(p.M, o.x),n),
        last_iterate = deepcopy(x0),
        vector_transport_method = ParallelTransport()
    )

Add average to a gradient problem, `n` determines the size of averaging and `gradients` can be prefilled with some history
Equivalently you can also use a `Manifold` `M` instead of the [`GradientProblem`](@ref) `p`.

    AverageGradient(
        p::StochasticGradientProblem
        x0
        n::Int=10
        s::DirectionUpdateRule=StochasticGradient();
        gradients = fill(zero_tangent_vector(p.M, o.x),n),
        last_iterate = deepcopy(x0),
        vector_transport_method = ParallelTransport()
    )

Add average to a stochastic gradient problem, `n` determines the size of averaging and `gradients` can be prefilled with some history
"""
mutable struct AverageGradient{P,T,VTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    gradients::AbstractVector{T}
    last_iterate::P
    direction::DirectionUpdateRule
    vector_transport_method::VTM
end
function AverageGradient(
    M::Manifold,
    x0::P,
    n::Int=10,
    s::DirectionUpdateRule=IdentityUpdateRule();
    gradients=fill(zero_tangent_vector(M, x0), n),
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM<:AbstractVectorTransportMethod}
    return AverageGradient{P,eltype(gradients),VTM}(
        gradients, x0, s, vector_transport_method
    )
end
function AverageGradient(
    p::GradientProblem,
    x0::P,
    n::Int=10,
    s::DirectionUpdateRule=IdentityUpdateRule();
    gradients=fill(zero_tangent_vector(p.M, x0), n),
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM}
    return AverageGradient{P,eltype(gradients),VTM}(
        gradients, deepcopy(x0), s, vector_transport_method
    )
end
function (a::AverageGradient)(p::Problem, o::AbstractGradientDescentOptions, i)
    pop!(a.gradients)
    s, d = a.direction(p, o, i) #get inner gradient and step
    a.gradients = vcat([d], a.gradients)
    for i in 1:(length(a.gradients) - 1) #transport & shift inplace
        vector_transport_to!(
            p.M,
            a.gradients[i],
            a.last_iterate,
            a.gradients[i + 1],
            o.x,
            a.vector_transport_method,
        )
    end
    a.gradients[1] = d
    a.last_iterate = deepcopy(o.x)
    return s, -1 / length(a.gradients) .* sum(a.gradients)
end

@doc raw"""
    Nesterov <: DirectionUpdateRule

## Fields
* `γ`
* `μ` the strong convexity coefficient
* `v` (=``=v_k``, ``v_0=x_0``) an interims point to compute the next gradient evaluation point `y_k`
* `shrinkage` (`= i -> 0.8`) a function to compute the shrinkage ``β_k`` per iterate.

Let's assume ``f`` is ``L``-Lipschitz and ``μ``-strongly convex.
Given
* a step size ``h_k<\frac{1}{L}`` (from the [`GradientDescentOptions`](@ref)
* a `shrinkage` parameter ``β_k``
* and a current iterate $x_k$
* as well as the interims values $γ_k`` and ``v_k`` from the previous iterate.

This compute a Nesterov type update using the following steps, see [^ZhangSra2018]

1. Copute the positive root, i.e. ``α_k∈(0,1)`` of ``α^2 = h_k\bigl((1-α_k)γ_k+α_k μ\bigr)``.
2. Set ``\bar γ_k+1 = (1-α_k)γ_k + α_kμ``
3. ``y_k = \operatorname{retr}_{x_k}\Bigl(\frac{α_kγ_k}{γ_k + α_kμ}\operatorname{retr}^{-1}_{x_k}v_k \Bigr)``
4. ``x_{k+1} = \operatorname{retr}_{y_k}(-h_k ∇f(y_k))``
5. ``v_{k+1} = `\operatorname{retr}_{y_k}\Bigl(\frac{(1-α_k)γ_k}{\barγ_k}\operatorname{retr}_{y_k}^{-1}(v_k) - \frac{α_k}{\bar γ_{k+1}}∇f(y_k) \Bigr)``
6. ``γ_{k+1} = \frac{1}{1+β_k}\bar γ_{k+1}``

Then the direction from ``x_k`` to ``x_k+1``, i.e. ``d = \operatorname{retr}^{-1}_{x_k}x_{k+1}`` is returned.

# Constructor
    Nesterov(x0::P, γ=0.001, μ=0.9, schrinkage = k -> 0.8;
        inverse_retraction_method=LogarithmicInverseRetraction())

Initialize the Nesterov acceleration, where `x0` initializes `v`.

[^ZhangSra2018]:
    > H. Zhang, S. Sra: _Towards Riemannian Accelerated Gradient Methods_,
    > Preprint, 2018, arXiv: [1806.02812](https://arxiv.org/abs/1806.02812)
"""
mutable struct Nesterov{P,T<:Real} <: DirectionUpdateRule
    γ::T
    μ::T
    v::P
    shrinkage::Function
    inverse_retraction_method::AbstractInverseRetractionMethod
end
function Nesterov(
    x0::P,
    γ::T=0.001,
    μ::T=0.9,
    shrinkage::Function=i -> 0.8;
    inverse_retraction_method::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
) where {P,T}
    return Nesterov{P,T}(γ, μ, deepcopy(x0), shrinkage, inverse_retraction_method)
end
function (s::Nesterov)(p::GradientProblem, o::AbstractGradientDescentOptions, i)
    h = get_stepsize(p, o, i)
    α = (h * (s.γ - s.μ) + sqrt(h^2 * (s.γ - s.μ)^2 + 4 * h * s.γ)) / 2
    γbar = (1 - α) * s.γ + α * s.μ
    y = retract(p.M, o.x, (α * s.γ) / (s.γ + α * s.μ) .* inverse_retract(p.M, o.x, s.v))
    gradf_yk = get_gradient(p, y)
    xn = retract(p.M, y, -h * gradf_yk)
    d =
        ((1 - α) * s.γ) / γbar .*
        inverse_retract(p.M, y, s.v, s.inverse_retraction_method) - α / γbar .* gradf_yk
    s.v = retract(p.M, y, d, o.retraction_method)
    s.γ = 1 / (1 + s.shrinkage(i)) * γbar
    return h, -1 / h .* inverse_retract(p.M, o.x, xn) # outer update
end

#
# Conjugate Gradient Descent
#

@doc raw"""
    ConjugateGradientOptions <: Options

specify options for a conjugate gradient descent algoritm, that solves a
[`GradientProblem`].

# Fields
* `x` – the current iterate, a point on a manifold
* `∇` – the current gradient
* `δ` – the current descent direction, i.e. also tangent vector
* `β` – the current update coefficient rule, see .
* `coefficient` – a [`DirectionUpdateRule`](@ref) function to determine the new `β`
* `stepsize` – a [`Stepsize`](@ref) function
* `stop` – a [`StoppingCriterion`](@ref)
* `retraction_method` – (`ExponentialRetraction()`) a type of retraction

# See also
[`conjugate_gradient_descent`](@ref), [`GradientProblem`](@ref), [`ArmijoLinesearch`](@ref)
"""
mutable struct ConjugateGradientDescentOptions{
    T,
    TCoeff<:DirectionUpdateRule,
    TStepsize<:Stepsize,
    TStop<:StoppingCriterion,
    TRetr<:AbstractRetractionMethod,
    TVTM<:AbstractVectorTransportMethod,
} <: Options
    x::T
    ∇::T
    δ::T
    β::Float64
    coefficient::TCoeff
    stepsize::TStepsize
    stop::TStop
    retraction_method::TRetr
    vector_transport_method::TVTM
    function ConjugateGradientDescentOptions{T}(
        x0::T,
        sC::StoppingCriterion,
        s::Stepsize,
        dC::DirectionUpdateRule,
        retr::AbstractRetractionMethod=ExponentialRetraction(),
        vtr::AbstractVectorTransportMethod=ParallelTransport(),
    ) where {T}
        o = new{T,typeof(dC),typeof(s),typeof(sC),typeof(retr),typeof(vtr)}()
        o.x = x0
        o.stop = sC
        o.retraction_method = retr
        o.stepsize = s
        o.coefficient = dC
        o.vector_transport_method = vtr
        return o
    end
end
function ConjugateGradientDescentOptions(
    x::T,
    sC::StoppingCriterion,
    s::Stepsize,
    dU::DirectionUpdateRule,
    retr::AbstractRetractionMethod=ExponentialRetraction(),
    vtr::AbstractVectorTransportMethod=ParallelTransport(),
) where {T}
    return ConjugateGradientDescentOptions{T}(x, sC, s, dU, retr, vtr)
end

@doc raw"""
    ConjugateDescentCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^Flethcer1987] adapted to manifolds:

$\beta_k =
\frac{ \lVert \xi_{k+1} \rVert_{x_{k+1}}^2 }
{\langle -\delta_k,\xi_k \rangle_{x_k}}.$

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    ConjugateDescentCoefficient(a::StoreOptionsAction=())

Construct the conjugate descnt coefficient update rule, a new storage is created by default.

[^Flethcer1987]:
    > R. Fletcher, __Practical Methods of Optimization vol. 1: Unconstrained Optimization__
    > John Wiley & Sons, New York, 1987. doi [10.1137/1024028](https://doi.org/10.1137/1024028)
"""
mutable struct ConjugateDescentCoefficient <: DirectionUpdateRule
    storage::StoreOptionsAction
    function ConjugateDescentCoefficient(a::StoreOptionsAction=StoreOptionsAction((:x, :∇)))
        return new(a)
    end
end
function (u::ConjugateDescentCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇]))
        update_storage!(u.storage, o) # if not given store current as old
        return 0.0
    end
    xOld, ∇Old = get_storage.(Ref(u.storage), [:x, :∇])
    update_storage!(u.storage, o)
    return inner(p.M, o.x, o.∇, o.∇) / inner(p.M, xOld, -o.δ, ∇Old)
end

@doc raw"""
    DaiYuanCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^DaiYuan1999]

adapted to manifolds: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$,
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

Then the coefficient reads

````math
\beta_k =
\frac{ \lVert \xi_{k+1} \rVert_{x_{k+1}}^2 }
{\langle P_{x_{k+1}\gets x_k}\delta_k, \nu_k \rangle_{x_{k+1}}}.
````

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    DaiYuanCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=(),
    )

Construct the Dai Yuan coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

[^DaiYuan1999]:
    > [Y. H. Dai and Y. Yuan, A nonlinear conjugate gradient method with a strong global convergence property,
    > SIAM J. Optim., 10 (1999), pp. 177–182.
    > doi: [10.1137/S1052623497318992](https://doi.org/10.1137/S1052623497318992)
"""
mutable struct DaiYuanCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function DaiYuanCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=StoreOptionsAction((:x, :∇, :δ)),
    )
        return new{typeof(t)}(t, a)
    end
end
function (u::DaiYuanCoefficient)(p::GradientProblem, o::ConjugateGradientDescentOptions, i)
    if !all(has_storage.(Ref(u.storage), [:x, :∇, :δ]))
        update_storage!(u.storage, o) # if not given store current as old
        return 0.0
    end
    xOld, ∇Old, δOld = get_storage.(Ref(u.storage), [:x, :∇, :δ])
    update_storage!(u.storage, o)

    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    ν = o.∇ - ∇tr #notation y from [HZ06]
    δtr = vector_transport_to(p.M, xOld, δOld, o.x, u.transport_method)
    return inner(p.M, o.x, o.∇, o.∇) / inner(p.M, xOld, δtr, ν)
end

@doc raw"""
    FletcherReevesCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^FletcherReeves1964] adapted to manifolds:

````math
\beta_k =
\frac{\lVert \xi_{k+1}\rVert_{x_{k+1}}^2}{\lVert \xi_{k}\rVert_{x_{k}}^2}.
````

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    FletcherReevesCoefficient(a::StoreOptionsAction=())

Construct the Fletcher Reeves coefficient update rule, a new storage is created by default.

[^FletcherReeves1964]:
    > R. Fletcher and C. Reeves, Function minimization by conjugate gradients,
    > Comput. J., 7 (1964), pp. 149–154.
    > doi: [10.1093/comjnl/7.2.149](http://dx.doi.org/10.1093/comjnl/7.2.149)
"""
mutable struct FletcherReevesCoefficient <: DirectionUpdateRule
    storage::StoreOptionsAction
    FletcherReevesCoefficient(a::StoreOptionsAction=StoreOptionsAction((:x, :∇))) = new(a)
end
function (u::FletcherReevesCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇]))
        update_storage!(u.storage, o) # if not given store current as old
    end
    xOld, ∇Old = get_storage.(Ref(u.storage), [:x, :∇])
    update_storage!(u.storage, o)
    return inner(p.M, o.x, o.∇, o.∇) / inner(p.M, xOld, ∇Old, ∇Old)
end

@doc raw"""
    HagerZhangCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the variables with
prequel `Old` based on [^HagerZhang2005]
adapted to manifolds: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$,
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

````math
\beta_k = \Bigl\langle\nu_k -
\frac{ 2\lVert \nu_k\rVert_{x_{k+1}}^2 }{ \langle P_{x_{k+1}\gets x_k}\delta_k, \nu_k \rangle_{x_{k+1}} }
P_{x_{k+1}\gets x_k}\delta_k,
\frac{\xi_{k+1}}{ \langle P_{x_{k+1}\gets x_k}\delta_k, \nu_k \rangle_{x_{k+1}} }
\Bigr\rangle_{x_{k+1}}.
````

This method includes a numerical stability proposed by those authors.

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    HagerZhangCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=(),
    )

Construct the Hager Zhang coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

[^HagerZhang2005]:
    > [W. W. Hager and H. Zhang, __A new conjugate gradient method with guaranteed descent and an efficient line search__,
    > SIAM J. Optim, (16), pp. 170-192, 2005.
    > doi: [10.1137/030601880](https://doi.org/10.1137/030601880)
"""
mutable struct HagerZhangCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function HagerZhangCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=StoreOptionsAction((:x, :∇, :δ)),
    )
        return new{typeof(t)}(t, a)
    end
end
function (u::HagerZhangCoefficient)(
    p::P, o::O, i
) where {P<:GradientProblem,O<:ConjugateGradientDescentOptions}
    if !all(has_storage.(Ref(u.storage), [:x, :∇, :δ]))
        update_storage!(u.storage, o) # if not given store current as old
        return 0.0
    end
    xOld, ∇Old, δOld = get_storage.(Ref(u.storage), [:x, :∇, :δ])
    update_storage!(u.storage, o)

    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    ν = o.∇ - ∇tr #notation y from [HZ06]
    δtr = vector_transport_to(p.M, xOld, δOld, o.x, u.transport_method)
    denom = inner(p.M, o.x, δtr, ν)
    νknormsq = inner(p.M, o.x, ν, ν)
    β = inner(p.M, o.x, ν, o.∇) / denom - 2 * νknormsq * inner(p.M, o.x, δtr, o.∇) / denom^2
    # Numerical stability from Manopt / Hager-Zhang paper
    ξn = norm(p.M, o.x, o.∇)
    η = -1 / (ξn * min(0.01, norm(p.M, xOld, ∇Old)))
    return max(β, η)
end

@doc raw"""
    HeestenesStiefelCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^HeestensStiefel1952]

adapted to manifolds as follows: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$.
Then the update reads

````math
\beta_k = \frac{\langle \xi_{k+1}, \nu_k \rangle_{x_{k+1}} }
    { \langle P_{x_{k+1}\gets x_k} \delta_k, \nu_k\rangle_{x_{k+1}} },
````
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

# Constructor
    HeestenesStiefelCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=()
    )

Construct the Heestens Stiefel coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

See also [`conjugate_gradient_descent`](@ref)

[^HeestensStiefel1952]:
    > M.R. Hestenes, E.L. Stiefel, Methods of conjugate gradients for solving linear systems,
    > J. Research Nat. Bur. Standards, 49 (1952), pp. 409–436.
    > doi: [10.6028/jres.049.044](http://dx.doi.org/10.6028/jres.049.044)
"""
mutable struct HeestenesStiefelCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function HeestenesStiefelCoefficient(
        transport_method::AbstractVectorTransportMethod=ParallelTransport(),
        storage_action::StoreOptionsAction=StoreOptionsAction((:x, :∇, :δ)),
    )
        return new{typeof(transport_method)}(transport_method, storage_action)
    end
end
function (u::HeestenesStiefelCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇, :δ]))
        update_storage!(u.storage, o) # if not given store current as old
        return 0.0
    end
    xOld, ∇Old, δOld = get_storage.(Ref(u.storage), [:x, :∇, :δ])
    update_storage!(u.storage, o)
    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    δtr = vector_transport_to(p.M, xOld, δOld, o.x, u.transport_method)
    ν = o.∇ - ∇tr #notation from [HZ06]
    β = inner(p.M, o.x, o.∇, ν) / inner(p.M, o.x, δtr, ν)
    return max(0, β)
end

@doc raw"""
    LiuStoreyCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^LuiStorey1991]
adapted to manifolds: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$,
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

Then the coefficient reads

````math
\beta_k = -
\frac{ \langle \xi_{k+1},\nu_k \rangle_{x_{k+1}} }
{\langle \delta_k,\xi_k \rangle_{x_k}}.
````

See also [`conjugate_gradient_descent`](@ref)

# Constructor
    LiuStoreyCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=()
    )

Construct the Lui Storey coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

[^LuiStorey1991]:
    > [Y. Liu and C. Storey, Efficient generalized conjugate gradient algorithms, Part 1: Theory
    > J. Optim. Theory Appl., 69 (1991), pp. 129–137.
    > doi: [10.1007/BF00940464](https://doi.org/10.1007/BF00940464)
"""
mutable struct LiuStoreyCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function LiuStoreyCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=StoreOptionsAction((:x, :∇, :δ)),
    )
        return new{typeof(t)}(t, a)
    end
end
function (u::LiuStoreyCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇, :δ]))
        update_storage!(u.storage, o) # if not given store current as old
    end
    xOld, ∇Old, δOld = get_storage.(Ref(u.storage), [:x, :∇, :δ])
    update_storage!(u.storage, o)
    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    ν = o.∇ - ∇tr # notation y from [HZ06]
    return inner(p.M, o.x, o.∇, ν) / inner(p.M, xOld, -δOld, ∇Old)
end

@doc raw"""
    PolakRibiereCoefficient <: DirectionUpdateRule

Computes an update coefficient for the conjugate gradient method, where
the [`ConjugateGradientDescentOptions`](@ref)` o` include the last iterates
$x_k,\xi_k$, the current iterates $x_{k+1},\xi_{k+1}$ and the last update
direction $\delta=\delta_k$, where the last three ones are stored in the
variables with prequel `Old` based on [^PolakRibiere1969][^Polyak1969]

adapted to manifolds: let $\nu_k = \xi_{k+1} - P_{x_{k+1}\gets x_k}\xi_k$,
where $P_{a\gets b}(\cdot)$ denotes a vector transport from the tangent space at $a$ to $b$.

Then the update reads
````math
\beta_k =
\frac{ \langle \xi_{k+1}, \nu_k \rangle_{x_{k+1}} }
{\lVert \xi_k \rVert_{x_k}^2 }.
````

# Constructor
    PolakRibiereCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=()
    )

Construct the PolakRibiere coefficient update rule, where the parallel transport is the
default vector transport and a new storage is created by default.

See also [`conjugate_gradient_descent`](@ref)

[^PolakRibiere1969]:
    > E. Polak, G. Ribiere, Note sur la convergence de méthodes de directions conjuguées
    > ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, Tome 3 (1969) no. R1, p. 35-43,
    > url: [http://www.numdam.org/item/?id=M2AN_1969__3_1_35_0](http://www.numdam.org/item/?id=M2AN_1969__3_1_35_0)

[^Polyak1969]:
    > B. T. Polyak, The conjugate gradient method in extreme problems,
    > USSR Comp. Math. Math. Phys., 9 (1969), pp. 94–112.
    > doi: [10.1016/0041-5553(69)90035-4](https://doi.org/10.1016/0041-5553(69)90035-4)
"""
mutable struct PolakRibiereCoefficient{TVTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    transport_method::TVTM
    storage::StoreOptionsAction
    function PolakRibiereCoefficient(
        t::AbstractVectorTransportMethod=ParallelTransport(),
        a::StoreOptionsAction=StoreOptionsAction((:x, :∇)),
    )
        return new{typeof(t)}(t, a)
    end
end
function (u::PolakRibiereCoefficient)(
    p::GradientProblem, o::ConjugateGradientDescentOptions, i
)
    if !all(has_storage.(Ref(u.storage), [:x, :∇]))
        update_storage!(u.storage, o) # if not given store current as old
    end
    xOld, ∇Old = get_storage.(Ref(u.storage), [:x, :∇])
    update_storage!(u.storage, o)

    ∇tr = vector_transport_to(p.M, xOld, ∇Old, o.x, u.transport_method)
    ν = o.∇ - ∇tr
    β = inner(p.M, o.x, o.∇, ν) / inner(p.M, xOld, ∇Old, ∇Old)
    return max(0, β)
end

@doc raw"""
    SteepestDirectionUpdateRule <: DirectionUpdateRule

The simplest rule to update is to have no influence of the last direction and
hence return an update $\beta = 0$ for all [`ConjugateGradientDescentOptions`](@ref)` o`

See also [`conjugate_gradient_descent`](@ref)
"""
struct SteepestDirectionUpdateRule <: DirectionUpdateRule end
function (u::SteepestDirectionUpdateRule)(
    ::GradientProblem, ::ConjugateGradientDescentOptions, i
)
    return 0.0
end

@doc raw"""
    DebugGradient <: DebugAction

debug for the gradient evaluated at the current iterate

# Constructors
    DebugGradient([long=false,p=print])

display the short (`false`) or long (`true`) default text for the gradient.

    DebugGradient(prefix[, p=print])

display the a `prefix` in front of the gradient.
"""
mutable struct DebugGradient <: DebugAction
    io::IO
    prefix::String
    function DebugGradient(long::Bool=false, io::IO=stdout)
        return new(io, long ? "Gradient: " : "∇F(x):")
    end
    DebugGradient(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugGradient)(::GradientProblem, o::GradientDescentOptions, i::Int)
    print(d.io, (i >= 0) ? d.prefix * "" * string(o.∇) : "")
    return nothing
end

@doc raw"""
    DebugGradientNorm <: DebugAction

debug for gradient evaluated at the current iterate.

# Constructors
    DebugGradientNorm([long=false,p=print])

display the short (`false`) or long (`true`) default text for the gradient norm.

    DebugGradientNorm(prefix[, p=print])

display the a `prefix` in front of the gradientnorm.
"""
mutable struct DebugGradientNorm <: DebugAction
    io::IO
    prefix::String
    function DebugGradientNorm(long::Bool=false, io::IO=stdout)
        return new(io, long ? "Norm of the Gradient: " : "|∇F(x)|:")
    end
    DebugGradientNorm(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugGradientNorm)(p::GradientProblem, o::Options, i::Int)
    print(d.io, (i >= 0) ? d.prefix * "$(norm(p.M,o.x,o.∇))" : "")
    return nothing
end

@doc raw"""
    DebugStepsize <: DebugAction

debug for the current step size.

# Constructors
    DebugStepsize([long=false,p=print])

display the short (`false`) or long (`true`) default text for the step size.

    DebugStepsize(prefix[, p=print])

display the a `prefix` in front of the step size.
"""
mutable struct DebugStepsize <: DebugAction
    io::IO
    prefix::String
    function DebugStepsize(long::Bool=false, io::IO=stdout)
        return new(io, long ? "step size:" : "s:")
    end
    DebugStepsize(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugStepsize)(p::GradientProblem, o::GradientDescentOptions, i::Int)
    print(d.io, (i > 0) ? d.prefix * "$(get_last_stepsize(p,o,i))" : "")
    return nothing
end

#
# Records
#
@doc raw"""
    RecordGradient <: RecordAction

record the gradient evaluated at the current iterate

# Constructors
    RecordGradient(ξ)

initialize the [`RecordAction`](@ref) to the corresponding type of the tangent vector.
"""
mutable struct RecordGradient{T} <: RecordAction
    recorded_values::Array{T,1}
    RecordGradient{T}() where {T} = new(Array{T,1}())
end
RecordGradient(ξ::T) where {T} = RecordGradient{T}()
function (r::RecordGradient{T})(
    ::GradientProblem, o::GradientDescentOptions, i::Int
) where {T}
    return record_or_reset!(r, o.∇, i)
end

@doc raw"""
    RecordGradientNorm <: RecordAction

record the norm of the current gradient
"""
mutable struct RecordGradientNorm <: RecordAction
    recorded_values::Array{Float64,1}
    RecordGradientNorm() = new(Array{Float64,1}())
end
function (r::RecordGradientNorm)(
    p::P, o::O, i::Int
) where {P<:GradientProblem,O<:GradientDescentOptions}
    return record_or_reset!(r, norm(p.M, o.x, o.∇), i)
end

@doc raw"""
    RecordStepsize <: RecordAction

record the step size
"""
mutable struct RecordStepsize <: RecordAction
    recorded_values::Array{Float64,1}
    RecordStepsize() = new(Array{Float64,1}())
end
function (r::RecordStepsize)(
    p::P, o::O, i::Int
) where {P<:GradientProblem,O<:GradientDescentOptions}
    return record_or_reset!(r, get_last_stepsize(p, o, i), i)
end

@doc raw"""
    AbstractQuasiNewtonDirectionUpdate

An abstract represresenation of an Quasi Newton Update rule to determine the next direction
given current [`QuasiNewtonOptions`](@ref).

All subtypes should be functors, i.e. one should be able to call them as `H(M,x,d)` to compute a new direction update.
"""
abstract type AbstractQuasiNewtonDirectionUpdate end

@doc raw"""
    AbstractQuasiNewtonUpdateRule

Specify a type for the different [`AbstractQuasiNewtonDirectionUpdate`](@ref)s.
"""
abstract type AbstractQuasiNewtonUpdateRule end

@doc raw"""
    BFGS <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian BFGS update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{H}_k^\mathrm{BFGS}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
H^\mathrm{BFGS}_{k+1} = \widetilde{H}^\mathrm{BFGS}_k  + \frac{y_k y^{\mathrm{T}}_k }{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{H}^\mathrm{BFGS}_k s_k s^{\mathrm{T}}_k \widetilde{H}^\mathrm{BFGS}_k }{s^{\mathrm{T}}_k \widetilde{H}^\mathrm{BFGS}_k s_k}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
∇f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(∇ f(x_k)) \in T_{x_{k+1}} \mathcal{M},
```

respectively.
"""
struct BFGS <: AbstractQuasiNewtonUpdateRule end
@doc raw"""
    InverseBFGS <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the inverse Riemanian BFGS update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{B}_k^\mathrm{BFGS}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
B^\mathrm{BFGS}_{k+1}  = \Bigl(
  \mathrm{id}_{T_{x_{k+1}} \mathcal{M}} - \frac{s_k y^{\mathrm{T}}_k }{s^{\mathrm{T}}_k y_k}
\Bigr)
\widetilde{B}^\mathrm{BFGS}_k
\Bigl(
  \mathrm{id}_{T_{x_{k+1}} \mathcal{M}} - \frac{y_k s^{\mathrm{T}}_k }{s^{\mathrm{T}}_k y_k}
\Bigr) + \frac{s_k s^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
∇f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(∇ f(x_k)) \in T_{x_{k+1}} \mathcal{M},
```

respectively.

"""
struct InverseBFGS <: AbstractQuasiNewtonUpdateRule end
@doc raw"""
    DFP <: AbstractQuasiNewtonUpdateRule

indicates in an [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian DFP update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{H}_k^\mathrm{DFP}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
H^\mathrm{DFP}_{k+1} = \Bigl(
  \mathrm{id}_{T_{x_{k+1}} \mathcal{M}} - \frac{y_k s^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
\Bigr)
\widetilde{H}^\mathrm{DFP}_k
\Bigl(
  \mathrm{id}_{T_{x_{k+1}} \mathcal{M}} - \frac{s_k y^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
\Bigr) + \frac{y_k y^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
∇f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(∇ f(x_k)) \in T_{x_{k+1}} \mathcal{M},
```

respectively.

"""
struct DFP <: AbstractQuasiNewtonUpdateRule end
@doc raw"""
    InverseDFP <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the inverse Riemanian DFP update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{B}_k^\mathrm{DFP}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
B^\mathrm{DFP}_{k+1} = \widetilde{B}^\mathrm{DFP}_k
+ \frac{s_k s^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
- \frac{\widetilde{B}^\mathrm{DFP}_k y_k y^{\mathrm{T}}_k \widetilde{B}^\mathrm{DFP}_k}{y^{\mathrm{T}}_k \widetilde{B}^\mathrm{DFP}_k y_k}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
∇f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(∇ f(x_k)) \in T_{x_{k+1}} \mathcal{M},
```

respectively.

"""
struct InverseDFP <: AbstractQuasiNewtonUpdateRule end
@doc raw"""
    SR1 <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian SR1 update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{H}_k^\mathrm{SR1}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
H^\mathrm{SR1}_{k+1} = \widetilde{H}^\mathrm{SR1}_k
+ \frac{
  (y_k - \widetilde{H}^\mathrm{SR1}_k s_k) (y_k - \widetilde{H}^\mathrm{SR1}_k s_k)^{\mathrm{T}}
}{
(y_k - \widetilde{H}^\mathrm{SR1}_k s_k)^{\mathrm{T}} s_k
}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
∇f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(∇ f(x_k)) \in T_{x_{k+1}} \mathcal{M},
```

respectively.

This method can be stabilized by only performing the update if denominator is larger than
``r\lVert s_k\rVert_{x_{k+1}}\lVert y_k - \widetilde{H}^\mathrm{SR1}_k s_k \rVert_{x_{k+1}}``
for some ``r>0``. For more details, see Section 6.2 in [^NocedalWright2006]

[^NocedalWright2006]:
    > Nocedal, J., Wright, S.: Numerical Optimization, Second Edition, Springer, 2006.
    > doi: [10.1007/978-0-387-40065-5](https://doi.org/10.1007/978-0-387-40065-5)

# Constructor
    SR1(r::Float64=-1.0)

Generate the `SR1` update, which by default does not include the check (since the default sets ``t<0```)

"""
struct SR1 <: AbstractQuasiNewtonUpdateRule
    r::Float64
    SR1(r::Float64=-1.0) = new(r)
end

@doc raw"""
    InverseSR1 <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the inverse Riemanian SR1 update is used in the Riemannian quasi-Newton method.

We denote by ``\widetilde{B}_k^\mathrm{SR1}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
B^\mathrm{SR1}_{k+1} = \widetilde{B}^\mathrm{SR1}_k
+ \frac{
  (s_k - \widetilde{B}^\mathrm{SR1}_k y_k) (s_k - \widetilde{B}^\mathrm{SR1}_k y_k)^{\mathrm{T}}
}{
  (s_k - \widetilde{B}^\mathrm{SR1}_k y_k)^{\mathrm{T}} y_k
}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
∇f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(∇ f(x_k)) \in T_{x_{k+1}} \mathcal{M},
```

respectively.

This method can be stabilized by only performing the update if denominator is larger than
``r\lVert y_k\rVert_{x_{k+1}}\lVert s_k - \widetilde{H}^\mathrm{SR1}_k y_k \rVert_{x_{k+1}}``
for some ``r>0``. For more details, see Section 6.2 in [^NocedalWright2006].

# Constructor
    InverseSR1(r::Float64=-1.0)

Generate the `InverseSR1` update, which by default does not include the check,
since the default sets ``t<0```.
"""
struct InverseSR1 <: AbstractQuasiNewtonUpdateRule
    r::Float64
    InverseSR1(r::Float64=-1.0) = new(r)
end

@doc raw"""
    Broyden <: AbstractQuasiNewtonUpdateRule

indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian Broyden update is used in the Riemannian quasi-Newton method, which is as a convex combination of [`BFGS`](@ref) and [`DFP`](@ref).

We denote by ``\widetilde{H}_k^\mathrm{Br}`` the operator concatenated with a vector transport and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
H^\mathrm{Br}_{k+1} = \widetilde{H}^\mathrm{Br}_k
  - \frac{\widetilde{H}^\mathrm{Br}_k s_k s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k}{s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k s_k} + \frac{y_k y^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
  + φ_k s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k s_k
  \Bigl(
        \frac{y_k}{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{H}^\mathrm{Br}_k s_k}{s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k s_k}
  \Bigr)
  \Bigl(
        \frac{y_k}{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{H}^\mathrm{Br}_k s_k}{s^{\mathrm{T}}_k \widetilde{H}^\mathrm{Br}_k s_k}
  \Bigr)^{\mathrm{T}}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
∇f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(∇ f(x_k)) \in T_{x_{k+1}} \mathcal{M},
```

respectively, and ``φ_k`` is the Broyden factor which is `:constant` by default but can also be set to `:Davidon`.

# Constructor
    Broyden(φ, update_rule::Symbol = :constant)

"""
mutable struct Broyden <: AbstractQuasiNewtonUpdateRule
    φ::Float64
    update_rule::Symbol
end
Broyden(φ::Float64) = Broyden(φ, :constant)

@doc raw"""
    InverseBroyden <: AbstractQuasiNewtonUpdateRule

Indicates in [`AbstractQuasiNewtonDirectionUpdate`](@ref) that the Riemanian Broyden update
is used in the Riemannian quasi-Newton method, which is as a convex combination
of [`InverseBFGS`](@ref) and [`InverseDFP`](@ref).

We denote by ``\widetilde{H}_k^\mathrm{Br}`` the operator concatenated with a vector transport
and its inverse before and after to act on ``x_{k+1} = R_{x_k}(α_k η_k)``.
Then the update formula reads

```math
B^\mathrm{Br}_{k+1} = \widetilde{B}^\mathrm{Br}_k
 - \frac{\widetilde{B}^\mathrm{Br}_k y_k y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k}{y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k y_k}
   + \frac{s_k s^{\mathrm{T}}_k}{s^{\mathrm{T}}_k y_k}
 + φ_k y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k y_k
 \Bigl(
     \frac{s_k}{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{B}^\mathrm{Br}_k y_k}{y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k y_k}
    \Bigr) \Bigl(
        \frac{s_k}{s^{\mathrm{T}}_k y_k} - \frac{\widetilde{B}^\mathrm{Br}_k y_k}{y^{\mathrm{T}}_k \widetilde{B}^\mathrm{Br}_k y_k}
 \Bigr)^{\mathrm{T}}
```

where ``s_k`` and ``y_k`` are the coordinate vectors with respect to the current basis (from [`QuasiNewtonOptions`](@ref)) of

```math
T^{S}_{x_k, α_k η_k}(α_k η_k) \quad\text{and}\quad
∇f(x_{k+1}) - T^{S}_{x_k, α_k η_k}(∇ f(x_k)) \in T_{x_{k+1}} \mathcal{M},
```

respectively, and ``φ_k`` is the Broyden factor which is `:constant` by default but can also be set to `:Davidon`.

# Constructor
    InverseBroyden(φ, update_rule::Symbol = :constant)

"""
mutable struct InverseBroyden <: AbstractQuasiNewtonUpdateRule
    φ::Float64
    update_rule::Symbol
end
InverseBroyden(φ::Float64) = InverseBroyden(φ, :constant)

@doc raw"""
    QuasiNewtonOptions <: Options

These Quasi Newton [`Options`](@ref) represent any quasi-Newton based method and can be
used with any update rule for the direction.

# Fields
* `x` – the current iterate, a point on a manifold
* `∇` – the current gradient
* `sk` – the current step
* `yk` the current gradient difference
* `direction_update` - an [`AbstractQuasiNewtonDirectionUpdate`](@ref) rule.
* `retraction_method` – an `AbstractRetractionMethod`
* `stop` – a [`StoppingCriterion`](@ref)

# See also
[`GradientProblem`](@ref)
"""
mutable struct QuasiNewtonOptions{
    P,
    T,
    U<:AbstractQuasiNewtonDirectionUpdate,
    SC<:StoppingCriterion,
    S<:Stepsize,
    RTR<:AbstractRetractionMethod,
    VT<:AbstractVectorTransportMethod,
} <: Options
    x::P
    ∇::T
    sk::T
    yk::T
    direction_update::U
    retraction_method::RTR
    stepsize::S
    stop::SC
    vector_transport_method::VT
end
function QuasiNewtonOptions(
    x::P,
    ∇::T,
    direction_update::U,
    stop::SC,
    stepsize::S;
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    vector_transport_method::AbstractVectorTransportMethod=ParallelTransport(),
) where {P,T,U<:AbstractQuasiNewtonDirectionUpdate,SC<:StoppingCriterion,S<:Stepsize}
    return QuasiNewtonOptions{
        P,T,U,SC,S,typeof(retraction_method),typeof(vector_transport_method)
    }(
        x,
        ∇,
        deepcopy(∇),
        deepcopy(∇),
        direction_update,
        retraction_method,
        stepsize,
        stop,
        vector_transport_method,
    )
end

@doc raw"""
    QuasiNewtonMatrixDirectionUpdate <: AbstractQuasiNewtonDirectionUpdate

These [`AbstractQuasiNewtonDirectionUpdate`](@ref)s represent any quasi-Newton update rule, where the operator is stored as a matrix. A distinction is made between the update of the approximation of the Hessian, ``H_k \mapsto H_{k+1}``, and the update of the approximation of the Hessian inverse, ``B_k \mapsto B_{k+1}``. For the first case, the coordinates of the search direction ``η_k`` with respect to a basis ``\{b_i\}^{n}_{i=1}`` are determined by solving a linear system of equations, i.e.

```math
\text{Solve} \quad \hat{η_k} = - H_k \widehat{∇f(x_k)}
```

where ``H_k`` is the matrix representing the operator with respect to the basis ``\{b_i\}^{n}_{i=1}`` and ``\widehat{∇f(x_k)}`` represents the coordinates of the gradient of the objective function ``f`` in ``x_k`` with respect to the basis ``\{b_i\}^{n}_{i=1}``.
If a method is chosen where Hessian inverse is approximated, the coordinates of the search direction ``η_k`` with respect to a basis ``\{b_i\}^{n}_{i=1}`` are obtained simply by matrix-vector multiplication, i.e.

```math
\hat{η_k} = - B_k \widehat{∇f(x_k)}
```

where ``B_k`` is the matrix representing the operator with respect to the basis ``\{b_i\}^{n}_{i=1}`` and ``\widehat{∇f(x_k)}`` as above. In the end, the search direction ``η_k`` is generated from the coordinates ``\hat{eta_k}`` and the vectors of the basis ``\{b_i\}^{n}_{i=1}`` in both variants.
The [`AbstractQuasiNewtonUpdateRule`](@ref) indicates which quasi-Newton update rule is used. In all of them, the Euclidean update formula is used to generate the matrix ``H_{k+1}`` and ``B_{k+1}``, and the basis ``\{b_i\}^{n}_{i=1}`` is transported into the upcoming tangent space ``T_{x_{k+1}} \mathcal{M}``, preferably with an isometric vector transport, or generated there.

# Fields
* `basis` – the basis.
* `matrix` – the matrix which represents the approximating operator.
* `scale` – indicates whether the initial matrix (= identity matrix) should be scaled before the first update.
* `update` – a [`AbstractQuasiNewtonUpdateRule`](@ref).
* `vector_transport_method` – an `AbstractVectorTransportMethod`

# See also
[`QuasiNewtonLimitedMemoryDirectionUpdate`](@ref)
[`QuasiNewtonCautiousDirectionUpdate`](@ref)
[`AbstractQuasiNewtonDirectionUpdate`](@ref)
"""
mutable struct QuasiNewtonMatrixDirectionUpdate{
    NT<:AbstractQuasiNewtonUpdateRule,
    B<:AbstractBasis,
    VT<:AbstractVectorTransportMethod,
    M<:AbstractMatrix,
} <: AbstractQuasiNewtonDirectionUpdate
    basis::B
    matrix::M
    scale::Bool
    update::NT
    vector_transport_method::VT
end
function QuasiNewtonMatrixDirectionUpdate(
    update::AbstractQuasiNewtonUpdateRule,
    basis::B,
    m::M,
    ;
    scale::Bool=true,
    vector_transport_method::AbstractVectorTransportMethod=ParallelTransport(),
) where {M<:AbstractMatrix,B<:AbstractBasis}
    return QuasiNewtonMatrixDirectionUpdate{
        typeof(update),B,typeof(vector_transport_method),M
    }(
        basis, m, scale, update, vector_transport_method
    )
end
function (d::QuasiNewtonMatrixDirectionUpdate{T})(
    p, o
) where {T<:Union{InverseBFGS,InverseDFP,InverseSR1,InverseBroyden}}
    return get_vector(
        p.M, o.x, -d.matrix * get_coordinates(p.M, o.x, o.∇, d.basis), d.basis
    )
end
function (d::QuasiNewtonMatrixDirectionUpdate{T})(
    p, o
) where {T<:Union{BFGS,DFP,SR1,Broyden}}
    return get_vector(
        p.M, o.x, -d.matrix \ get_coordinates(p.M, o.x, o.∇, d.basis), d.basis
    )
end

@doc raw"""
    QuasiNewtonLimitedMemoryDirectionUpdate <: AbstractQuasiNewtonDirectionUpdate

This [`AbstractQuasiNewtonDirectionUpdate`](@ref) represents the limited-memory Riemanian BFGS update, where the approximating  oprator is represented by ``m`` stored pairs of tangent vectors ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}`` in the ``k``-th iteration.
For the calculation of the search direction ``η_k``, the generalisation of the two-loop recursion is used (see [^HuangGallivanAbsil2015]), since it only requires inner products and linear combinations of tangent vectors in ``T_{x_k} \mathcal{M}``. For that the stored pairs of tangent vectors ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}``, the gradient ``∇f(x_k)`` of the objective function ``f`` in ``x_k`` and the positive definite self-adjoint operator

```math
\mathcal{B}^{(0)}_k[\cdot] = \frac{g_{x_k}(s_{k-1}, y_{k-1})}{g_{x_k}(y_{k-1}, y_{k-1})} \; \mathrm{id}_{T_{x_k} \mathcal{M}}[\cdot]
```

are used. The two-loop recursion can be understood as that the [`InverseBFGS`](@ref) update is executed ``m`` times in a row on ``\mathcal{B}^{(0)}_k[\cdot]`` using the tangent vectors ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}``, and in the same time the resulting operator ``\mathcal{B}^{LRBFGS}_k [\cdot]`` is directly applied on ``∇f(x_k)``.
When updating there are two cases: if there is still free memory, i.e. ``k < m``, the previously stored vector pairs ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}`` have to be transported into the upcoming tangent space ``T_{x_{k+1}} \mathcal{M}``; if there is no free memory, the oldest pair ``\{ \widetilde{s}_{k−m}, \widetilde{y}_{k−m}\}`` has to be discarded and then all the remaining vector pairs ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m+1}^{k-1}`` are transported into the tangent space ``T_{x_{k+1}} \mathcal{M}``. After that we calculate and store ``s_k = \widetilde{s}_k = T^{S}_{x_k, α_k η_k}(α_k η_k)`` and ``y_k = \widetilde{y}_k``. This process ensures that new information about the objective function is always included and the old, probably no longer relevant, information is discarded.

# Fields
* `method` – the maximum number of vector pairs stored.
* `memory_s` – the set of the stored (and transported) search directions times step size ``\{ \widetilde{s}_i\}_{i=k-m}^{k-1}``.
* `memory_y` – set of the stored gradient differences ``\{ \widetilde{y}_i\}_{i=k-m}^{k-1}``.
* `ξ` – a variable used in the two-loop recursion.
* `ρ` – a variable used in the two-loop recursion.
* `scale` –
* `vector_transport_method` – a `AbstractVectorTransportMethod`

# See also
[`InverseBFGS`](@ref)
[`QuasiNewtonCautiousDirectionUpdate`](@ref)
[`AbstractQuasiNewtonDirectionUpdate`](@ref)

[^HuangGallivanAbsil2015]:
    > Huang, Wen and Gallivan, K. A. and Absil, P.-A., A Broyden Class of Quasi-Newton Methods for Riemannian Optimization,
    > SIAM J. Optim., 25 (2015), pp. 1660-1685.
    > doi: [10.1137/140955483](https://doi.org/10.1137/140955483)
"""
mutable struct QuasiNewtonLimitedMemoryDirectionUpdate{
    NT<:AbstractQuasiNewtonUpdateRule,T,VT<:AbstractVectorTransportMethod
} <: AbstractQuasiNewtonDirectionUpdate
    method::NT
    memory_s::CircularBuffer{T}
    memory_y::CircularBuffer{T}
    ξ::Vector{Float64}
    ρ::Vector{Float64}
    scale::Float64
    vector_transport_method::VT
end
function QuasiNewtonLimitedMemoryDirectionUpdate(
    method::NT,
    ::T,
    memory_size::Int;
    scale::Bool=true,
    vector_transport_method::AbstractVectorTransportMethod=ParallelTransport(),
) where {NT<:AbstractQuasiNewtonUpdateRule,T,VT<:AbstractVectorTransportMethod}
    return QuasiNewtonLimitedMemoryDirectionUpdate{NT,T,typeof(vector_transport_method)}(
        method,
        CircularBuffer{T}(memory_size),
        CircularBuffer{T}(memory_size),
        zeros(memory_size),
        zeros(memory_size),
        scale,
        vector_transport_method,
    )
end
function (d::QuasiNewtonLimitedMemoryDirectionUpdate{InverseBFGS})(p, o)
    r = deepcopy(o.∇)
    m = length(d.memory_s)
    m == 0 && return -r
    for i in m:-1:1
        d.ρ[i] = 1 / inner(p.M, o.x, d.memory_s[i], d.memory_y[i]) # 1 sk 2 yk
        d.ξ[i] = inner(p.M, o.x, d.memory_s[i], r) * d.ρ[i]
        r .= r .- d.ξ[i] .* d.memory_y[i]
    end
    r .= 1 / (d.ρ[m] * norm(p.M, o.x, last(d.memory_y))^2) .* r
    for i in 1:m
        r .= r .+ (d.ξ[i] - d.ρ[i] * inner(p.M, o.x, d.memory_y[i], r)) .* d.memory_s[i]
    end
    return -project(p.M, o.x, r)
end

@doc raw"""
    QuasiNewtonCautiousDirectionUpdate <: AbstractQuasiNewtonDirectionUpdate

These [`AbstractQuasiNewtonDirectionUpdate`](@ref)s represent any quasi-Newton update rule, which are based on the idea of a so-called cautious update. The search direction is calculated as given in [`QuasiNewtonMatrixDirectionUpdate`](@ref) or [`LimitedMemoryQuasiNewctionDirectionUpdate`]. But the update given in [`QuasiNewtonMatrixDirectionUpdate`](@ref) or [`LimitedMemoryQuasiNewctionDirectionUpdate`] is only executed if

```math
\frac{g_{x_{k+1}}(y_k,s_k)}{\lVert s_k \rVert^{2}_{x_{k+1}}} \geq \theta(\lVert ∇f(x_k) \rVert_{x_k}),
```

is satisfied, where ``\theta`` is a monotone increasing function satisfying ``\theta(0) = 0`` and ``\theta`` is strictly increasing at ``0``. If this is not the case, the corresponding update will be skipped, which means that for [`QuasiNewtonMatrixDirectionUpdate`](@ref) the matrix ``H_k`` or ``B_k`` is not updated, but the basis ``\{b_i\}^{n}_{i=1}`` is nevertheless transported into the upcoming tangent space ``T_{x_{k+1}} \mathcal{M}``, and for [`LimitedMemoryQuasiNewctionDirectionUpdate`] neither the oldest vector pair ``\{ \widetilde{s}_{k−m}, \widetilde{y}_{k−m}\}`` is discarded nor the newest vector pair ``\{ \widetilde{s}_{k}, \widetilde{y}_{k}\}`` is added into storage, but all stored vector pairs ``\{ \widetilde{s}_i, \widetilde{y}_i\}_{i=k-m}^{k-1}`` are transported into the tangent space ``T_{x_{k+1}} \mathcal{M}``.
If [`InverseBFGS`](@ref) or [`InverseBFGS`](@ref) is chosen as update, then the resulting method follows the method of [^HuangAbsilGallivan2018], taking into account that the corresponding step size is chosen.

# Fields
* `update` – an [`AbstractQuasiNewtonDirectionUpdate`](@ref)
* `θ` – a monotone increasing function satisfying ``θ(0) = 0`` and ``θ`` is strictly increasing at ``0``.

# See also
[`QuasiNewtonMatrixDirectionUpdate`](@ref)
[`QuasiNewtonLimitedMemoryDirectionUpdate`](@ref)

[^HuangAbsilGallivan2018]:
    > Huang, Wen and Absil, P.-A and Gallivan, Kyle, A Riemannian BFGS Method Without Differentiated Retraction for Nonconvex Optimization Problems,
    > SIAM J. Optim., 28 (2018), pp. 470-495.
    > doi: [10.1137/17M1127582](https://doi.org/10.1137/17M1127582)
"""
mutable struct QuasiNewtonCautiousDirectionUpdate{U} <:
               AbstractQuasiNewtonDirectionUpdate where {
    U<:Union{
        QuasiNewtonMatrixDirectionUpdate,QuasiNewtonLimitedMemoryDirectionUpdate{T}
    },
} where {T<:AbstractQuasiNewtonUpdateRule}
    update::U
    θ::Function
end
function QuasiNewtonCautiousDirectionUpdate(
    update::U; θ::Function=x -> x
) where {
    U<:Union{
        QuasiNewtonMatrixDirectionUpdate,QuasiNewtonLimitedMemoryDirectionUpdate{T}
    },
} where {T<:AbstractQuasiNewtonUpdateRule}
    return QuasiNewtonCautiousDirectionUpdate{U}(update, θ)
end
(d::QuasiNewtonCautiousDirectionUpdate)(p, o) = d.update(p, o)

# access the inner vector transport method
function get_update_vector_transport(u::AbstractQuasiNewtonDirectionUpdate)
    return u.vector_transport_method
end
function get_update_vector_transport(u::QuasiNewtonCautiousDirectionUpdate)
    return get_update_vector_transport(u.update)
end


text:
# To use, you must clone directly from github:
# https://github.com/rgiordan/LinearResponseVariationalBayes.jl

module LinearResponseVariationalBayes

push!(LOAD_PATH, joinpath(Pkg.dir("LinearResponseVariationalBayes"), "src"))

# include(joinpath(Pkg.dir("LinearResponseVariationalBayes"), "src", "VariationalModelIndices.jl"))
# include(joinpath(Pkg.dir("LinearResponseVariationalBayes"), "src", "ExponentialFamilies.jl"))
# include(joinpath(Pkg.dir("LinearResponseVariationalBayes"), "src", "SubModels.jl"))
# include(joinpath(Pkg.dir("LinearResponseVariationalBayes"), "src", "HessianReparameterization.jl"))

import VariationalModelIndices
import ExponentialFamilies
import SubModels
import HessianReparameterization

end # module


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# # Ocean Color Data & Model Recipes
#
# Here we set out to compare [CBIOMES model output](https://cbiomes.readthedocs.io/en/latest/) with [ocean color data](https://www.oceancolour.org).
#
# <img src="../figs/cbiomes-01.png" alt="Drawing" style="height: 100px;"/>
#
# _Notebooks are written in [Julia](https://julialang.org) and reproducible via [binder](https://mybinder.org/v2/gh/gaelforget/Cbiomes2019Notebooks/master)._

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Activate packages for later use
#
# It is assumed that listed packages have aleary been installed using `julia`'s package manager (documentation available [here](https://docs.julialang.org/en/)).
# -

using OceanColorData, Plots, Distributions

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Model and data wavebands
#
# Currently, the `OC-CCI` [satellite data set](https://esa-oceancolour-cci.org) provides remotely sensed reflectance at 6 wavelengths (`wv_cci` in `nm`) while the `CBIOMES-global` [ocean model](https://cbiomes.readthedocs.io/) outputs irradiance reflectance at 13 wavelengths (`wv_drwn3` in `nm`).
# -

wv_cci=[412, 443, 490, 510, 555, 670]
wv_drwn3=[400,425,450,475,500,525,550,575,600,625,650,675,700];

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Define interpolation factors
#
# Later on, `jj` and `ww` are used to interpolate, in wavelength space, model output from `wv_drwn3` to `wv_cci`.

# + {"slideshow": {"slide_type": "-"}}
jj=Array{Int64,1}(undef,6)
ww=Array{Float64,1}(undef,6)
for ii=1:6
    tmp=wv_cci[ii].-wv_drwn3
    kk=maximum(findall(tmp.>=0))
    jj[ii]=kk
    ww[ii]=tmp[kk]/(wv_drwn3[kk+1]-wv_drwn3[kk])
end

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Define a test case
#
# A vector of 13 irradiance reflectances (`Rirr`) is used to represent one space-time location in the model. Later on, we derive `Rrs0` and `Rrs` from `Rirr`. The expected results of the recipe are `ref_Rrs0` and `ref_Rrs`.

# + {"slideshow": {"slide_type": "-"}}
siz=[1,1]

Rirr=Array{Float32,3}(undef,(siz[1],siz[2],13))
Rirr[1,1,:]= 1e-3*[23.7641,26.5037,27.9743,30.4914,28.1356,
    21.9385,18.6545,13.5100,5.6338,3.9272,2.9621,2.1865,1.8015]

ref_Rrs0=Array{Float32,3}(undef,(siz[1],siz[2],13))
ref_Rrs0=1e-3*[4.1753, 4.6640, 4.9270, 5.3781, 4.9558, 3.8505
    , 3.2680, 2.3598, 0.9796, 0.6822, 0.5143, 0.3795, 0.3126]

ref_Rrs=Array{Float32,3}(undef,(siz[1],siz[2],6))
ref_Rrs=1e-3*[4.4099, 4.8533, 5.1247, 4.5137, 3.0864, 0.4064];

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Convert to remotely sensed reflectance
#
# The following recipe is from [Dutkiewicz et al 2018](https://doi.org/10.5194/bg-15-613-2018).
# -

tmp=Rirr/3
Rrs0=(0.52*tmp)./(1.0 .-1.7*tmp);

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Interpolate in wavelength space
#
# Interpolating model output from `wv_drwn3` to `wv_cci` allows for direct comparison with satellite data.

# + {"slideshow": {"slide_type": "-"}}
Rrs=Array{Float32,3}(undef,(siz[1],siz[2],6))
for vv=1:6
    tmp0=Rrs0[:,:,jj[vv]]
    tmp1=Rrs0[:,:,jj[vv]+1]
    Rrs[:,:,vv]=tmp0.*(1-ww[vv])+tmp1.*ww[vv]
end

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Verify result
#
# Let's visualize using the `Plots.jl` package that the resulting `Rrs` matches `ref_Rrs`.

# +
plot(vec(Rrs),linewidth=4,lab="recomputed Rrs")
plot!(ref_Rrs,linewidth=4,ls=:dash,lab="reference result")

#Or, equivalentaly, via OceanColorData.RemotelySensedReflectance :
RrsBis=RemotelySensedReflectance(Rirr,wv_drwn3,wv_cci)
scatter!(vec(RrsBis),lab="OceanColorData.jl")

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Estimate chlorophyll  from reflectances
#
# Satellite `Chl_a` estimates provided by `OC-CCI` derive from `Rrs` using the blue/green reflectance ratio method as done with model output in the next code bloc (See [Dutkiewicz et al 2018](https://doi.org/10.5194/bg-15-613-2018) for details).

# +
RRSB=max.(Rrs[:,:,2],Rrs[:,:,3]) #blue
RRSG=Rrs[:,:,5] #green
X = log10.(RRSB./RRSG) #ratio of blue to green

C=[0.3272, -2.9940, +2.7218, -1.2259, -0.5683] #OC4 algorithms (SeaWifs, CCI)
a0=C[1]; a1=C[2]; a2=C[3]; a3=C[4]; a4=C[5];
chld=exp10.(a0.+a1*X+a2*X.^2+a3*X.^3+a4*X.^4); #apply polynomial recipe

#Or, equivalently, via OceanColorData.RrsToChla
[chld[1] RrsToChla(vec(Rrs))]

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Optical classification using reflectances
#
# `Fuzzy logic` classifiers defined in [Moore et al 2009](https://doi.org/10.1016/j.rse.2009.07.016) and [Jackson et al 2017](http://dx.doi.org/10.1016/j.rse.2017.03.036) can be used to assign optical class memberships from an `Rrs` vector. While Moore et al define `n=8` classes using an in-situ database, Jackson et al instead define `n=14` classes using a satellite database. The latter benefits from better data coverage across all of the ecological provinces of the global ocean and is used in `OC-CCI`.
#
# In both cases the classifier is encoded in a mean reflectance spectra (`M[i][1:6]`) and a covariance matrix (`S[i][1:6,1:6]`) provided for each optical class (`i` in `1:n`). Class memberships are then derived by computing the squared Mahalanobis distance to each `M[i]` and passing the result to cumulative chi-squared distribution function (Equations 11 and 12 in [Moore et al 2011](https://doi.org/10.1109/36.942555)).

# + {"slideshow": {"slide_type": "subslide"}}
(M,Sinv)=Moore2009()
M09=Dict("M" => M, "Sinv" => Sinv)
plot(wv_cci,M,w=3); xlabel!("nm"); ylabel!("Rrs")

# + {"slideshow": {"slide_type": "subslide"}}
(M,Sinv)=Jackson2017()
J17=Dict("M" => M, "Sinv" => Sinv)
plot(wv_cci,M,w=3); xlabel!("nm"); ylabel!("Rrs")

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Compute class memberships
#
# The `fcm` function below takes in a classifier matrix (`M` + `Sinv`) along with a vector of remotely sensed reflectcances (`Rrs`) as input and returns a membership vector (values between 0 and 1).

# + {"slideshow": {"slide_type": "subslide"}}
M09["Membership"]=FuzzyClassification(M09["M"],M09["Sinv"],Rrs)
bar(M09["Membership"]); xlabel!("class"); ylabel!("M09 membership")

# + {"slideshow": {"slide_type": "subslide"}}
J17["Membership"]=FuzzyClassification(J17["M"],J17["Sinv"],Rrs)
bar(J17["Membership"]); xlabel!("class"); ylabel!("J17 membership")

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
using FractionalCalculus
using Plots

tspan=collect(0:0.01:3)

exampleintegral = []

for i in range(0.1, 0.9, step=0.1)
    push!(exampleintegral, fracint(x->x, i, tspan, 0.001, RLIntApprox()))
end

plot(tspan, exampleintegral, title="Different order of integral", linewidth=1, label=["α=0.1" "α=0.2" "α=0.3" "α=0.4" "α=0.5" "α=0.6" "α=0.7" "α=0.8" "α=0.9" "α=1"], legend=:bottomright)


text:
# ImagePFM.jl
# FileIO.jl integration for PFM image files.
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

"""
    InvalidPfmFileFormat <: Exception

The parsed PFM file is not conform to standard or is corrupted. Read `msg` for root cause.
"""
struct InvalidPfmFileFormat <: Exception
    msg::String
end

Base.showerror(io::IO, e::InvalidPfmFileFormat) = print(io, typeof(e), ": ", e.msg)

"""
    InvalidRGBEltype <: Exception

The eltype of RGB can only be Float32.
"""
struct InvalidRGBEltype <: Exception end

Base.showerror(io::IO, e::InvalidRGBEltype) = print(io, typeof(e), ": eltype of RGB can only be Float32")


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

@recipe function f(c::Composition{D}) where {D}
  p = parts(c)
  w = components(c)

  # handle missing components
  s = string.(p)
  m = ismissing.(w)
  x = [m[i] ? s[i]*" (missing)" : s[i] for i in 1:D]
  y = [m[i] ? NaN : w[i] for i in 1:D]

  seriestype --> :bar
  legend --> false
  xticks --> (1:length(x), x)
  xrotation --> 30
  yguide --> "components"

  y
end


text:
function collect_subnet_results(results_dict, objective, gen_areas, gen_zones, output_folder, plot_best_x)
   # upscale = 8 #upscaling in resolution
   # fntsm = font("serif", pointsize=round(8.0*upscale))
   # fntlg = font("serif", pointsize=round(12.0*upscale))
   # default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
   # default(size=(800*upscale,600*upscale)) #Plot canvas size
   # default(left_margin = 15mm, right_margin = 15mm)
   # default(dpi=300) #Only for PyPlot - presently broken

   # pyplot(guidefont=font, xtickfont=font, ytickfont=font, legendfont=font)
   # plotlyjs()

   # gr(titlefont=font("serif", pointsize=round(10.0*upscale)), guidefont=font("serif", pointsize=round(10.0*upscale)), tickfont=font("serif", pointsize=round(10.0*upscale)), legendfont=font("serif", pointsize=round(10.0*upscale)4))
   # pyplot()
   # gr()
   # Plot the primary label first
   if objective == "mincost"
      first_label = "cost (USD)"
   elseif objective == "areagen"
      n_areas = length(gen_areas)
      if n_areas > 1
         gen_areas_sorted = sort(gen_areas)
         first_label = "generation in areas"
         for i in 1:(n_areas-2)
            first_label *= " $(gen_areas_sorted[i]),"
         end
         first_label *=  " $(gen_areas_sorted[end-1]) and $(gen_areas_sorted[end]) (p.u.)"
      else
         first_label = "generation in areas $(gen_areas[1]) (p.u.)"
      end
   elseif objective == "zonegen"
      n_zones = length(gen_zones)
      if n_zones > 1
         gen_zones_sorted = sort(gen_zones)
         first_label = "generation in zones"
         for i in 1:(n_zones-2)
            first_label *= " $(gen_zones_sorted[i]),"
         end
         first_label *=  " $(gen_zones_sorted[end-1]) and $(gen_zones_sorted[end]) (p.u.)"
      else
         first_label = "generation in zones $(gen_zones[1]) (p.u.)"
      end
   end

   label = first_label
   results = results_dict[label]
   # Remove base value to plot separately
   br_results = Array{Float64,1}()
   br_indices = Array{String,1}()
   br_results_fixf = Array{Float64,1}()
   br_indices_fixf = Array{String,1}()
   br_results_indirPQ = Array{Float64,1}()
   br_indices_indirPQ = Array{String,1}()
   base_res = -9999.0
   subnet_array = []
   n_subnets = 0
   for (folder,res) in results
      println("folder: $folder")
      if results_dict["status"][folder] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
         subnet_array = union(subnet_array, results_dict["subnet"][folder]) # Get maximal set with subnets represented in any folders
         # n_subnets = max(n_subnets,length(res))
         # if subnet_array == nothing
         #    subnet_array = results_dict["subnet"][folder]
         #    println("subnet_array: $subnet_array")
         # else
         #    @assert (subnet_array == results_dict["subnet"][folder]) "Subnet order does not match in results files: $subnet_array vs. $(results_dict["subnet"][folder])"
         # end
         if occursin("br", folder)
            if occursin("_fixf", folder)
               push!(br_results_fixf, sum(res))
               push!(br_indices_fixf, folder[3:(end-5)])
            elseif occursin("_indirPQ", folder)
               push!(br_results_indirPQ, sum(res))
               push!(br_indices_indirPQ, folder[3:(end-8)])
            else
               push!(br_results, sum(res))
               push!(br_indices, folder[3:end])
            end
         # elseif folder == "base"
         #    base_res = res
         elseif folder == "base_unconstrPS"
            base_res = res
         end
      else
         subnet_array = union(subnet_array, results_dict["subnet"][folder]) # Get maximal set with subnets represented in any folders
         if occursin("br", folder)
            if occursin("_fixf", folder)
               push!(br_results_fixf, -9999.0)
               push!(br_indices_fixf, folder[3:(end-5)])
            elseif occursin("_indirPQ", folder)
               push!(br_results_indirPQ, -9999.0)
               push!(br_indices_indirPQ, folder[3:(end-8)])
            else
               push!(br_results, -9999.0)
               push!(br_indices, folder[3:end])
            end
         # elseif folder == "base"
         #    base_res = -9999.0
         elseif folder == "base_unconstrPS"
            base_res = -9999.0
         end
      end
      n_subnets = length(subnet_array)
   end
   # println("br_results: $br_results")
   # Sort values
   br_results_tosort = copy(br_results)
   br_results_tosort[br_results_tosort.==-9999.0] .= Inf
   perm = sortperm(br_results_tosort)
   idx_sorted = br_indices[perm]
   # println("perm: $perm")
   # println("idx_sorted: $idx_sorted")

   # plot all results
   plot_results_dict_bar(results_dict,n_subnets,subnet_array,idx_sorted,output_folder,plot_best_x)

end


text:
#+Tag L:target
# This is an enum example [[target][]]
@enum Alpha A  B=1 C