wrapper.jl

using LBFGSB
using Test

# define a function
function f(x)
    n = length(x)
    y = 0.25 * (x[1] - 1)^2
    for i = 2:n
        y += (x[i] - x[i-1]^2)^2
    end
    4y
end

# and its gradient function that maps a vector x to a vector z
function g!(z, x)
    n = length(x)
    t₁ = x[2] - x[1]^2
    z[1] = 2 * (x[1] - 1) - 1.6e1 * x[1] * t₁
    for i = 2:n-1
        t₂ = t₁
        t₁ = x[i+1] - x[i]^2
        z[i] = 8 * t₂ - 1.6e1 * x[i] * t₁
    end
    z[n] = 8 * t₁
end

# and a combined version (which can be useful if f and g! share some computations)
function fg!(z, x)
    g!(z, x)
    return f(x)
end

@testset "wrapper" begin
    # the first argument is the dimension of the largest problem to be solved
    # the second argument is the maximum number of limited memory corrections
    optimizer = L_BFGS_B(1024, 17)
    n = 25  # the dimension of the problem
    x = fill(Cdouble(3e0), n)  # the initial guess
    # set up bounds
    bounds = zeros(3, n)
    for i = 1:n
        bounds[1,i] = 2  # represents the type of bounds imposed on the variables:
                         #  0->unbounded, 1->only lower bound, 2-> both lower and upper bounds, 3->only upper bound
        bounds[2,i] = isodd(i) ? 1e0 : -1e2  #  the lower bound on x, of length n.
        bounds[3,i] = 1e2  #  the upper bound on x, of length n.
    end
    # call the optimizer
    # - m: the maximum number of variable metric corrections used to define the limited memory matrix
    # - factr: the iteration will stop when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch,
    #     where epsmch is the machine precision, which is automatically generated by the code. Typical values for factr:
    #     1e12 for low accuracy, 1e7 for moderate accuracy, 1e1 for extremely high accuracy
    # - pgtol: the iteration will stop when max{|proj g_i | i = 1, ..., n} <= pgtol where pg_i is the ith component of the projected gradient
    # - iprint: controls the frequency of output. iprint < 0 means no output:
    #     iprint = 0 print only one line at the last iteration
    #     0 < iprint < 99 print also f and |proj g| every iprint iterations
    #     iprint = 99 print details of every iteration except n-vectors
    #     iprint = 100 print also the changes of active set and final x
    #     iprint > 100 print details of every iteration including x and g
    # - maxfun: the maximum number of function evaluations
    # - maxiter: the maximum number of iterations
    fout, xout = optimizer(f, g!, x, bounds, m=5, factr=1e7, pgtol=1e-5, iprint=-1, maxfun=15000, maxiter=15000)
    @test fout ≈ 1.083490083518441e-9
    
    #testing the original wrapper:
    function g(x)
        n = length(x)
        z = zeros(n)
        t₁ = x[2] - x[1]^2
        z[1] = 2 * (x[1] - 1) - 1.6e1 * x[1] * t₁
        for i = 2:n-1
            t₂ = t₁
            t₁ = x[i+1] - x[i]^2
            z[i] = 8 * t₂ - 1.6e1 * x[i] * t₁
        end
        z[n] = 8 * t₁
        z
    end
    func(x) = (f(x), g(x))
    fout2, xout2 = optimizer(func, x, bounds, m=5, factr=1e7, pgtol=1e-5, iprint=-1, maxfun=15000, maxiter=15000)
    @test fout2 ≈ 1.083490083518441e-9

    # testing the combined in-place version
    fout3, xout3 = optimizer(LBFGSB.InPlace(fg!), x, bounds, m=5, factr=1e7, pgtol=1e-5, iprint=-1, maxfun=15000, maxiter=15000)
    @test fout3 ≈ 1.083490083518441e-9
end

@testset "lbfgsb" begin
    #  0->unbounded, 1->only lower bound, 2-> both lower and upper bounds, 3->only upper bound
    lb=[-Inf,-Inf,-2.,6f0]
    ub=[Inf,3.,2,Inf]
    btyps=[0,3,2,1]
    
    @test LBFGSB.typ_bnd.(lb,ub) == btyps
    
    @test_throws ErrorException LBFGSB.typ_bnd(6,-1.)
    
    n,m=4,6
    o,b=LBFGSB._opt_bounds(n,m,lb,ub)
    
    @test size(b)==(3,n)
    @test b[2,:] == lb
    @test b[3,:] == ub
    @test b[1,:] == btyps
    
    od=L_BFGS_B(n,m)
    @test all(getproperty(od,i)==getproperty(o,i) for i in fieldnames(L_BFGS_B))
    
    fout,xout=lbfgsb(x->0.5sum(abs,x),(g,x)->g.=x,100*rand(3),lb=[0.5,-Inf,-0.1],ub=[Inf,0.2,Inf])
    @test fout ≈ 0.25
    @test xout ≈ [0.5,0,0]

    fout2,xout2=lbfgsb(LBFGSB.InPlace((g,x)->(g.=x; 0.5sum(abs,x))),100*rand(3),lb=[0.5,-Inf,-0.1],ub=[Inf,0.2,Inf])
    @test fout2 ≈ 0.25
    @test xout2 ≈ [0.5,0,0]
 end