# fix eta k cchoice
# fix quadstep?

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


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

"""
   _U3

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

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

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

"""
  bracket

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

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

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

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