ActiveBox fails when reaching maximum number of iterations about nlsolvers.jl HOT 8 CLOSED

function asolve(prob::OptimizationProblem, x0, scheme::ActiveBox, options::OptimizationOptions)
    t0 = time()

    x0, B0 = x0, (false*x0*x0' +I)
    lower, upper = NLSolvers.bounds(prob)
    ϵbounds = mapreduce(b->(b[2] - b[1])/2, min, zip(lower, upper)) # [1, pp. 100: 5.41]

    !any(clamp.(x0, lower, upper) .!= x0) || error("Initial guess not in the feasible region")

    linesearch = NLSolvers.ArmijoBertsekas()
    mstyle = NLSolvers.OutOfPlace()
    objvars = NLSolvers.prepare_variables(prob, scheme, x0, copy(x0), B0)
    f0, ∇f0 = objvars.fz, norm(objvars.∇fz, Inf) # use user norm
    fz, ∇fz = objvars.fz, objvars.∇fz # use user norm
    fx, ∇fx = fz, copy(∇fz)
    B = B0
    x, z = copy(x0), copy(x0)
    Tf = typeof(fz)
    is_first=false
    Ix = Diagonal(z.*0 .+ 1)
    for iter = 1:options.maxiter
        x = copy(z)
        fx = copy(fz)
        ∇fx = copy(∇fz)
        
        ϵ = min(norm(clamp.(x.-∇fx, lower, upper).-x), ϵbounds) # Kelley 5.41 and just after (83) in [1]
        activeset = NLSolvers.is_ϵ_active.(x, lower, upper, ∇fx, ϵ)
        Hhat = NLSolvers.diagrestrict.(B, activeset, activeset', Ix)
        # Update current gradient and calculate the search direction
        unconstrained_d = -Hhat\∇fx  # solve Bd = -∇fx  #use find_direction here
        d = clamp.(x .+ unconstrained_d, lower, upper).-x #clamp unconstrained direction
        d_scale = minimum(d./unconstrained_d) 
        #scale the direction, the original direction is the same.
        !iszero(d_scale) && (d = unconstrained_d .* d_scale)
        #if the scale is zero, we are actually in a boundary, and we can only grow in one direction.
        
        φ = NLSolvers._lineobjective(mstyle, prob, ∇fz, z, x, d, fz, dot(∇fz, d))
        
        # Perform line search along d
        # Also returns final step vector and update the state
        α, f_α, ls_success = NLSolvers.find_steplength(mstyle, linesearch, φ, Tf(1), ∇fz, activeset, lower, upper, x, d, ∇fx, activeset)

        # # Calculate final step vector and update the state
        if !ls_success
            return NLSolvers.ConvergenceInfo(scheme, (prob=prob, B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=iter, time=time()-t0), options)
        end
        s = @. α * d
        z = @. x + s
        s = clamp.(z, lower, upper) - x
        z = x + s
        # Update approximation
        fz, ∇fz, B, s, y = NLSolvers.update_obj(prob.objective, s, ∇fx, z, ∇fz, B, Newton(), is_first)
        normx = norm(x.-clamp.(x.-∇fz, lower, upper), Inf)
        if normx < options.g_abstol
            return NLSolvers.ConvergenceInfo(scheme, (prob=prob, B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=iter, time=time()-t0), options)
        end
    end
    #@show z.-min.(upper, max.(z.-∇fz, lower))
  z, fz, options.maxiter
  iter =  options.maxiter
  return NLSolvers.ConvergenceInfo(scheme, (prob=prob, B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=iter, time=time()-t0), options)

        Hhat = diagrestrict.(B, activeset, activeset', Ix)
        # Update current gradient and calculate the search direction
        HhatChol = cholesky(Positive, Hhat)
        d = clamp.(x.-HhatChol\∇fx, lower, upper).-x # solve Bd = -∇fx  #use find_direction here