text:
# SPDX-License-Identifier: MIT

struct NormalWishart{T1<:AbstractVector{<:AbstractFloat},T2<:AbstractFloat,
                     T3<:AbstractMatrix{<:AbstractFloat},T4<:AbstractFloat} <: ExponentialFamilyDistribution
	μ::T1
	β::T2
	W::T3
	ν::T4
end

function η(p::NormalWishart)
	D = length(p.μ)
    μ, β, W, ν = p.μ, p.β, p.W, p.ν

	M =
	η₁ = β * μ
	η₂ = -β/2
	η₃ = -(1/2)*vec(inv(W) + β * μ*μ')
	η₄ =  (ν - D) / 2
	η = vcat(η₁, η₂, η₃, η₄)
end

function ξ(p::NormalWishart, η)
	D = length(p.μ)
	β = -2*η[D+1]
	μ = (1/β)*η[1:D]
	ν = 2*η[end] + D
    W = inv(Symmetric(reshape(-2*η[D+2:end-1], D, D) - β*(μ*μ')))
	vcat(μ, β, vec(W), ν)
end

function unpack(p::NormalWishart, μ)
	D = length(p.μ)
	Yx = μ[1:D]
	xᵀYx = μ[D+1]
    Y = reshape(μ[D+2:end-1], D, D)
	logdetY = μ[end]
	(Yx, xᵀYx, Y, logdetY)
end

function A(p::NormalWishart, η)
	D = length(p.μ)

	# We don't use `p.W` or `p.ν` to be able to calculate the gradient w.r.t. η.
    β = -2*η[D+1]
	μ = (1/β)*η[1:D]
	ν = 2*η[end] + D
    W⁻¹ = reshape(-2*η[D+2:end-1], D, D) - β*μ*μ'

	(
		- (D/2)*log(β)
		+ (ν/2)*(-logdet(W⁻¹) + D*log(2))
		+ sum(loggamma.((ν + 1 .- (1:D))/2))
	)
end

function sample(p::NormalWishart)
    μ, β, W, ν = p.μ, p.β, p.W, p.ν
    Λ = rand(Distributions.Wishart(ν, W))
    m = rand(Distributions.MultivariateNormal(m, inv(β*Λ)))
	m, Λ
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
using MultipleScattering
using Plots; pyplot()

function box_size_convergence(m=4, volfrac = 0.05,
    radius = 1.0, times = [20.0,30.0,40.0,50.0,60.0,70.0], k_arr=collect(LinRange(0.01,1.0,100)) )

    listener_position = [-10.0,0.0]
    bigshape = TimeOfFlight(listener_position,maximum(times))

    seed = MersenneTwister(1).seed
    allparticles = random_particles(volfrac, radius, bigshape; seed = seed)

    simulations = map(times) do t
        println("Calculating response with particles at a maximum of $t away")
        shape = TimeOfFlight(listener_position,t)
        particles = filter(p->p⊆shape, allparticles)
        return FrequencySimulation(particles, k_arr; seed=seed, hankel_order=m)
    end

    return map(m->TimeSimulation(m;time_arr=LinRange(0.0,100.0,201)),simulations)
end

function plot_box_size_convergence(simulations;times = [20,30,40,50,60,70])

    colors = LinRange(RGB(0.6,1,0.6),RGB(0,0.4,0),length(times))

    plot(xguide ="Time (t)",yguide ="Response")
    for s in eachindex(simulations)
        plot!(simulations[s],color=colors[s],label="Box cut at t=$(times[s])")
    end
    plot!(title="Time response from random particles from increasingly large boxes")
end

simulations = box_size_convergence()
plot_box_size_convergence(simulations)

ext: jl


text:
@inline function mulby2(a::Float120)
    Float120(a.hi*2.0, a.lo*2.0)
end

@inline function divby2(a::Float120)
    Float120(a.hi*0.5, a.lo*0.5)
end

@inline function mulbypow2(a::Float120,p::Float64)
    Float120(a.hi*p, a.lo*p)
end

@inline function divbypow2(a::Float120,p::Float64)
    fr,xp = frexp(p)
    mulbypow2(a, ldexp(fr,-xp))
end

function div{T<:Float120}(a::T,b::T)
    if (b.hi == zero(Float64))
        throw(DomainError("denominator must be nonzero"))
    end
    temp = a/b
    trunc(temp)
end
div(a::Float120,b::Float64) = div(a,Float120(b))
div(a::Float64,b::Float120) = div(Float120(a),b)

function fld{T<:Float120}(a::T,b::T)
    if (b.hi == zero(Float64))
        throw(DomainError("denominator must be nonzero"))
    end
    floor( a/b )
end
fld(a::Float120,b::Float64) = fld(a,Float120(b))
fld(a::Float64,b::Float120) = fld(Float120(a),b)

function cld{T<:Float120}(a::T,b::T)
    if (b.hi == zero(Float64))
        throw(DomainError("denominator must be nonzero"))
    end
    ceil( a/b )
end
cld(a::Float120,b::Float64) = cld(a,Float120(b))
cld(a::Float64,b::Float120) = cld(Float120(a),b)

function rem{T<:Float120}(a::T,b::T)
    if (b.hi == zero(Float64))
        throw(DomainError("denominator must be nonzero"))
    end
    a - div(a,b)*b
end

function mod{T<:Float120}(a::T,b::T)
    if (b.hi == zero(Float64))
        throw(DomainError("modulus must be nonzero"))
    end
    if signbit(a.hi)==signbit(b.hi)
        rem(a,b)
    else
        d = floor(a/b)
        a - d*b
    end
end
mod(a::Float120,b::Float64) = mod(a,convert(Float120,b))
mod(a::Float64,b::Float120) = mod(convert(Float120,a),b)

function divrem{T<:Float120}(a::T,b::T)
    if (b.hi == zero(Float64))
        throw(DomainError("denominator must be nonzero"))
    end
    temp = a/b
    d = trunc(temp)
    r = a - d*b
    d,r
end
divrem(a::Float120, b::Float64) = divrem(a,convert(Float120,b))
divrem(a::Float64, b::Float120) = divrem(convert(Float120,a),b)

function fldmod{T<:Float120}(a::T,b::T)
    d = floor(a/b)
    a - d*b
    d,a
end
fldmod(a::Float120, b::Float64) = fldmod(a,convert(Float120,b))
fldmod(a::Float64, b::Float120) = fldmod(convert(Float120,a),b)

#=
   This well-behaved bounded modulo implementation is from
   The pitfalls of verifying floating-point computations
   by David Monniaux, 2008
   http://arxiv.org/abs/cs/0701192v5
=#
function modulo{T<:Float120}(a::T, lowerbound::T, upperbound::T)
    delta = upperbound - lowerbound
    a - (floor((a - lowerbound)/delta) * delta)
end
modulo(a::Float120, lowerbound::Float64, upperbound::Float64) = modulo(a,convert(Float120,lowerbound),convert(Float120,upperbound))
modulo(a::Float64, lowerbound::Float120, upperbound::Float120) = modulo(convert(Float120,a),lowerbound,upperbound)


text:
module TrigPolys

using FFTW

export TrigPoly, random_trig_poly, evaluate, evaluateT, interpolatev, interpolate
export pad_to, truncate

"""
    struct TrigPoly{T<:AbstractFloat, VT<:AbstractVector{T}}
        a0::T    # Constant coefficient
        ac::VT   # cos coefficients
        as::VT   # sin coefficients
    end

Represents a trigonometric polynomial by its coefficients.
"""
struct TrigPoly{T<:AbstractFloat, VT<:AbstractVector{T}}
    a0::T    # Constant coefficient
    ac::VT   # cos coefficients
    as::VT   # sin coefficients
    function TrigPoly(a0::S, ac::VS, as::VS) where {S<:AbstractFloat, VS<:AbstractVector{S}}
        @assert length(ac) == length(as) "sin and cos coefficients must have same length!"
        new{S, VS}(a0, ac, as)
    end
end

TrigPoly(a0::S, ac::VT, as::VU) where {
    S<:Number, T<:Number, U<:Number, VT<:AbstractVector{T}, VU<:AbstractVector{U}
} = TrigPoly(float(a0), float(ac), float(as))

TrigPoly(a::AbstractFloat) = TrigPoly(a, typeof(a)[], typeof(a)[])

"""
    TrigPoly(a::Number)

Converts a number `a` into a TrigPoly with a constant coefficient.
`a` will be converted to a float.
"""
TrigPoly(a::Number) = TrigPoly(float(a))

"""
    a0(p::TrigPoly)

Equivalent to `p.a0`, used for incorporation into automatic differentiation libraries.
"""
a0(p::TrigPoly) = p.a0

"""
    ac(p::TrigPoly)

Equivalent to `p.ac`, used for incorporation into automatic differentiation libraries.
"""
ac(p::TrigPoly) = p.ac

"""
    as(p::TrigPoly)

Equivalent to `p.as`, used for incorporation into automatic differentiation libraries.
"""
as(p::TrigPoly) = p.as

function Base.getproperty(p::TrigPoly, v::Symbol)
    if v == :n
        return length(p.ac)
    else
        return getfield(p, v)
    end
end

"""
    n(p::TrigPoly)

Equivalent to `p.n`. Returns the length of cosine/sine coefficients,
which is `length(p.ac) == length(p.as)`.
"""
n(p::TrigPoly) = p.n

"""
    degree(p::TrigPoly)

Returns the number of coefficients of `p`, which is `2*p.n + 1`.
"""
degree(p::TrigPoly) = 2 * p.n + 1

function get_mn(u::AbstractVector)
    m = length(u)
    @assert m % 2 == 1 "Only odd-length vectors supported!"
    n = (m - 1)÷2
    (m,n)
end

"""
    TrigPoly(u::AbstractVector{T}) where {T<:AbstractFloat}

Constructs a `TrigPoly` given a vector of coefficients of length `2n+1`. The
first entry of `u` is `a0`, the next `n` entries are `ac` and the last `n`
entries are `as`.
"""
function TrigPoly(u::AbstractVector{T}) where {T<:AbstractFloat}
    m, n = get_mn(u)
    TrigPoly(u[1], u[2:n+1], u[n+2:m])
end

"""
    Base.vec(p::TrigPoly)

Converts a TrigPoly into a vector, inverse operation of
[`TrigPolys.TrigPoly`](@ref TrigPolys.TrigPoly(u::AbstractVector{T}) where {T<:AbstractFloat})
"""
Base.vec(p::TrigPoly) = [a0(p); ac(p); as(p)]

"""
    Base.:(==)(p1::TrigPoly, p2::TrigPoly)

Returns true if the coefficients of `p1` and `p2` are all equal.
"""
function Base.:(==)(p1::TrigPoly, p2::TrigPoly)
    p1.n == p2.n && a0(p1) == a0(p2) && ac(p1) == ac(p2) && as(p1) == as(p2)
end

"""
    Base.:+(p1::TrigPoly, p2::TrigPoly)

Adds `p1` to `p2`. The degree of the resulting polynomial will be the maximum
degree of `p1` and `p2`: `(p1 + p2).n == max(p1.n, p2.n)`
"""
function Base.:+(p1::TrigPoly, p2::TrigPoly)
    n = max(p1.n, p2.n)
    p1 = pad_to(p1, n)
    p2 = pad_to(p2, n)
    TrigPoly(a0(p1)+a0(p2), ac(p1)+ac(p2), as(p1)+as(p2))
end
Base.:+(p::TrigPoly, a::Number)     = +(promote(p, a)...)
Base.:+(a::Number, p::TrigPoly)     = p + a

Base.:-(p::TrigPoly)                = TrigPoly(-a0(p), -ac(p), -as(p))
Base.:-(p1::TrigPoly, p2::TrigPoly) = p1 + (-p2)
Base.:-(p::TrigPoly, a::Number)     = -(promote(p, a)...)
Base.:-(a::Number, p::TrigPoly)     = -(promote(p, a)...)

"""
    Base.:*(p1::TrigPoly, p2::TrigPoly)

Multiplies `p1` and `p2`. The degree of the resulting polynomial will be the sum of
degrees of `p1` and `p2`: `(p1 * p2).n == p1.n + p2.n`
"""
function Base.:*(p1::TrigPoly, p2::TrigPoly)
    n = p1.n + p2.n
    interpolate(evaluate(pad_to(p1, n)) .* evaluate(pad_to(p2, n)))
end
Base.:*(a::Number, p::TrigPoly)     = *(promote(p, a)...)
Base.:*(p::TrigPoly, a::Number)     = *(promote(p, a)...)

"""
    Base.:/(p::TrigPoly, a::Number)

Divides `p` by a scalar.
"""
Base.:/(p::TrigPoly, a::Number)     = p * (1/a)

Base.promote_rule(::Type{TrigPoly{S, VS}}, ::Type{T}) where {S<:AbstractFloat,VS<:AbstractVector{S},T<:Number} = TrigPoly{S, VS}
Base.convert(::Type{TrigPoly{S, VS}}, x::T) where {S<:AbstractFloat,VS<:AbstractVector{S},T<:Number} = TrigPoly(S(x))

"""
    pad_to(p::TrigPoly{T,VT}, n::Integer) where {T<:AbstractFloat, VT<:AbstractVector{T}}

Increase the length of `p.ac` and `p.as` to `n`, padding with zeros if necessary.
`n` cannot be smaller than `p.n`.
"""
function pad_to(p::TrigPoly{T,VT}, n::Integer) where {T<:AbstractFloat, VT<:AbstractVector{T}}
    @assert n >= p.n "Cannot pad to smaller degree!"
    TrigPoly(a0(p), VT([ac(p); zeros(n-p.n)]), VT([as(p); zeros(n-p.n)]))
end

function pad_to(v, n::Integer)
    xn = length(v)
    xs = (xn - 1) ÷ 2
    [v[1:xs+1]; zeros(n - xs); v[xs+2:xn]; zeros(n - xs)]
end

"""
    pad_by(p::TrigPoly, n::Integer)

Increase the length of `p.ac` and `p.as` by `n`, padding with `n` zeros.
"""
pad_by(p::TrigPoly, n::Integer) = pad_to(p, p.n + n)

"""
    truncate(p::TrigPoly, n::Integer)

Reduce the length of `p.ac` and `p.as` to `n` by deleting coefficients.
`n` cannot be larger than `p.n`.
"""
function truncate(p::TrigPoly, n::Integer)
    @assert n <= p.n "Cannot truncate to higher degree!"
    TrigPoly(a0(p), ac(p)[1:n], as(p)[1:n])
end

"""
    (p::TrigPoly)(x::Number)

Evaluate `p(x)` at a single point. See [TrigPolys.evaluate](@red) for faster
evaluation for special points on the unit circle.
"""
function (p::TrigPoly)(x::Number)
    a0(p) + sum([ac(p)[k]*cos(k*x) + as(p)[k]*sin(k*x) for k in 1:p.n])
end

"""
    random_trig_poly(n)

Generate a `TrigPoly` with random coefficients.
"""
function random_trig_poly(n)
    TrigPoly(randn(), randn(n), randn(n))
end

"""
    basis(n, x)

Generate the basis for `TrigPoly` expressed as a vector, so that
`basis(p.n, x)'*vec(p) = p(x)`.
"""
basis(n, x) = vcat([1], [cos(k*x) for k in 1:n], [sin(k*x) for k in 1:n])

"""
    evaluate(u::AbstractVector)

Evaluates an array representation of `p::TrigPoly`, returned by `vec(p)`.
"""
function evaluate(u::AbstractVector)
    m, n = get_mn(u)
    out = similar(u, Complex{Float64})
    out[1] = u[1]
    for (i, (up, un)) in enumerate(zip(2:n+1, n+2:m))
       out[1+i]   = (u[up]+im*u[un])/2
       out[m+1-i] = (u[up]-im*u[un])/2
    end
    fft!(out)
    real(out)
    # Slower non-inplace version
    #   u0 = u[1]
    #   up = u[2:n+1]
    #   un = u[n+2:m]
    #   pos = (up-im*un)./2
    #   neg = (up+im*un)./2
    #   real(fft([u0; neg; reverse(pos)]))'
end

"""
    evaluate(p::TrigPoly)

Evaluates `p` on `degree(p)` uniformly-spaced points on the circle. See
[here](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Trigonometric_interpolation_polynomial)
for more details.
"""
evaluate(p::TrigPoly) = evaluate(vec(p))

"""
    evaluateT(u::AbstractVector)

Adjoint of evaluate operator on an array representation of `p::TrigPoly`
, returned by `vec(p)`.
"""
function evaluateT(u::AbstractVector)
    m, n = get_mn(u)
    x = fft(u)
    [real(x[1:n+1]); reverse(imag(x[n+2:m]))]
end

"""
    evaluateT(p::TrigPoly)

Adjoint of evaluate operator. Returns a vector.
"""
evaluateT(p::TrigPoly) = evaluateT(vec(p))

"""
    interpolatev(u::Vector)

Inverse of [`TrigPolys.evaluate(p::AbstractVector)`](@ref). Returns a vector.
"""
function interpolatev(u::Vector)
    m, n = get_mn(u)
    out = evaluateT(u::Vector)
    for i = 1:m
        out[i] /= n+.5
    end
    out[1] /= 2
    out
end

"""
    interpolate(u::Vector)

Inverse of [`TrigPolys.evaluate(p::TrigPoly)`](@ref). Returns a `TrigPoly`.
"""
interpolate(u::Vector) = TrigPoly(interpolatev(u))

end


text:
module TestGraphs

using Test

@testset "BasicGraphs" begin
  include("BasicGraphs.jl")
end

@testset "BipartiteGraphs" begin
  include("BipartiteGraphs.jl")
end

@testset "PropertyGraphs" begin
  include("PropertyGraphs.jl")
end

@testset "GraphAlgorithms" begin
  include("GraphAlgorithms.jl")
end

@testset "GraphGenerators" begin
  include("GraphGenerators.jl")
end

end


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

module PkgToMaxSumInterface

using ...Types, ...Query, ..VersionWeights

export Interface, compute_output_dict, greedysolver,
       verify_solution, enforce_optimality!

# A collection of objects which allow interfacing external (Pkg) and
# internal (MaxSum) representation
mutable struct Interface
    # requirements and dependencies, in external representation
    reqs::Requires
    deps::Dict{String,Dict{VersionNumber,Available}}

    # packages list
    pkgs::Vector{String}

    # number of packages
    np::Int

    # states per package: one per version + uninstalled
    spp::Vector{Int}

    # pakage dict: associates an index to each package name
    pdict::Dict{String,Int}

    # package versions: for each package, keep the list of the
    #                   possible version numbers; this defines a
    #                   mapping from version numbers of a package
    #                   to indices
    pvers::Vector{Vector{VersionNumber}}

    # versions dict: associates a version index to each package
    #                version; such that
    #                  pvers[p0][vdict[p0][vn]] = vn
    vdict::Vector{Dict{VersionNumber,Int}}

    # version weights: the weight for each version of each package
    #                  (versions include the uninstalled state; the
    #                   higher the weight, the more favored the version)
    vweight::Vector{Vector{VersionWeight}}

    function Interface(reqs::Requires, deps::Dict{String,Dict{VersionNumber,Available}})
        # generate pkgs
        pkgs = sort!(String[keys(deps)...])

        np = length(pkgs)

        # generate pdict
        pdict = Dict{String,Int}(pkgs[i] => i for i = 1:np)

        # generate spp and pvers
        spp = Vector{Int}(uninitialized, np)

        pvers = [VersionNumber[] for i = 1:np]

        for (p,depsp) in deps, vn in keys(depsp)
            p0 = pdict[p]
            push!(pvers[p0], vn)
        end
        for p0 = 1:np
            sort!(pvers[p0])
            spp[p0] = length(pvers[p0]) + 1
        end

        # generate vdict
        vdict = [Dict{VersionNumber,Int}() for p0 = 1:np]
        for (p,depsp) in deps
            p0 = pdict[p]
            vdict0 = vdict[p0]
            pvers0 = pvers[p0]
            for vn in keys(depsp)
                for v0 in 1:length(pvers0)
                    if pvers0[v0] == vn
                        vdict0[vn] = v0
                        break
                    end
                end
            end
        end

        ## generate wveights:
        vweight = Vector{Vector{VersionWeight}}(uninitialized, np)
        for p0 = 1:np
            pvers0 = pvers[p0]
            spp0 = spp[p0]
            vweight0 = vweight[p0] = Vector{VersionWeight}(uninitialized, spp0)
            for v0 = 1:spp0-1
                vweight0[v0] = VersionWeight(pvers0[v0])
            end
            vweight0[spp0] = VersionWeight(v"0") # last version means uninstalled
        end

        return new(reqs, deps, pkgs, np, spp, pdict, pvers, vdict, vweight)
    end
end

# The output format is a Dict which associates a VersionNumber to each installed package name
function compute_output_dict(sol::Vector{Int}, interface::Interface)
    pkgs = interface.pkgs
    np = interface.np
    pvers = interface.pvers
    spp = interface.spp

    want = Dict{String,VersionNumber}()
    for p0 = 1:np
        p = pkgs[p0]
        s = sol[p0]
        if s != spp[p0]
            v = pvers[p0][s]
            want[p] = v
        end
    end

    return want
end

# Produce a trivial solution: try to maximize each version;
# bail out as soon as some non-trivial requirements are detected.
function greedysolver(interface::Interface)
    reqs = interface.reqs
    deps = interface.deps
    spp = interface.spp
    pdict = interface.pdict
    pvers = interface.pvers
    np = interface.np

    # initialize solution: all uninstalled
    sol = [spp[p0] for p0 = 1:np]

    # set up required packages to their highest allowed versions
    for (rp,rvs) in reqs
        rp0 = pdict[rp]
        # look for the highest version which satisfies the requirements
        rv = spp[rp0] - 1
        while rv > 0
            rvn = pvers[rp0][rv]
            rvn ∈ rvs && break
            rv -= 1
        end
        @assert rv > 0
        sol[rp0] = rv
    end

    # we start from required packages and explore the graph
    # following dependencies
    staged = Set{String}(keys(reqs))
    seen = copy(staged)

    while !isempty(staged)
        staged_next = Set{String}()
        for p in staged
            p0 = pdict[p]
            @assert sol[p0] < spp[p0]
            vn = pvers[p0][sol[p0]]
            a = deps[p][vn]

            # scan dependencies
            for (rp,rvs) in a.requires
                rp0 = pdict[rp]
                # look for the highest version which satisfies the requirements
                rv = spp[rp0] - 1
                while rv > 0
                    rvn = pvers[rp0][rv]
                    rvn ∈ rvs && break
                    rv -= 1
                end
                # if we found a version, and the package was uninstalled
                # or the same version was already selected, we're ok;
                # otherwise we can't be sure what the optimal configuration is
                # and we bail out
                if rv > 0 && (sol[rp0] == spp[rp0] || sol[rp0] == rv)
                    sol[rp0] = rv
                else
                    return (false, Int[])
                end

                rp ∈ seen || push!(staged_next, rp)
            end
        end
        union!(seen, staged_next)
        staged = staged_next
    end

    @assert verify_solution(sol, interface)

    return true, sol
end

# verifies that the solution fulfills all hard constraints
# (requirements and dependencies)
function verify_solution(sol::Vector{Int}, interface::Interface)
    reqs = interface.reqs
    deps = interface.deps
    spp = interface.spp
    pdict = interface.pdict
    pvers = interface.pvers
    vdict = interface.vdict

    # verify requirements
    for (p,vs) in reqs
        p0 = pdict[p]
        sol[p0] != spp[p0] || return false
        vn = pvers[p0][sol[p0]]
        vn ∈ vs || return false
    end

    # verify dependencies
    for (p,d) in deps
        p0 = pdict[p]
        vdict0 = vdict[p0]
        for (vn,a) in d
            v0 = vdict0[vn]
            if sol[p0] == v0
                for (rp, rvs) in a.requires
                    p1 = pdict[rp]
                    if sol[p1] == spp[p1]
                        println("""
                                VERIFICATION ERROR: REQUIRED DEPENDENCY NOT INSTALLED
                                    package p=$p (p0=$p0) version=$vn (v0=$v0) requires package rp=$rp in version set rvs=$rvs
                                    but package $rp is not being installed (p1=$p1 sol[p1]=$(sol[p1]) == spp[p1]=$(spp[p1]))
                                """)
                        return false
                    end
                    vn1 = pvers[p1][sol[p1]]
                    if vn1 ∉ rvs
                        println("""
                                VERIFICATION ERROR: INVALID VERSION
                                    package p=$p (p0=$p0) version=$vn (v0=$v0) requires package rp=$rp in version set rvs=$rvs
                                    but package $rp version is being set to $vn1 (p1=$p1 sol[p1]=$(sol[p1]) spp[p1]=$(spp[p1]))
                                """)
                        return false
                    end
                end
            end
        end
    end
    return true
end

# Push the given solution to a local optimium if needed
function enforce_optimality!(sol::Vector{Int}, interface::Interface)
    np = interface.np

    reqs = interface.reqs
    deps = interface.deps
    pkgs = interface.pkgs
    spp = interface.spp
    pdict = interface.pdict
    pvers = interface.pvers
    vdict = interface.vdict

    # prepare some useful structures
    # pdeps[p0][v0] has all dependencies of package p0 version v0
    pdeps = [Vector{Requires}(uninitialized, spp[p0]-1) for p0 = 1:np]
    # prevdeps[p1][p0][v0] is the VersionSet of package p1 which package p0 version v0
    # depends upon
    prevdeps = [Dict{Int,Dict{Int,VersionSet}}() for p0 = 1:np]

    for (p,d) in deps
        p0 = pdict[p]
        vdict0 = vdict[p0]
        for (vn,a) in d
            v0 = vdict0[vn]
            pdeps[p0][v0] = a.requires
            for (rp, rvs) in a.requires
                p1 = pdict[rp]
                if !haskey(prevdeps[p1], p0)
                    prevdeps[p1][p0] = Dict{Int,VersionSet}()
                end
                prevdeps[p1][p0][v0] = rvs
            end
        end
    end

    restart = true
    while restart
        restart = false
        for p0 = 1:np
            s0 = sol[p0]
            if s0 >= spp[p0] - 1
                # either the package is not installed,
                # or it's already at the maximum version
                continue
            end
            viol = false
            # check if the higher version has a depencency which
            # would be violated by the state of the remaining packages
            for (p,vs) in pdeps[p0][s0+1]
                p1 = pdict[p]
                if sol[p1] == spp[p1]
                    # the dependency is violated because
                    # the other package is not being installed
                    viol = true
                    break
                end
                vn = pvers[p1][sol[p1]]
                if vn ∉ vs
                    # the dependency is violated because
                    # the other package version is invalid
                    viol = true
                    break
                end
            end
            viol && continue

            # check if bumping the version would violate some
            # dependency of another package
            for (p1,d) in prevdeps[p0]
                vs = get(d, sol[p1], nothing)
                vs === nothing && continue
                vn = pvers[p0][s0+1]
                if vn ∉ vs
                    # bumping the version would violate
                    # the dependency
                    viol = true
                    break
                end
            end
            viol && continue
            # So the solution is non-optimal: we bump it manually
            #warn("nonoptimal solution for package $(interface.pkgs[p0]): sol=$s0")
            sol[p0] += 1
            restart = true
        end
    end

    # Finally uninstall unneeded packages:
    # start from the required ones and keep only
    # the packages reachable from them along the graph
    uninst = trues(np)
    staged = Set{String}(keys(reqs))
    seen = copy(staged)

    while !isempty(staged)
        staged_next = Set{String}()
        for p in staged
            p0 = pdict[p]
            uninst[p0] = false
            @assert sol[p0] < spp[p0]
            vn = pvers[p0][sol[p0]]
            a = deps[p][vn]

            # scan dependencies
            for (rp,rvs) in a.requires
                rp0 = pdict[rp]
                @assert sol[rp0] < spp[rp0] && pvers[rp0][sol[rp0]] ∈ rvs
                rp ∈ seen || push!(staged_next, rp)
            end
        end
        union!(seen, staged_next)
        staged = staged_next
    end

    for p0 in find(uninst)
        sol[p0] = spp[p0]
    end

    return
end

end

max_stars_repo_name: JuliaLangEs/julieta

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using DataFrames, DataFramesMeta, CSV
using PGFPlots
using Statistics

# Function for results of Fig. A.6
function getfittedparamvalues_multipleruns(;
        filepath =  joinpath(dirname(pathof(SinergiafMRI_datafit)),"mcmc_rl_fit/projects/fmri/someresults/fits_multipleruns/"),
        filestr = "",
        )

    savepath = makehomesavepath("getfittedparamvalues")
    mkdir(savepath)

    prefix = split(filepath, "/")[end-1]
    runs = getsubfolders(filepath)
    algorithms = getsubfolders(joinpath(filepath, runs[1]))

    for i in algorithms
        dfresults = DataFrame()

        for irun in 1:length(runs)
            if isfile(joinpath(filepath, runs[irun], i, "maxLL_fit/MAP_params.csv"))
                dfalgo = CSV.read(joinpath(filepath, runs[irun], i, "maxLL_fit/MAP_params.csv"), header=false)
                dfresults[:, Symbol("run"*string(irun))] = dfalgo[:, :Column1]
            else
                println("File: $(joinpath(filepath, runs[irun], i, "maxLL_fit/MAP_params.csv")) does not exist!")
            end
        end
        CSV.write(joinpath(savepath, prefix*"_paramvalues_acrossruns_"*i*filestr*".csv"), dfresults,
                delim = " ");
    end
end


text:
getobs(data; obsdim = default_obsdim(data)) =
    __getobs(data, convert(LearnBase.ObsDimension,obsdim))

# These "__" methods serve the sole purpose of avoiding ambiguities
# when a user defines "getobs" for his/her type.
# Basically if obsdim is undefined, we assume there is no
# getobs(::MyType, ::ObsDimension) to care about.
# All others are fed back into "getobs" dispatch
@inline __getobs(data, obsdim::ObsDim.Undefined) = _getobs(data, obsdim)
@inline __getobs(data, obsdim::ObsDimension) = getobs(data, obsdim)
@inline __getobs(data::Tuple, obsdim::Tuple) = getobs(data, obsdim)

getobs!(buffer, data) = getobs(data)
getobs!(buffer, data, idx, obsdim) = getobs(data, idx, obsdim)
getobs!(buffer, data, idx; obsdim = default_obsdim(data)) =
    getobs!(buffer, data, idx, convert(LearnBase.ObsDimension,obsdim))
# NOTE: default to not use buffer since copy! may not be defined
# getobs!(buffer, data) = copy!(buffer, getobs(data))
# getobs!(buffer, data, idx, obsdim) = copy!(buffer, getobs(data, idx, obsdim))

# fallback methods discards unused obsdim
nobs(data, ::ObsDim.Undefined)::Int = nobs(data)
getobs(data, idx, ::ObsDim.Undefined) = getobs(data, idx)

# to accumulate indices as views instead of copies
_view(indices::AbstractRange, i::Int) = indices[i]
_view(indices::AbstractRange, i::AbstractRange) = indices[i]
_view(indices, i::Int) = indices[i] # to throw error in case
_view(indices, i) = view(indices, i)

# can be used to prevent specific data container to be used
# in (specific) functions. For example passing a BatchView to
# oversample or undersample should throw an error
allowcontainer(fun, data) = true

"""
    nobs(data, [obsdim]) -> Int

Return the total number of observations contained in `data`.

The optional parameter `obsdim` can be used to specify which
dimension denotes the observations, if that concept makes sense for
the type of `data`. See `?LearnBase.ObsDim` for more information.
"""
function nobs(data; obsdim = default_obsdim(data))::Int
    nobsdim = convert(LearnBase.ObsDimension,obsdim)
    # make sure we don't bounce between fallback methods
    typeof(nobsdim) <: ObsDim.Undefined && throw(MethodError(nobs, (data,)))
    nobs(data, nobsdim)
end

"""
    getobs(data, [idx], [obsdim])

Return the observation(s) in `data` that correspond to the given
index/indices in `idx`. Note that `idx` can be of type `Int` or
`AbstractVector`.

The returned observation(s) should be in the form intended to be
passed as-is to some learning algorithm. There is no strict
requirement that dictates what form or type that is. We do,
however, expect it to be consistent for `idx` being an integer,
as well as `idx` being an abstract vector, respectively.

If it makes sense for the type of `data`, `obsdim` can be used to
specify which dimension of `data` denotes the observations. It
can be specified in a type-stable manner as a positional argument
(see `?LearnBase.ObsDim`), or more conveniently as a smart
keyword argument.
"""
getobs(data, arg, args...; kw...) = _getobs(data, arg, args...; kw...)
getobs(data, arg) = _getobs(data, arg)

# These "_" methods serve the sole purpose of avoiding ambiguities
# when a user defines "getobs" for his/her type.
# The problem would be getobs(::Any, ::Union{ObsDimension,Tuple}),
# which would likely collide with getobs(::MyType, idx)
@inline function _getobs(data, idx; obsdim = default_obsdim(data))
    nobsdim = convert(LearnBase.ObsDimension,obsdim)
    # make sure we don't bounce between fallback methods
    typeof(nobsdim) <: ObsDim.Undefined && throw(MethodError(getobs, (data,idx)))
    getobs(data, idx, nobsdim)
end

@inline _getobs(data, obsdim::Union{ObsDimension,Tuple}) =
    getobs(data, 1:nobs(data,obsdim), obsdim)

# --------------------------------------------------------------------
# Arrays

nobs(A::AbstractArray, ::ObsDim.Constant{DIM}) where {DIM} = size(A, DIM)::Int
nobs(A::AbstractArray{T,N}, ::ObsDim.Last) where {T,N} = size(A, N)::Int
nobs(A::AbstractArray{T,0}, ::ObsDim.Last) where {T} = 1

getobs(A::Array, ::ObsDimension=default_obsdim(A)) = A
getobs(A::AbstractSparseArray, ::ObsDimension=default_obsdim(A)) = A
getobs(A::SubArray, ::ObsDimension=default_obsdim(A)) = copy(A)
getobs(A::SubArray{T,0}, ::ObsDimension=default_obsdim(A)) where {T} = A[1]
getobs(A::AbstractArray, ::ObsDimension=default_obsdim(A)) = collect(A)

"""
    getobs!(buffer, data, [idx], [obsdim]) -> buffer

Write the observation(s) from `data` that correspond to the given
index/indices in `idx` into `buffer`. Note that `idx` can be of
type `Int` or `AbstractVector`.

Unless explicitly implemented for `data` it defaults to returning
`getobs(data, idx, obsdim)` in which case `buffer` is ignored.

If it makes sense for the type of `data`, `obsdim` can be used to
specify which dimension of `data` denotes the observations. It
can be specified in a type-stable manner as a positional argument
(see `?LearnBase.ObsDim`), or more conveniently as a smart
keyword argument.
"""
getobs!(buffer, A::AbstractSparseArray, idx, obsdim) = getobs(A, idx, obsdim)
getobs!(buffer, A::AbstractSparseArray) = getobs(A)

getobs!(buffer, A::AbstractArray, idx, obsdim) = getobs!(buffer, datasubset(A, idx, obsdim))
getobs!(buffer, A::AbstractArray) = copyto!(buffer, A)

getobs!(buffer, A::SubArray{T,0}) where {T} = A[1]

# catch the undefined setting for consistency.
# should never happen by accident
getobs(A::AbstractArray, idx, obsdim::ObsDim.Undefined) =
    throw(MethodError(getobs, (A, idx, obsdim)))

@generated function getobs(A::AbstractArray{T,N}, idx, obsdim::ObsDimension) where {T,N}
    if N <= 1 && idx <: Integer
        :(A[idx])
    elseif N == 0
        throw(ArgumentError("attempting to index 0-dimensional array using indices of type $idx"))
    elseif obsdim <: ObsDim.First
        :(getindex(A, idx, $(fill(:(:),N-1)...)))
    elseif obsdim <: ObsDim.Last || (obsdim <: ObsDim.Constant && obsdim.parameters[1] == N)
        :(getindex(A, $(fill(:(:),N-1)...), idx))
    else # obsdim <: ObsDim.Constant
        DIM = obsdim.parameters[1]
        DIM > N && throw(DimensionMismatch("the given obsdim=$DIM is greater than the number of available dimensions N=$N"))
        :(getindex(A, $(fill(:(:),DIM-1)...), idx, $(fill(:(:),N-DIM)...)))
    end
end

# --------------------------------------------------------------------
# Tuples

_check_nobs_error() = throw(DimensionMismatch("all data container must have the same number of observations"))

function _check_nobs(tup::Tuple)
    length(tup) == 0 && return
    n1 = nobs(tup[1])
    for i=2:length(tup)
        nobs(tup[i]) != n1 && _check_nobs_error()
    end
end

function _check_nobs(tup::Tuple, obsdim::ObsDimension)
    length(tup) == 0 && return
    n1 = nobs(tup[1], obsdim)
    for i=2:length(tup)
        nobs(tup[i], obsdim) != n1 && _check_nobs_error()
    end
end

function _check_nobs(tup::Tuple, obsdims::Tuple)
    length(tup) == 0 && return
    length(tup) == length(obsdims) || throw(DimensionMismatch("number of elements in obsdim doesn't match data"))
    all(map(x-> typeof(x) <: Union{ObsDimension,Tuple}, obsdims)) || throw(MethodError(_check_nobs, (tup, obsdims)))
    n1 = nobs(tup[1], obsdims[1])
    for i=2:length(tup)
        nobs(tup[i], obsdims[i]) != n1 && _check_nobs_error()
    end
end

function nobs(tup::Tuple, ::ObsDim.Undefined = ObsDim.Undefined())::Int
    _check_nobs(tup)
    length(tup) == 0 ? 0 : nobs(tup[1])
end

function nobs(tup::Tuple, obsdim::ObsDimension)::Int
    _check_nobs(tup, obsdim)
    length(tup) == 0 ? 0 : nobs(tup[1], obsdim)
end

function nobs(tup::Tuple, obsdims::Tuple)::Int
    _check_nobs(tup, obsdims)
    length(tup) == 0 ? 0 : nobs(tup[1], obsdims[1])
end

function getobs(tup::Tuple, obsdim::Tuple)
    _check_nobs(tup)
    map(getobs, tup, obsdim)
end

function getobs(tup::Tuple, obsdim::ObsDimension)
    _check_nobs(tup, obsdim)
    map(data -> getobs(data, obsdim), tup)
end

function getobs(tup::Tuple, indices)
    _check_nobs(tup)
    map(data -> getobs(data, indices), tup)
end

function getobs(tup::Tuple, indices, obsdim::ObsDimension)
    _check_nobs(tup, obsdim)
    map(data -> getobs(data, indices, obsdim), tup)
end

@generated function getobs(tup::Tuple, indices, obsdims::Tuple)
    N = length(obsdims.types)
    quote
        _check_nobs(tup, obsdims)
        # This line generates a tuple of N elements:
        # (getobs(tup[1], indices, obsdims[1]), getobs(tup[2], indi...
        $(Expr(:tuple, (:(getobs(tup[$i], indices, obsdims[$i])) for i in 1:N)...))
    end
end

_getobs_error() = throw(DimensionMismatch("the first argument (tuple with the buffers) must have the same length as the second argument (tuple with the data container)"))

@generated function getobs!(buffer::Tuple, tup::Tuple)
    N = length(buffer.types)
    N == length(tup.types) || _getobs_error()
    quote
        # _check_nobs(tup) # don't check because of single obs
        $(Expr(:tuple, (:(getobs!(buffer[$i],tup[$i])) for i in 1:N)...))
    end
end

@generated function getobs!(buffer::Tuple, tup::Tuple, indices, obsdim)
    N = length(buffer.types)
    N == length(tup.types) || _getobs_error()
    expr = if obsdim <: ObsDimension
        Expr(:tuple, (:(getobs!(buffer[$i], tup[$i], indices, obsdim)) for i in 1:N)...)
    else
        Expr(:tuple, (:(getobs!(buffer[$i], tup[$i], indices, obsdim[$i])) for i in 1:N)...)
    end
    quote
        # _check_nobs(tup, obsdim) # don't check because of single obs
        $expr
    end
end


text:
using FunctionalCollections
using Random: MersenneTwister
using UUIDs

const newid = let rng = MersenneTwister()
    function(prefix)
        string(prefix, '-', uuid4(rng))
    end
end

_pvec(x::PersistentVector) = x
_pvec(x::AbstractArray) = pvec(x)

# b can be array of pairs / kwargs etc.
function recmerge!(a, b, f=recmerge!)
    for (k, v) in b
        if isa(v, AbstractDict) && haskey(a, k) && isa(a[k], AbstractDict)
            a[k] = f(a[k], v)
        else
            a[k] = v
        end
    end
    a
end

recmerge(a, b) = recmerge!(Dict{Any, Any}(a), b, recmerge)

function kebab2camel(str)
    strs = split(str, "-")
    titlefrom2 = titlecase.(strs[2:end])
    join([strs[1]; titlefrom2], "")
end

function camel2kebab(str)
    lowstr = lowercase(str)
    lowercaseidxs = Vector{Char}(lowstr) .== Vector{Char}(str)
    sb = IOBuffer()
    for (i,c) in enumerate(lowercase(str))
        !lowercaseidxs[i] && print(sb, '-')
        print(sb, c)
    end
    String(take!(sb))
end


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    html_div("hello", id="input"),
    html_div(id="side-effect"),
    html_div("output", id="output")
end

callback!(
    ClientsideFunction("clientside", "side_effect_and_return_a_promise"),
    app, Output("output","children"), Input("input","children")
    )

run_server(app)


text:
## @lc code=start
using LeetCode
using DataStructures
function number_of_good_subsets(nums::Vector{Int})::Int
    table, MOD = counter(nums), 10 ^ 9 + 7
    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
    good_part = [i for i in 2:30 if haskey(table, i) && all(i % p ^ 2 != 0 for p in primes)]
    good_nums = DefaultDict{Int, Int}(0)
    good_nums[1] = (1 << table[1]) % MOD
    for num in good_part, x in collect(keys(good_nums))
        if gcd(num, x) == 1
            good_nums[num * x] += table[num] * good_nums[x]
            good_nums[num * x] %= MOD
        end
    end
    return (sum(values(good_nums)) - good_nums[1]) % MOD
end
## @lc code=end

## @lc test=start
@testset "1994.the-number-of-good-subsets.jl" begin
    @test number_of_good_subsets([1,2,3,4]) == 6
    @test number_of_good_subsets([4,2,3,15]) == 5
    @test number_of_good_subsets([5,10,1,26,24,21,24,23,11,12,27,4,17,16,2,6,1,1,6,8,13,30,24,20,2,19]) == 5368
end
## @lc test=end

## @lc info=start
# link: [solution to 1994]()
## @lc info=end


text:
module TestZygote

include("preamble.jl")
using Zygote

include("../examples/zygote.jl")
@test g isa Matrix
@test size(g) == (d, d)

end  # module


text:
module Approximators

using Random
using StatsBase: mean, pweights

using Ramnet.Models
using Ramnet.Encoders

include("PolicyApproximators/BinaryActionPolicy.jl")
export BinaryActionPolicy

include("PolicyApproximators/RegularizedBinaryPolicy.jl")
export RegularizedBinaryPolicy

include("PolicyApproximators/AdaptiveBinaryPolicy.jl")
export AdaptiveBinaryPolicy

include("PolicyApproximators/DiscreteActionPolicy.jl")
export DiscreteActionPolicy

include("PolicyApproximators/ContinousActionPolicy.jl")
export ContinousActionPolicy

include("PolicyApproximators/RegularizedContinousPolicy.jl")
export RegularizedContinousPolicy

include("PolicyApproximators/AltContinousActionPolicy.jl")
export AltContinousActionPolicy

include("PolicyApproximators/AdaptiveContinousPolicy.jl")
export AdaptiveContinousPolicy

include("ValueApproximators/VDiscriminator.jl")
export VDiscriminator

include("ValueApproximators/QDiscriminator.jl")
export QDiscriminator

struct DiscreteQDiscriminator{O <: AbstractVector,A <: Real,E <: AbstractEncoder}
    action_range::UnitRange{A}
    discount::Float64
    encoder::E
    Q::Dict{A,RegressionDiscriminator}

    # TODO: It's annoying to have to specify the input size. It requires knowledge
    #       of the encoder's resolution and the length of the state of the environment.
    #       Maybe it would be better to have a type that specifies the characteristics
    #       of the environment, such as action and observation spaces.
    function DiscreteQDiscriminator{O,A,E}(action_range::UnitRange{A}, size::Int, n::Int, discount::Float64, encoder::E) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
        n < 1 && throw(DomainError(n, "Tuple size must be at least 1"))
        size < n && throw(DomainError(size, "Input size may not be smaller then tuple size"))
        !(0 ≤ discount ≤ 1) && throw(DomainError(discount, "Discount must lie in the [0,1] interval"))

        Q = Dict(a => RegressionDiscriminator(size, n; γ=discount) for a in action_range)

        new(action_range, discount, encoder, Q)
    end
end

# TODO: Properly check if observation and action are valid
function (M::DiscreteQDiscriminator{O,A,E})(observation::O, action::A) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
    return predict(M.Q[action], encode(M.encoder, observation))
end

# TODO: Properly check if observation is valid
function q_values!(M::DiscreteQDiscriminator{O,A,E}, observation::O, q_values::Vector{Float64}) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
    length(q_values) != length(M.action_range) && throw(
      DimensionMismatch("Expected vector of length $(length(M.action_range)), got $(length(q_values))"))

    q_max = typemin(Float64)
    q = 0.0
    encoded_obs = encode(M.encoder, observation) # TODO: Don't have to create a bunch of temp arrays here
    for (i, action) in Iterators.zip(eachindex(q_values), M.action_range)
        q = predict(M.QQ[action], encoded_obs)
        dest[i] = q
        q_max = q > q_max ? q : q_max
    end

    return q_max
end

function (M::DiscreteQDiscriminator{O,A,E})(observation::O) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
    values = similar(M.action_range, Float64)
    q_values!(M, observation, values)

return values
end

function expected_q_value(M::DiscreteQDiscriminator{O,A,E}, observation::O; ϵ::Float64=0.0) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
    values = similar(M.action_range, Float64)
    q_max = q_values!(M, observation, values)

    return mean(values, pweights((values .== q_max) .+ ϵ))
end

# function select_greedy_action(M::DiscreteQDiscriminator{O,A,E}, observation, rng::MersenneTwister; ϵ::Float64=0.0 ) where {A <: AbstractAgent}
#     # TODO: I could avoid creating these temporary arrays here
#     values = similar(M.action_range, Float64)
#     q_max = q_values!(M, observation, values)

#     if rand(agent.rng) < ϵ # Random action taken with ϵ probability
#         agent.action = rand(agent.rng, agent.actions)
#     else
#       # Selecting action that maximizes Q̂, with ties broken randomly
#         agent.action = sample(agent.rng, agent.actions, pweights(values .== q_max))
#     end

#     return agent.action
# end

end


text:
using Mimi
using OptiMimi
using Test

# Easier basic model to optimize:
# Each period, you have 1 unit to divide between consuming and staving off disaster.
# If disaster happens, you get 0 consumption for that period.
# Utility is defined as sqrt(c).  The probability of disaster is 1-c in each period.
# Assume 5% discounting and 10 periods.

## Analytical solution in 1 period

## max_c (1 - c) u(c) => c = 1/3
gamma = exp(-.05) # 5% discount rate

# OptiMimi solution

# Create simple model
@defcomp Bellmano begin
    # The x-value to evaluate the quadratic
    consumption = Parameter(index=[time])

    # The y-value of the quadratic at the x-value
    utility = Variable(index=[time])

    function run_timestep(p, v, d, t)
        #v.utility[tt] = sqrt(p.consumption[tt]) * (1 - p.consumption[tt])

        if rand() < p.consumption[tt]
            v.utility[tt] = 0
        else
            v.utility[tt] = sqrt(max(0, p.consumption[tt]))
        end
    end
end

m = Model()
set_dimension!(m, :time, collect(1:10))

bellmano = add_comp!(m, Bellmano)
bellmano[:consumption] = repmat([.25], 10)

run(m)

import OptiMimi.uncertainproblem
using DataFrames

df = DataFrame(algorithm=[], mcperlife=[], cons1=[], cons2=[], cons3=[], cons4=[], cons5=[], cons6=[], cons7=[], cons8=[], cons9=[], cons10=[], cons1serr=[], cons2serr=[], cons3serr=[], cons4serr=[], cons5serr=[], cons6serr=[], cons7serr=[], cons8serr=[], cons9serr=[], cons10serr=[], time=[])
push!(df, [:exact; Inf; repmat([1/3], 10); repmat([0], 10); 0])

for samplemult in [1, 2, 4]
    for mcperlife in [5, 20, 50]
        println(mcperlife)
        prob = uncertainproblem(m, [:Bellmano], [:consumption], [0.], [1.], m -> sum(m[:Bellmano, :utility] .* exp(-(0:9) * .05)), (model) -> nothing)

        tic()
        soln = solution(prob, () -> repmat([.25], 10), :social, mcperlife, samplemult)
        time = toc()
        push!(df, [:social; mcperlife; soln.xmean; soln.xserr; time])

        println(df)

        tic()
        soln = solution(prob, () -> repmat([.25], 10), :biological, mcperlife, 120000*samplemult)
        time = toc()
        push!(df, [:biological; mcperlife; soln.xmean; soln.xserr; time])

        println(df)

        tic()
        soln = solution(prob, () -> repmat([.25], 10), :sampled, mcperlife, 300*samplemult)
        time = toc()
        push!(df, [:sampled; mcperlife; soln.xmean; soln.xserr; time])

        println(df)
    end

    writetable("attempts-compare2.csv", df)
end


text:
module TKEBudgetTerms

export TurbulentKineticEnergy, KineticEnergy
export IsotropicViscousDissipationRate, IsotropicPseudoViscousDissipationRate
export AnisotropicPseudoViscousDissipationRate
export XPressureRedistribution, YPressureRedistribution, ZPressureRedistribution
export XShearProduction, YShearProduction, ZShearProduction

using Oceananigans
using Oceananigans.Operators
using Oceananigans.AbstractOperations
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.Grids: Center, Face
using Oceananigans.Fields: ZeroField

# Some useful operators
@inline ψ²(i, j, k, grid, ψ) = @inbounds ψ[i, j, k]^2

@inline ψ′²(i, j, k, grid, ψ, Ψ) = @inbounds (ψ[i, j, k] - Ψ[i, j, k])^2
@inline ψ′²(i, j, k, grid, ψ, Ψ::Number) = @inbounds (ψ[i, j, k] - Ψ)^2

@inline fψ²(i, j, k, grid, f, ψ) = @inbounds f(i, j, k, grid, ψ)^2

@inline fψ_plus_gφ²(i, j, k, grid, f, ψ, g, φ) = @inbounds (f(i, j, k, grid, ψ) + g(i, j, k, grid, φ))^2

#++++ Turbulent kinetic energy
@inline function turbulent_kinetic_energy_ccc(i, j, k, grid, u, v, w, U, V, W)
    return (ℑxᶜᵃᵃ(i, j, k, grid, ψ′², u, U) +
            ℑyᵃᶜᵃ(i, j, k, grid, ψ′², v, V) +
            ℑzᵃᵃᶜ(i, j, k, grid, ψ′², w, W)) / 2
end

function TurbulentKineticEnergy(model, u, v, w;
                                U = 0,
                                V = 0,
                                W = 0,
                                location = (Center, Center, Center))

    if location == (Center, Center, Center)
        return KernelFunctionOperation{Center, Center, Center}(turbulent_kinetic_energy_ccc, model.grid,
                                       computed_dependencies=(u, v, w, U, V, W))
    else
        error("TurbulentKineticEnergy only supports location = (Center, Center, Center) for now.")
    end
end

KineticEnergy(model, u, v, w; location = (Center, Center, Center), kwargs...) =
    TurbulentKineticEnergy(model, u, v, w; location, kwargs...)

TurbulentKineticEnergy(model; kwargs...) = TurbulentKineticEnergy(model, model.velocities...; kwargs...)
KineticEnergy(model; kwargs...) = KineticEnergy(model, model.velocities...; kwargs...)
#------

#++++ Energy dissipation rate for a fluid with isotropic viscosity
function isotropic_viscous_dissipation_rate_les_ccc(i, j, k, grid, u, v, w, ν)

    Σˣˣ² = ∂xᶜᵃᵃ(i, j, k, grid, u)^2
    Σʸʸ² = ∂yᵃᶜᵃ(i, j, k, grid, v)^2
    Σᶻᶻ² = ∂zᵃᵃᶜ(i, j, k, grid, w)^2

    Σˣʸ² = ℑxyᶜᶜᵃ(i, j, k, grid, fψ_plus_gφ², ∂yᵃᶠᵃ, u, ∂xᶠᵃᵃ, v) / 4
    Σˣᶻ² = ℑxzᶜᵃᶜ(i, j, k, grid, fψ_plus_gφ², ∂zᵃᵃᶠ, u, ∂xᶠᵃᵃ, w) / 4
    Σʸᶻ² = ℑyzᵃᶜᶜ(i, j, k, grid, fψ_plus_gφ², ∂zᵃᵃᶠ, v, ∂yᵃᶠᵃ, w) / 4

    return ν[i, j, k] * 2 * (Σˣˣ² + Σʸʸ² + Σᶻᶻ² + 2 * (Σˣʸ² + Σˣᶻ² + Σʸᶻ²))
end

function isotropic_viscous_dissipation_rate_dns_ccc(i, j, k, grid, u, v, w, ν)

    Σˣˣ² = ∂xᶜᵃᵃ(i, j, k, grid, u)^2
    Σʸʸ² = ∂yᵃᶜᵃ(i, j, k, grid, v)^2
    Σᶻᶻ² = ∂zᵃᵃᶜ(i, j, k, grid, w)^2

    Σˣʸ² = ℑxyᶜᶜᵃ(i, j, k, grid, fψ_plus_gφ², ∂yᵃᶠᵃ, u, ∂xᶠᵃᵃ, v) / 4
    Σˣᶻ² = ℑxzᶜᵃᶜ(i, j, k, grid, fψ_plus_gφ², ∂zᵃᵃᶠ, u, ∂xᶠᵃᵃ, w) / 4
    Σʸᶻ² = ℑyzᵃᶜᶜ(i, j, k, grid, fψ_plus_gφ², ∂zᵃᵃᶠ, v, ∂yᵃᶠᵃ, w) / 4

    return ν * 2 * (Σˣˣ² + Σʸʸ² + Σᶻᶻ² + 2 * (Σˣʸ² + Σˣᶻ² + Σʸᶻ²))
end

"""
Calculates the Viscous Dissipation Rate for a fluid with an isotropic turbulence closure (i.e., a
turbulence closure where ν (eddy or not) is the same for all directions.
"""
function IsotropicViscousDissipationRate(model; U=nothing, V=nothing, W=nothing,
                                         location = (Center, Center, Center))
    if location != (Center, Center, Center)
        error("IsotropicViscousDissipationRate only supports location = (Center, Center, Center) for now.")
    end

    u, v, w = model.velocities
    if U != nothing
        u -= U
    end

    if V != nothing
        v -= V
    end

    if W != nothing
        w -= W
    end

    if model.closure isa Oceananigans.TurbulenceClosures.AbstractEddyViscosityClosure
        νₑ = model.diffusivity_fields.νₑ
        return KernelFunctionOperation{Center, Center, Center}(isotropic_viscous_dissipation_rate_les_ccc, model.grid;
                                                               computed_dependencies=(u, v, w, νₑ))

    elseif model.closure == nothing
        error("Trying to calculate a TKE viscous dissipation rate with `model.closure==nothing`.")

    else
        ν = model.closure.ν
        return KernelFunctionOperation{Center, Center, Center}(isotropic_viscous_dissipation_rate_dns_ccc, model.grid;
                                                               computed_dependencies=(u, v, w), parameters=ν)
    end

end

function isotropic_pseudo_viscous_dissipation_rate_les_ccc(i, j, k, grid, u, v, w, ν)
    ddx² = ∂xᶜᵃᵃ(i, j, k, grid, ψ², u) + ℑxyᶜᶜᵃ(i, j, k, grid, fψ², ∂xᶠᵃᵃ, v) + ℑxzᶜᵃᶜ(i, j, k, grid, fψ², ∂xᶠᵃᵃ, w)
    ddy² = ℑxyᶜᶜᵃ(i, j, k, grid, fψ², ∂yᵃᶠᵃ, u) + ∂yᵃᶜᵃ(i, j, k, grid, ψ², v) + ℑyzᵃᶜᶜ(i, j, k, grid, fψ², ∂yᵃᶠᵃ, w)
    ddz² = ℑxzᶜᵃᶜ(i, j, k, grid, fψ², ∂zᵃᵃᶠ, u) + ℑyzᵃᶜᶜ(i, j, k, grid, fψ², ∂zᵃᵃᶠ, v) + ∂zᵃᵃᶜ(i, j, k, grid, ψ², w)
    return ν[i,j,k] * (ddx² + ddy² + ddz²)
end

function isotropic_pseudo_viscous_dissipation_rate_dns_ccc(i, j, k, grid, u, v, w, ν)
    ddx² = ∂xᶜᵃᵃ(i, j, k, grid, ψ², u) + ℑxyᶜᶜᵃ(i, j, k, grid, fψ², ∂xᶠᵃᵃ, v) + ℑxzᶜᵃᶜ(i, j, k, grid, fψ², ∂xᶠᵃᵃ, w)
    ddy² = ℑxyᶜᶜᵃ(i, j, k, grid, fψ², ∂yᵃᶠᵃ, u) + ∂yᵃᶜᵃ(i, j, k, grid, ψ², v) + ℑyzᵃᶜᶜ(i, j, k, grid, fψ², ∂yᵃᶠᵃ, w)
    ddz² = ℑxzᶜᵃᶜ(i, j, k, grid, fψ², ∂zᵃᵃᶠ, u) + ℑyzᵃᶜᶜ(i, j, k, grid, fψ², ∂zᵃᵃᶠ, v) + ∂zᵃᵃᶜ(i, j, k, grid, ψ², w)
    return ν * (ddx² + ddy² + ddz²)
end

function IsotropicPseudoViscousDissipationRate(model; U=nothing, V=nothing, W=nothing,
                                               location = (Center, Center, Center))
    if location != (Center, Center, Center)
        error("IsotropicPseudoViscousDissipationRate only supports location = (Center, Center, Center) for now.")
    end

    u, v, w = model.velocities
    if U != nothing
        u -= U
    end

    if V != nothing
        v -= V
    end

    if W != nothing
        w -= W
    end

    if model.closure isa Oceananigans.TurbulenceClosures.AbstractEddyViscosityClosure
        νₑ = model.diffusivity_fields.νₑ
        return KernelFunctionOperation{Center, Center, Center}(isotropic_pseudo_viscous_dissipation_rate_les_ccc, model.grid;
                                       computed_dependencies=(u, v, w, νₑ))
    elseif model.closure == nothing
        error("Trying to calculate a TKE pseudo viscous dissipation rate with `model.closure==nothing`.")

    else
        ν = model.closure.ν
        return KernelFunctionOperation{Center, Center, Center}(isotropic_pseudo_viscous_dissipation_rate_dns_ccc, model.grid;
                                                               computed_dependencies=(u, v, w), parameters=ν)
    end

end
#------

#+++++ Energy dissipation rate for a fluid with constant anisotropic viscosity (closure)
function anisotropic_pseudo_viscous_dissipation_rate_ccc(i, j, k, grid, u, v, w, params)

    ddx² = ∂xᶜᵃᵃ(i, j, k, grid, ψ², u) + ℑxyᶜᶜᵃ(i, j, k, grid, fψ², ∂xᶠᵃᵃ, v) + ℑxzᶜᵃᶜ(i, j, k, grid, fψ², ∂xᶠᵃᵃ, w)
    ddy² = ℑxyᶜᶜᵃ(i, j, k, grid, fψ², ∂yᵃᶠᵃ, u) + ∂yᵃᶜᵃ(i, j, k, grid, ψ², v) + ℑyzᵃᶜᶜ(i, j, k, grid, fψ², ∂yᵃᶠᵃ, w)
    ddz² = ℑxzᶜᵃᶜ(i, j, k, grid, fψ², ∂zᵃᵃᶠ, u) + ℑyzᵃᶜᶜ(i, j, k, grid, fψ², ∂zᵃᵃᶠ, v) + ∂zᵃᵃᶜ(i, j, k, grid, ψ², w)

    return params.νx*ddx² + params.νy*ddy² + params.νz*ddz²
end

function AnisotropicPseudoViscousDissipationRate(model; U=nothing, V=nothing, W=nothing,
                                                 location = (Center, Center, Center))

    if location != (Center, Center, Center)
        error("AnisotropicPseudoViscousDissipationRate only supports location = (Center, Center, Center) for now.")
    end

    u, v, w = model.velocities
    if U != nothing
        u -= U
    end

    if V != nothing
        v -= V
    end

    if W != nothing
        w -= W
    end

    νx = model.closure.νx
    νy = model.closure.νy
    νz = model.closure.νz
    return KernelFunctionOperation{Center, Center, Center}(anisotropic_pseudo_viscous_dissipation_rate_ccc, model.grid;
                                                           computed_dependencies=(u, v, w),
                                                           parameters=(νx=νx, νy=νy, νz=νz,))
end
#-----

#++++ Pressure redistribution terms
function XPressureRedistribution(model)
    u, v, w = model.velocities
    p = sum(model.pressures)
    return ∂x(u*p) # p is the total kinematic pressure (there's no need to ρ₀)
end

function YPressureRedistribution(model)
    u, v, w = model.velocities
    p = sum(model.pressures)
    return ∂y(v*p) # p is the total kinematic pressure (there's no need to ρ₀)
end

function ZPressureRedistribution(model)
    u, v, w = model.velocities
    p = sum(model.pressures)
    return ∂z(w*p) # p is the total kinematic pressure (there's no need to ρ₀)
end
#----

#++++ Shear production terms
function shear_production_x_ccc(i, j, k, grid, u, v, w, U, V, W)
    u_int = ℑxᶜᵃᵃ(i, j, k, grid, u) # F, C, C  → C, C, C

    ∂xU = ∂xᶜᵃᵃ(i, j, k, grid, U) # F, C, C  → C, C, C
    uu = ℑxᶜᵃᵃ(i, j, k, grid, ψ², u)
    uu∂xU = uu * ∂xU

    ∂xV = ℑxyᶜᶜᵃ(i, j, k, grid, ∂xᶠᵃᵃ, V) # C, F, C  → F, F, C  → C, C, C
    vu = ℑyᵃᶜᵃ(i, j, k, grid, v) * u_int
    vu∂xV = vu * ∂xV

    ∂xW = ℑxzᶜᵃᶜ(i, j, k, grid, ∂xᶠᵃᵃ, W) # C, C, F  → F, C, F  → C, C, C
    wu = ℑzᵃᵃᶜ(i, j, k, grid, w) * u_int
    wu∂xW = wu * ∂xW

    return -(uu∂xU + vu∂xV + wu∂xW)
end

function XShearProduction(model, u, v, w, U, V, W; location = (Center, Center, Center))
    if location == (Center, Center, Center)
        return KernelFunctionOperation{Center, Center, Center}(shear_production_x_ccc, model.grid;
                                       computed_dependencies=(u, v, w, U, V, W))
    else
        error("XShearProduction only supports location = (Center, Center, Center) for now.")
    end
end

function shear_production_y_ccc(i, j, k, grid, u, v, w, U, V, W)
    v_int = ℑyᵃᶜᵃ(i, j, k, grid, v) # C, F, C  → C, C, C

    ∂yU = ℑxyᶜᶜᵃ(i, j, k, grid, ∂yᵃᶠᵃ, U) # F, C, C  → F, F, C  → C, C, C
    uv = ℑxᶜᵃᵃ(i, j, k, grid, u) * v_int
    uv∂yU = uv * ∂yU

    ∂yV = ∂yᵃᶜᵃ(i, j, k, grid, V)
    vv = ℑyᵃᶜᵃ(i, j, k, grid, ψ², v) # C, F, C  → C, C, C
    vv∂yV = vv * ∂yV

    ∂yW = ℑyzᵃᶜᶜ(i, j, k, grid, ∂yᵃᶠᵃ, W) # C, C, F  → C, F, F  → C, C, C
    wv = ℑzᵃᵃᶜ(i, j, k, grid, w) * v_int
    wv∂yW = wv * ∂yW

    return -(uv∂yU + vv∂yV + wv∂yW)
end

function YShearProduction(model, u, v, w, U, V, W; location = (Center, Center, Center))
    if location == (Center, Center, Center)
        return KernelFunctionOperation{Center, Center, Center}(shear_production_y_ccc, model.grid;
                                       computed_dependencies=(u, v, w, U, V, W))
    else
        error("YShearProduction only supports location = (Center, Center, Center) for now.")
    end
end

function shear_production_z_ccc(i, j, k, grid, u, v, w, U, V, W)
    w_int = ℑzᵃᵃᶜ(i, j, k, grid, w) # C, C, F  → C, C, C

    ∂zU = ℑxzᶜᵃᶜ(i, j, k, grid, ∂zᵃᵃᶠ, U) # F, C, C  → F, C, F  → C, C, C
    uw = ℑxᶜᵃᵃ(i, j, k, grid, u) * w_int
    uw∂zU = uw * ∂zU

    ∂zV = ℑyzᵃᶜᶜ(i, j, k, grid, ∂zᵃᵃᶠ, V) # C, F, C  → C, F, F  → C, C, C
    vw = ℑyᵃᶜᵃ(i, j, k, grid, v) * w_int
    vw∂zV = vw * ∂zV

    ∂zW = ∂zᵃᵃᶜ(i, j, k, grid, W)
    ww = ℑzᵃᵃᶜ(i, j, k, grid, ψ², w) # C, C, F  → C, C, C
    ww∂zW = ww * ∂zW

    return - (uw∂zU + vw∂zV + ww∂zW)
end

function ZShearProduction(model, u, v, w, U, V, W; location = (Center, Center, Center))
    if location == (Center, Center, Center)
        return KernelFunctionOperation{Center, Center, Center}(shear_production_z_ccc, model.grid;
                                       computed_dependencies=(u, v, w, U, V, W))
    else
        error("ZShearProduction only supports location = (Center, Center, Center) for now.")
    end
end

ZShearProduction(model; U=ZeroField(), V=ZeroField(), W=ZeroField(), kwargs...) =
    ZShearProduction(model, model.velocities..., U, V, W; kwargs...)
#----

end # module


text:
module TestVsSVector
include("preamble.jl")

using Bifurcations: resolved_points
using Bifurcations.Continuations: find_errors, print_errors
using Bifurcations: examples
using Bifurcations.Examples.Calcium
using Bifurcations.Examples.PredatorPrey

use_array(ex) =
    ex.make_prob(
        f = generalized_f(ex),
        u0 = Array(ex.u0),
    )

EXAMPLES = [
    # name => (
    #     prob1_svec,   # SVECtor
    #     prob1_stda,   # STanDard Array
    #     sweep_ids,    # sweeps to be tested
    # )
    :Calcium => (
        Calcium.prob,
        use_array(Calcium),
        :,
        ),
    # PredatorPrey* tests fail at sweeps[4].u[8].  Probably this kind
    # of test is not a good idea after all.  But let's keep it as-is
    # for now to "freeze" the implementation.
    # TODO: Find out how PredatorPrey* tests diverge.
    :PredatorPreyImmutableState => (
        PredatorPrey.prob,
        use_array(PredatorPrey),
        1:3,
    ),
    :PredatorPreyMutableState => (
        PredatorPrey.prob,
        PredatorPrey.make_prob(u0=Array(PredatorPrey.u0)),
        1:3,
    ),
]

@testset "example $name" for (name, probs) in EXAMPLES
    prob1_svec, prob1_stda, sweep_ids = probs
    solver1_svec = init(prob1_svec)
    solver1_stda = init(prob1_stda)

    solve!(solver1_svec)
    solve!(solver1_stda)
    sweeps_svec = solver1_svec.super.sol.sweeps
    sweeps_stda = solver1_stda.super.sol.sweeps

    @test length(sweeps_svec) == length(sweeps_stda)
    for (i, (sweep_svec, sweep_stda)) in enumerate(zip(sweeps_svec[sweep_ids],
                                                       sweeps_stda[sweep_ids]))
        @test length(sweep_svec.u) == length(sweep_stda.u)
        for (j, (u_svec, u_stda)) in enumerate(zip(sweep_svec.u,
                                                   sweep_stda.u))
            # Manually abort on the first test failure:
            if ! (u_svec ≈ u_stda)
                println("sweeps[$i].u[$j]")
                @test u_svec ≈ u_stda
                break
            else
                @test u_svec ≈ u_stda
            end
        end
    end
end

end  # module


text:
using CSV
using DataFrames
using Plots
using LinearAlgebra
using Statistics
using ProgressMeter
using DataStructures
using EcologicalNetworks


text:
"""
    ConstantStepsize <: Stepsize

A functor that always returns a fixed step size.

# Fields
* `length` – constant value for the step size.

# Constructors

    ConstantStepsize(s)

initialize the stepsize to a constant `s` (deprecated)

    ConstantStepsize(M::AbstractManifold=DefaultManifold(2); stepsize=injectivity_radius(M)/2)

initialize the stepsize to a constant `stepsize`, which by default is half the injectivity
radius, unless the radius is infinity, then the default step size is `1`.
"""
mutable struct ConstantStepsize{T} <: Stepsize
    length::T
end
@deprecate ConstantStepsize(s::Real) ConstantStepsize(; stepsize=s)
function ConstantStepsize(
    M::AbstractManifold=DefaultManifold(2);
    stepsize=isinf(injectivity_radius(M)) ? 1.0 : injectivity_radius(M) / 2,
)
    return ConstantStepsize{typeof(stepsize)}(stepsize)
end
function (s::ConstantStepsize)(
    p::P, o::O, i::Int, args...; kwargs...
) where {P<:Problem,O<:Options}
    return s.length
end
get_initial_stepsize(s::ConstantStepsize) = s.length

@doc raw"""
    DecreasingStepsize()

A functor that represents several decreasing step sizes

# Fields
* `length` – (`1`) the initial step size ``l``.
* `factor` – (`1`) a value ``f`` to multiply the initial step size with every iteration
* `subtrahend` – (`0`) a value ``a`` that is subtracted every iteration
* `exponent` – (`1`) a value ``e`` the current iteration numbers ``e``th exponential
  is taken of
* `shift` – (`0`) shift the denominator iterator ``i`` by ``s```.

In total the complete formulae reads for the ``i``th iterate as

````math
s_i = \frac{(l - i a)f^i}{(i+s)^e}
````

and hence the default simplifies to just ``s_i = \frac{l}{i}``

# Constructor

    DecreasingStepsize(l=1,f=1,a=0,e=1,s=0)

Alternatively one can also use the following keyword.

    DecreasingStepsize(
        M::AbstractManifold=DefaultManifold(3);
        length=injectivity_radius(M)/2, multiplier=1.0, subtrahend=0.0, exponent=1.0, shift=0)

initialiszes all fields above, where none of them is mandatory and the length is set to
half and to $1$ if the injectivity radius is infinite.
"""
mutable struct DecreasingStepsize <: Stepsize
    length::Float64
    factor::Float64
    subtrahend::Float64
    exponent::Float64
    shift::Int
    function DecreasingStepsize(l::Real, f::Real=1.0, a::Real=0.0, e::Real=1.0, s::Int=0)
        return new(l, f, a, e, s)
    end
end
function DecreasingStepsize(
    M::AbstractManifold=DefaultManifold(3);
    length=isinf(manifold_dimension(M)) ? 1.0 : manifold_dimension(M) / 2,
    factor=1.0,
    subtrahend=0.0,
    exponent=1.0,
    shift=0,
)
    return DecreasingStepsize(length, factor, subtrahend, exponent, shift)
end
function (s::DecreasingStepsize)(
    ::P, ::O, i::Int, args...; kwargs...
) where {P<:Problem,O<:Options}
    return (s.length - i * s.subtrahend) * (s.factor^i) / ((i + s.shift)^(s.exponent))
end
get_initial_stepsize(s::DecreasingStepsize) = s.length

"""
    Linesearch <: Stepsize

An abstract functor to represent line search type step size deteminations, see
[`Stepsize`](@ref) for details. One example is the [`ArmijoLinesearch`](@ref)
functor.

Compared to simple step sizes, the linesearch functors provide an interface of
the form `(p,o,i,η) -> s` with an additional (but optional) fourth parameter to
provide a search direction; this should default to something reasonable, e.g.
the negative gradient.
"""
abstract type Linesearch <: Stepsize end

@doc raw"""
    ArmijoLinesearch <: Linesearch

A functor representing Armijo line search including the last runs state, i.e. a
last step size.

# Fields
* `initialStepsize` – (`1.0`) and initial step size
* `retraction_method` – (`ExponentialRetraction()`) the rectraction to use, defaults to
  the exponential map
* `contractionFactor` – (`0.95`) exponent for line search reduction
* `sufficientDecrease` – (`0.1`) gain within Armijo's rule
* `last_stepsize` – (`initialstepsize`) the last step size we start the search with
* `linesearch_stopsize` - (`0.0`) a safeguard when to stop the line search
    before the step is numerically zero. This should be combined with [`StopWhenStepSizeLess`](@ref)

# Constructor

    ArmijoLineSearch()

with the Fields above in their order as optional arguments (deprecated).

    ArmijoLineSearch(M)

with the Fields above as keyword arguments and the retraction is set to the default retraciton on `M`.

The constructors return the functor to perform Armijo line search, where two inter
faces are available:
* based on a tuple `(p,o,i)` of a [`GradientProblem`](@ref) `p`, [`Options`](@ref) `o`
  and a current iterate `i`.
* with `(M, x, F, gradFx[,η=-gradFx]) -> s` where [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#Manifold) `M`, a current
  point `x` a function `F`, that maps from the manifold to the reals,
  its gradient (a tangent vector) `gradFx```=\operatorname{grad}F(x)`` at  `x` and an optional
  search direction tangent vector `η=-gradFx` are the arguments.
"""
mutable struct ArmijoLinesearch{TRM<:AbstractRetractionMethod} <: Linesearch
    initialStepsize::Float64
    retraction_method::TRM
    contractionFactor::Float64
    sufficientDecrease::Float64
    last_stepsize::Float64
    linesearch_stopsize::Float64
    @deprecate ArmijoLinesearch(
        s::Float64=1.0,
        r::AbstractRetractionMethod=ExponentialRetraction(),
        contractionFactor::Float64=0.95,
        sufficientDecrease::Float64=0.1,
        linesearch_stopsize::Float64=0.0,
    ) ArmijoLinesearch(
        DefaultManifold(2);
        initial_stepsize=s,
        retraction_method=r,
        contractionFactor=contractionFactor,
        sufficientDecrease=sufficientDecrease,
        linesearch_stopsize=linesearch_stopsize,
    )
    function ArmijoLinesearch(
        M;
        initial_stepsize::Float64=1.0,
        retraction_method::AbstractRetractionMethod=default_retraction_method(M),
        contractionFactor::Float64=0.95,
        sufficientDecrease::Float64=0.1,
        linesearch_stopsize::Float64=0.0,
    )
        return new{typeof(retraction_method)}(
            initial_stepsize,
            retraction_method,
            contractionFactor,
            sufficientDecrease,
            initial_stepsize,
            linesearch_stopsize,
        )
    end
end
function (a::ArmijoLinesearch)(
    p::GradientProblem, o::Options, ::Int, η=-get_gradient(p, o.x); kwargs...
)
    a.last_stepsize = linesearch_backtrack(
        p.M,
        x -> p.cost(p.M, x),
        o.x,
        get_gradient!(p, o.gradient, o.x),
        a.last_stepsize,
        a.sufficientDecrease,
        a.contractionFactor,
        a.retraction_method,
        η;
        stop_step=a.linesearch_stopsize,
    )
    return a.last_stepsize
end
get_initial_stepsize(a::ArmijoLinesearch) = a.initialStepsize

@doc raw"""
    linesearch_backtrack(M, F, x, gradFx, s, decrease, contract, retr, η = -gradFx, f0 = F(x); stop_step=0.)

perform a linesearch for
* a manifold `M`
* a cost function `F`,
* an iterate `x`
* the gradient ``\operatorname{grad}F(x)``
* an initial stepsize `s` usually called ``γ``
* a sufficient `decrease`
* a `contract`ion factor ``σ``
* a `retr`action, which defaults to the `ExponentialRetraction()`
* a search direction ``η = -\operatorname{grad}F(x)``
* an offset, ``f_0 = F(x)``
* a keyword `stop_step` as a minimal step size when to stop
"""
function linesearch_backtrack(
    M::AbstractManifold,
    F::TF,
    x,
    gradFx::T,
    s,
    decrease,
    contract,
    retr::AbstractRetractionMethod=ExponentialRetraction(),
    η::T=-gradFx,
    f0=F(x);
    stop_step=0.0,
) where {TF,T}
    xNew = retract(M, x, s * η, retr)
    fNew = F(xNew)
    search_dir_inner = inner(M, x, η, gradFx)
    extended = false
    while fNew < f0 + decrease * s * search_dir_inner # increase
        extended = true
        s = s / contract
        retract!(M, xNew, x, s * η, retr)
        fNew = F(xNew)
    end
    if extended
        s *= contract  # undo last increase
        retract!(M, xNew, x, s * η, retr)
        fNew = F(xNew)
    end
    while fNew > f0 + decrease * s * search_dir_inner # decrease
        s = contract * s
        retract!(M, xNew, x, s * η, retr)
        fNew = F(xNew)
        (s < stop_step) && break
    end
    return s
end

@doc raw"""
    NonmonotoneLinesearch <: Linesearch

A functor representing a nonmonotone line search using the Barzilai-Borwein step size[^Iannazzo2018]. Together with a gradient descent algorithm
this line search represents the Riemannian Barzilai-Borwein with nonmonotone line-search (RBBNMLS) algorithm. We shifted the order of the algorithm steps from the paper
by Iannazzo and Porcelli so that in each iteration we first find

```math
y_{k} = \operatorname{grad}F(x_{k}) - \operatorname{T}_{x_{k-1} → x_k}(\operatorname{grad}F(x_{k-1}))
```

and

```math
s_{k} = - α_{k-1} * \operatorname{T}_{x_{k-1} → x_k}(\operatorname{grad}F(x_{k-1})),
```

where ``α_{k-1}`` is the step size computed in the last iteration and ``\operatorname{T}`` is a vector transport.
We then find the Barzilai–Borwein step size

```math
α_k^{\text{BB}} = \begin{cases}
\min(α_{\text{max}}, \max(α_{\text{min}}, τ_{k})),  & \text{if } ⟨s_{k}, y_{k}⟩_{x_k} > 0,\\
α_{\text{max}}, & \text{else,}
\end{cases}
```

where

```math
τ_{k} = \frac{⟨s_{k}, s_{k}⟩_{x_k}}{⟨s_{k}, y_{k}⟩_{x_k}},
```

if the direct strategy is chosen,

```math
τ_{k} = \frac{⟨s_{k}, y_{k}⟩_{x_k}}{⟨y_{k}, y_{k}⟩_{x_k}},
```

in case of the inverse strategy and an alternation between the two in case of the
alternating strategy. Then we find the smallest ``h = 0, 1, 2, …`` such that

```math
F(\operatorname{retr}_{x_k}(- σ^h α_k^{\text{BB}} \operatorname{grad}F(x_k)))
\leq
\max_{1 ≤ j ≤ \min(k+1,m)} F(x_{k+1-j}) - γ σ^h α_k^{\text{BB}} ⟨\operatorname{grad}F(x_k), \operatorname{grad}F(x_k)⟩_{x_k},
```

where ``σ`` is a step length reduction factor ``∈ (0,1)``, ``m`` is the number of iterations
after which the function value has to be lower than the current one
and ``γ`` is the sufficient decrease parameter ``∈(0,1)``. We can then find the new stepsize by

```math
α_k = σ^h α_k^{\text{BB}}.
```

[^Iannazzo2018]:
    > B. Iannazzo, M. Porcelli, __The Riemannian Barzilai–Borwein Method with Nonmonotone Line Search and the Matrix Geometric Mean Computation__,
    > In: IMA Journal of Numerical Analysis. Volume 38, Issue 1, January 2018, Pages 495–517,
    > doi [10.1093/imanum/drx015](https://doi.org/10.1093/imanum/drx015)

# Fields
* `initial_stepsize` – (`1.0`) the step size we start the search with
* `linesearch_stopsize` - (`0.0`) a safeguard when to stop the line search
    before the step is numerically zero. This should be combined with [`StopWhenStepSizeLess`](@ref)
* `memory_size` – (`10`) number of iterations after which the cost value needs to be lower than the current one
* `min_stepsize` – (`1e-3`) lower bound for the Barzilai-Borwein step size greater than zero
* `max_stepsize` – (`1e3`) upper bound for the Barzilai-Borwein step size greater than min_stepsize
* `retraction_method` – (`ExponentialRetraction()`) the rectraction to use
* `strategy` – (`direct`) defines if the new step size is computed using the direct, indirect or alternating strategy
* `storage` – (`x`, `gradient`) a [`StoreOptionsAction`](@ref) to store `old_x` and `old_gradient`, the x-value and corresponding gradient of the previous iteration
* `stepsize_reduction` – (`0.5`) step size reduction factor contained in the interval (0,1)
* `sufficient_decrease` – (`1e-4`) sufficient decrease parameter contained in the interval (0,1)
* `vector_transport_method` – (`ParallelTransport()`) the vector transport method to use

# Constructor

    NonmonotoneLinesearch()

with the Fields above in their order as optional arguments (deprecated).

    NonmonotoneLinesearch(M)

with the Fields above in their order as keyword arguments and where the retraction
and vector transport are set to the default ones on `M`, repsectively.

The constructors return the functor to perform nonmonotone line search.
"""
mutable struct NonmonotoneLinesearch{
    TRM<:AbstractRetractionMethod,VTM<:AbstractVectorTransportMethod,T<:AbstractVector
} <: Linesearch
    retraction_method::TRM
    vector_transport_method::VTM
    stepsize_reduction::Float64
    sufficient_decrease::Float64
    min_stepsize::Float64
    max_stepsize::Float64
    initial_stepsize::Float64
    old_costs::T
    strategy::Symbol
    storage::StoreOptionsAction
    linesearch_stopsize::Float64
    function NonmonotoneLinesearch(
        M::AbstractManifold=DefaultManifold(2);
        initial_stepsize::Float64=1.0,
        retraction_method::AbstractRetractionMethod=default_retraction_method(M),
        vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(
            M
        ),
        stepsize_reduction::Float64=0.5,
        sufficient_decrease::Float64=1e-4,
        memory_size::Int=10,
        min_stepsize::Float64=1e-3,
        max_stepsize::Float64=1e3,
        strategy::Symbol=:direct,
        storage::StoreOptionsAction=StoreOptionsAction((:x, :gradient)),
        linesearch_stopsize::Float64=0.0,
    )
        if strategy ∉ [:direct, :inverse, :alternating]
            @warn string(
                "The strategy '",
                strategy,
                "' is not defined. The 'direct' strategy is used instead.",
            )
            strategy = :direct
        end
        if min_stepsize <= 0.0
            throw(
                DomainError(
                    min_stepsize,
                    "The lower bound for the step size min_stepsize has to be greater than zero.",
                ),
            )
        end
        if max_stepsize <= min_stepsize
            throw(
                DomainError(
                    max_stepsize,
                    "The upper bound for the step size max_stepsize has to be greater its lower bound min_stepsize.",
                ),
            )
        end
        if memory_size <= 0
            throw(DomainError(memory_size, "The memory_size has to be greater than zero."))
        end
        return new{
            typeof(retraction_method),typeof(vector_transport_method),Vector{Float64}
        }(
            retraction_method,
            vector_transport_method,
            stepsize_reduction,
            sufficient_decrease,
            min_stepsize,
            max_stepsize,
            initial_stepsize,
            zeros(memory_size),
            strategy,
            storage,
            linesearch_stopsize,
        )
    end
    @deprecate NonmonotoneLinesearch(
        initial_stepsize::Float64,
        retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
        vector_transport_method::AbstractVectorTransportMethod=ParallelTransport(),
        stepsize_reduction::Float64=0.5,
        sufficient_decrease::Float64=1e-4,
        memory_size::Int=10,
        min_stepsize::Float64=1e-3,
        max_stepsize::Float64=1e3,
        strategy::Symbol=:direct,
        storage::StoreOptionsAction=StoreOptionsAction((:x, :gradient)),
        linesearch_stopsize::Float64=0.0,
    ) NonmonotoneLinesearch(
        DefaultManifold(3);
        initial_stepsize=initial_stepsize,
        retraction_method=retraction_method,
        vector_transport_method=vector_transport_method,
        stepsize_reduction=stepsize_reduction,
        sufficient_decrease=sufficient_decrease,
        memory_size=memory_size,
        min_stepsize=min_stepsize,
        max_stepsize=max_stepsize,
        strategy=strategy,
        storage=storage,
        linesearch_stopsize=linesearch_stopsize,
    )
end
function (a::NonmonotoneLinesearch)(
    p::GradientProblem, o::Options, i::Int, η=-get_gradient(p, o.x); kwargs...
)
    if !all(has_storage.(Ref(a.storage), [:x, :gradient]))
        old_x = o.x
        old_gradient = get_gradient(p, o.x)
    else
        old_x, old_gradient = get_storage.(Ref(a.storage), [:x, :gradient])
    end
    update_storage!(a.storage, o)
    return a(p.M, o.x, x -> get_cost(p, x), get_gradient(p, o.x), η, old_x, old_gradient, i)
end
function (a::NonmonotoneLinesearch)(
    M::mT, x, F::TF, gradFx::T, η::T, old_x, old_gradient, iter::Int; kwargs...
) where {mT<:AbstractManifold,TF,T}
    #find the difference between the current and previous gardient after the previous gradient is transported to the current tangent space
    grad_diff =
        gradFx - vector_transport_to(M, old_x, old_gradient, x, a.vector_transport_method)
    #transport the previous step into the tangent space of the current manifold point
    x_diff =
        -a.initial_stepsize *
        vector_transport_to(M, old_x, old_gradient, x, a.vector_transport_method)

    #compute the new Barzilai-Borwein step size
    s1 = inner(M, x, x_diff, grad_diff)
    s2 = inner(M, x, grad_diff, grad_diff)
    s2 = s2 == 0 ? 1.0 : s2
    s3 = inner(M, x, x_diff, x_diff)
    #indirect strategy
    if a.strategy == :inverse
        if s1 > 0
            BarzilaiBorwein_stepsize = min(a.max_stepsize, max(a.min_stepsize, s1 / s2))
        else
            BarzilaiBorwein_stepsize = a.max_stepsize
        end
        #alternating strategy
    elseif a.strategy == :alternating
        if s1 > 0
            if iter % 2 == 0
                BarzilaiBorwein_stepsize = min(a.max_stepsize, max(a.min_stepsize, s1 / s2))
            else
                BarzilaiBorwein_stepsize = min(a.max_stepsize, max(a.min_stepsize, s3 / s1))
            end
        else
            BarzilaiBorwein_stepsize = a.max_stepsize
        end
        #direct strategy
    else
        if s1 > 0
            BarzilaiBorwein_stepsize = min(a.max_stepsize, max(a.min_stepsize, s2 / s1))
        else
            BarzilaiBorwein_stepsize = a.max_stepsize
        end
    end

    memory_size = length(a.old_costs)
    if iter <= memory_size
        a.old_costs[iter] = F(x)
    else
        a.old_costs[1:(memory_size - 1)] = a.old_costs[2:memory_size]
        a.old_costs[memory_size] = F(x)
    end

    #compute the new step size with the help of the Barzilai-Borwein step size
    a.initial_stepsize = linesearch_backtrack(
        M,
        F,
        x,
        gradFx,
        BarzilaiBorwein_stepsize,
        a.sufficient_decrease,
        a.stepsize_reduction,
        a.retraction_method,
        η,
        maximum([a.old_costs[j] for j in 1:min(iter, memory_size)]);
        stop_step=a.linesearch_stopsize,
    )
    return a.initial_stepsize
end

@doc raw"""
    WolfePowellLinesearch <: Linesearch

Do a backtracking linesearch to find a step size ``α`` that fulfils the
Wolfe conditions along a search direction ``η`` starting from ``x``, i.e.

```math
f\bigl( \operatorname{retr}_x(αη) \bigr) ≤ f(x_k) + c_1 α_k ⟨\operatorname{grad}f(x), η⟩_x
\quad\text{and}\quad
\frac{\mathrm{d}}{\mathrm{d}t} f\bigr(\operatorname{retr}_x(tη)\bigr)
\Big\vert_{t=α}
≥ c_2 \frac{\mathrm{d}}{\mathrm{d}t} f\bigl(\operatorname{retr}_x(tη)\bigr)\Big\vert_{t=0}.
```

# Constructors

There exist two constructors, where, when prodivind the manifold `M` as a first (optional)
parameter, its default retraction and vector transport are the default.
In this case the retraction and the vector transport are also keyword arguments for ease of use.
The other constructor is kept for backward compatibility.
Note that the `linesearch_stopsize` to stop for too small stepsizes is only available in the
new signature including `M`.
For the old (deprecated) signature the `linesearch_stopsize` is set to the old hard-coded default of  `1e-12`

    WolfePowellLinesearch(
        retr::AbstractRetractionMethod=ExponentialRetraction(),
        vtr::AbstractVectorTransportMethod=ParallelTransport(),
        c1::Float64=10^(-4),
        c2::Float64=0.999
    )

    WolfePowellLinesearch(
        M,
        c1::Float64=10^(-4),
        c2::Float64=0.999;
        retraction_method = default_retraction_method(M),
        vector_transport_method = default_vector_transport(M),
        linesearch_stopsize = 0.0
    )
"""
mutable struct WolfePowellLinesearch <: Linesearch
    retraction_method::AbstractRetractionMethod
    vector_transport_method::AbstractVectorTransportMethod
    c1::Float64
    c2::Float64
    last_stepsize::Float64
    linesearch_stopsize::Float64

    @deprecate WolfePowellLinesearch(
        retr::AbstractRetractionMethod=ExponentialRetraction(),
        vtr::AbstractVectorTransportMethod=ParallelTransport(),
        c1::Float64=10^(-4),
        c2::Float64=0.999,
    ) WolfePowellLinesearch(
        DefaultManifold(3),
        c1,
        c2,
        retraction_method=retr,
        vector_transport_method=vtr,
        linesearch_stopsize=1e-12,
    )
    function WolfePowellLinesearch(
        M::AbstractManifold,
        c1::Float64=10^(-4),
        c2::Float64=0.999;
        retraction_method::AbstractRetractionMethod=default_retraction_method(M),
        vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(
            M
        ),
        linesearch_stopsize::Float64=0.0,
    )
        return new(
            retraction_method, vector_transport_method, c1, c2, 0.0, linesearch_stopsize
        )
    end
end

function (a::WolfePowellLinesearch)(
    p::P, o::O, ::Int, η=-get_gradient(p, o.x); kwargs...
) where {P<:GradientProblem{T,mT} where {T,mT<:AbstractManifold},O<:Options}
    s = 1.0
    s_plus = 1.0
    s_minus = 1.0
    f0 = get_cost(p, o.x)
    xNew = retract(p.M, o.x, s * η, a.retraction_method)
    fNew = get_cost(p, xNew)
    η_xNew = vector_transport_to(p.M, o.x, η, xNew, a.vector_transport_method)
    if fNew > f0 + a.c1 * s * inner(p.M, o.x, η, o.gradient)
        while (fNew > f0 + a.c1 * s * inner(p.M, o.x, η, o.gradient)) && (s_minus > 10^(-9)) # decrease
            s_minus = s_minus * 0.5
            s = s_minus
            retract!(p.M, xNew, o.x, s * η, a.retraction_method)
            fNew = get_cost(p, xNew)
        end
        s_plus = 2.0 * s_minus
    else
        vector_transport_to!(p.M, η_xNew, o.x, η, xNew, a.vector_transport_method)
        if inner(p.M, xNew, get_gradient(p, xNew), η_xNew) <
            a.c2 * inner(p.M, o.x, η, o.gradient)
            while fNew <= f0 + a.c1 * s * inner(p.M, o.x, η, o.gradient) &&
                (s_plus < 10^(9))# increase
                s_plus = s_plus * 2.0
                s = s_plus
                retract!(p.M, xNew, o.x, s * η, a.retraction_method)
                fNew = get_cost(p, xNew)
            end
            s_minus = s_plus / 2.0
        end
    end
    retract!(p.M, xNew, o.x, s_minus * η, a.retraction_method)
    vector_transport_to!(p.M, η_xNew, o.x, η, xNew, a.vector_transport_method)
    while inner(p.M, xNew, get_gradient(p, xNew), η_xNew) <
          a.c2 * inner(p.M, o.x, η, o.gradient)
        s = (s_minus + s_plus) / 2
        retract!(p.M, xNew, o.x, s * η, a.retraction_method)
        fNew = get_cost(p, xNew)
        if fNew <= f0 + a.c1 * s * inner(p.M, o.x, η, o.gradient)
            s_minus = s
        else
            s_plus = s
        end
        abs(s_plus - s_minus) <= a.linesearch_stopsize && break
        retract!(p.M, xNew, o.x, s_minus * η, a.retraction_method)
        vector_transport_to!(p.M, η_xNew, o.x, η, xNew, a.vector_transport_method)
    end
    s = s_minus
    a.last_stepsize = s
    return s
end

@doc raw"""
    WolfePowellBinaryLinesearch <: Linesearch

A [`Linesearch`](@ref) method that determines a step size `t` fulfilling the Wolfe conditions

based on a binary chop. Let ``η`` be a search direction and ``c1,c_2>0`` be two constants.
Then with

```math
A(t) = f(x_+) ≤ c1 t ⟨\operatorname{grad}f(x), η⟩_{x}
\quad\text{and}\quad
W(t) = ⟨\operatorname{grad}f(x_+), \text{V}_{x_+\gets x}η⟩_{x_+} ≥ c_2 ⟨η, \operatorname{grad}f(x)⟩_x,
```

where ``x_+ = \operatorname{retr}_x(tη)`` is the current trial point, and ``\text{V}`` is a
vector transport, we perform the following Algorithm similar to Algorithm 7 from [^Huang2014]
1. set ``α=0``, ``β=∞`` and ``t=1``.
2. While either ``A(t)`` does not hold or ``W(t)`` does not hold do steps 3-5.
3. If ``A(t)`` fails, set ``β=t``.
4. If ``A(t)`` holds but ``W(t)`` fails, set ``α=t``.
5. If ``β<∞`` set ``t=\frac{α+β}{2}``, otherwise set ``t=2α``.

# Constructors

There exist two constructors, where, when prodivind the manifold `M` as a first (optional)
parameter, its default retraction and vector transport are the default.
In this case the retraction and the vector transport are also keyword arguments for ease of use.
The other constructor is kept for backward compatibility.
Note that the `linesearch_stopsize` to stop for too small stepsizes is only available in the
new signature including `M`, for the first it is set to the old default of `1e-9`.

    WolfePowellBinaryLinesearch(
        retr::AbstractRetractionMethod=ExponentialRetraction(),
        vtr::AbstractVectorTransportMethod=ParallelTransport(),
        c1::Float64=10^(-4),
        c2::Float64=0.999
    )

    WolfePowellLinesearch(
        M,
        c1::Float64=10^(-4),
        c2::Float64=0.999;
        retraction_method = default_retraction_method(M),
        vector_transport_method = default_vector_transport(M),
        linesearch_stopsize = 0.0
    )

[^Huang2014]:
    > Huang, W.: _Optimization algorithms on Riemannian manifolds with applications_,
    > Dissertation, Flordia State University, 2014.
    > [pdf](https://www.math.fsu.edu/~whuang2/pdf/Huang_W_Dissertation_2013.pdf)
"""
mutable struct WolfePowellBinaryLinesearch <: Linesearch
    retraction_method::AbstractRetractionMethod
    vector_transport_method::AbstractVectorTransportMethod
    c1::Float64
    c2::Float64
    last_stepsize::Float64
    linesearch_stopsize::Float64

    @deprecate WolfePowellBinaryLinesearch(
        retr::AbstractRetractionMethod=ExponentialRetraction(),
        vtr::AbstractVectorTransportMethod=ParallelTransport(),
        c1::Float64=10^(-4),
        c2::Float64=0.999,
    ) WolfePowellBinaryLinesearch(
        DefaultManifold(3),
        c1,
        c2,
        retraction_method=retr,
        vector_transport_method=vtr,
        linesearch_stopsize=1e-9,
    )
    function WolfePowellBinaryLinesearch(
        M::AbstractManifold,
        c1::Float64=10^(-4),
        c2::Float64=0.999;
        retraction_method::AbstractRetractionMethod=default_retraction_method(M),
        vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(
            M
        ),
        linesearch_stopsize::Float64=0.0,
    )
        return new(
            retraction_method, vector_transport_method, c1, c2, 0.0, linesearch_stopsize
        )
    end
end

function (a::WolfePowellBinaryLinesearch)(
    p::GradientProblem, o::Options, ::Int, η=-get_gradient(p, o.x); kwargs...
)
    α = 0.0
    β = Inf
    t = 1.0
    f0 = get_cost(p, o.x)
    xNew = retract(p.M, o.x, t * η, a.retraction_method)
    fNew = get_cost(p, xNew)
    η_xNew = vector_transport_to(p.M, o.x, η, xNew, a.vector_transport_method)
    gradient_new = get_gradient(p, xNew)
    nAt = fNew > f0 + a.c1 * t * inner(p.M, o.x, η, o.gradient)
    nWt = inner(p.M, xNew, gradient_new, η_xNew) < a.c2 * inner(p.M, o.x, η, o.gradient)
    while (nAt || nWt) &&
              (t > a.linesearch_stopsize) &&
              ((α + β) / 2 - 1 > a.linesearch_stopsize)
        nAt && (β = t)            # A(t) fails
        (!nAt && nWt) && (α = t)  # A(t) holds but W(t) fails
        t = isinf(β) ? 2 * α : (α + β) / 2
        # Update trial point
        retract!(p.M, xNew, o.x, t * η, a.retraction_method)
        fNew = get_cost(p, xNew)
        gradient_new = get_gradient(p, xNew)
        vector_transport_to!(p.M, η_xNew, o.x, η, xNew, a.vector_transport_method)
        # Update conditions
        nAt = fNew > f0 + a.c1 * t * inner(p.M, o.x, η, o.gradient)
        nWt = inner(p.M, xNew, gradient_new, η_xNew) < a.c2 * inner(p.M, o.x, η, o.gradient)
    end
    a.last_stepsize = t
    return t
end

@doc raw"""
    get_stepsize(p::Problem, o::Options, vars...)

return the stepsize stored within [`Options`](@ref) `o` when solving [`Problem`](@ref) `p`.
This method also works for decorated options and the [`Stepsize`](@ref) function within
the options, by default stored in `o.stepsize`.
"""
function get_stepsize(p::Problem, o::Options, vars...; kwargs...)
    return get_stepsize(p, o, dispatch_options_decorator(o), vars...; kwargs...)
end
function get_stepsize(p::Problem, o::Options, ::Val{true}, vars...; kwargs...)
    return get_stepsize(p, o.options, vars...; kwargs...)
end
function get_stepsize(p::Problem, o::Options, ::Val{false}, vars...; kwargs...)
    return o.stepsize(p, o, vars...; kwargs...)
end

function get_initial_stepsize(p::Problem, o::Options, vars...; kwargs...)
    return get_initial_stepsize(
        p::Problem, o::Options, dispatch_options_decorator(o), vars...; kwargs...
    )
end
function get_initial_stepsize(p::Problem, o::Options, ::Val{true}, vars...)
    return get_initial_stepsize(p, o.options)
end
function get_initial_stepsize(p::Problem, o::Options, ::Val{false}, vars...)
    return get_initial_stepsize(o.stepsize)
end

function get_last_stepsize(p::Problem, o::Options, vars...)
    return get_last_stepsize(p, o, dispatch_options_decorator(o), vars...)
end
function get_last_stepsize(p::Problem, o::Options, ::Val{true}, vars...)
    return get_last_stepsize(p, o.options, vars...)
end
function get_last_stepsize(p::Problem, o::Options, ::Val{false}, vars...)
    return get_last_stepsize(p, o, o.stepsize, vars...)
end
#
# dispatch on stepsize
get_last_stepsize(p::Problem, o::Options, s::Stepsize, vars...) = s(p, o, vars...)
get_last_stepsize(::Problem, ::Options, s::ArmijoLinesearch, ::Any...) = s.last_stepsize
function get_last_stepsize(::Problem, ::Options, s::WolfePowellLinesearch, ::Any...)
    return s.last_stepsize
end
function get_last_stepsize(::Problem, ::Options, s::WolfePowellBinaryLinesearch, ::Any...)
    return s.last_stepsize
end