Add wrapper for combined fg! functions

README.md

@@ -72,6 +72,30 @@ julia> fout, xout = optimizer(f, g!, x, bounds, m=5, factr=1e7, pgtol=1e-5, ipri
 
 # if your func needs extra arguments, a closure can be used to do the trick:
 # optimizer(λ->func(λ, extra_args...), x, bounds, m=5, factr=1e7, pgtol=1e-5, iprint=-1, maxfun=15000, maxiter=15000)
+
+# in some problems, the function and gradient share many compuations
+# in which case it can be more efficient to compute them together.
+# for the example above, we can get a bit of re-use by combining `f` and `g!`
+# as follows
+function fg!(z, x)
+    n = length(x)
+    y = 0.25 * (x[1] - 1)^2
+
+    t₁ = x[2] - x[1]^2
+    y += t₁^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
+        y += t₁^2
+        z[i] = 8 * t₂ - 1.6e1 * x[i] * t₁
+    end
+    z[n] = 8 * t₁
+    return 4y
+end
+# to use this function, wrap it in the `LBFGSB.InPlace` type like so:
+julia> fout, xout = optimizer(LBFGSB.InPlace(fg!), x, bounds, m=5, factr=1e7, pgtol=1e-5, iprint=-1, maxfun=15000, maxiter=15000)
+(1.0834900835178617e-9, [1.000002325877184, 1.0000030446388726, 1.0000045143926763, 1.0000050772587235, 1.00000581930634, 1.0000070253663256, 1.0000083438107845, 1.0000120681556417, 1.0000196904590828, 1.0000374605368185  …  1.002601944378863, 1.0052112931969248, 1.0104498202560672, 1.0210095629751155, 1.0424600004275002, 1.0867234743980112, 1.1809686991764257, 1.394691851289783, 1.9451626718243045, 3.783662657064461])
 ```
 This package also provides `lbfgsb`, a convenience function that computes the `bounds` input matrix and the `optimizer` internally. The inputs of `lbfgs` are the same as above, except for lower/upper bounds specified optionally as keyword parameters `lb`/`ub` (defaulting to `-Inf/Inf`, i.e. unbounded problem).  With the example above, an equivalent call to `lbfgsb` would be:
 ```julia

src/wrapper.jl

@@ -116,6 +116,53 @@ function (obj::L_BFGS_B)(func, x0::AbstractVector, bounds::AbstractMatrix;
     end
 end
 
+struct InPlace{F}
+    func_grad!::F
+end
+
+function (obj::L_BFGS_B)(func::InPlace, x0::AbstractVector, bounds::AbstractMatrix;
+    m=10, factr=1e7, pgtol=1e-5, iprint=-1, maxfun=15000, maxiter=15000)
+    x = copy(x0)
+    n = length(x)
+    f = 0.0
+    # clean up
+    fill!(obj.task, Cuchar(' '))
+    fill!(obj.csave, Cuchar(' '))
+    fill!(obj.lsave, zero(Cint))
+    fill!(obj.isave, zero(Cint))
+    fill!(obj.dsave, zero(Cdouble))
+    fill!(obj.wa, zero(Cdouble))
+    fill!(obj.iwa, zero(Cint))
+    fill!(obj.g, zero(Cdouble))
+    fill!(obj.nbd, zero(Cint))
+    fill!(obj.l, zero(Cdouble))
+    fill!(obj.u, zero(Cdouble))
+    # set bounds
+    for i = 1:n
+        obj.nbd[i] = bounds[1,i]
+        obj.l[i] = bounds[2,i]
+        obj.u[i] = bounds[3,i]
+    end
+    # start
+    obj.task[1:5] = b"START"
+    while true
+        iprint == 1 && flush(stdout)  # PR 12
+        setulb(n, m, x, obj.l, obj.u, obj.nbd, f, obj.g, factr, pgtol, obj.wa,
+               obj.iwa, obj.task, iprint, obj.csave, obj.lsave, obj.isave, obj.dsave)
+        if obj.task[1:2] == b"FG"
+            f = func.func_grad!(obj.g, x)
+        elseif obj.task[1:5] == b"NEW_X"
+            if obj.isave[30] ≥ maxiter
+                obj.task[1:43] = b"STOP: TOTAL NO. of ITERATIONS REACHED LIMIT"
+            elseif obj.isave[34] ≥ maxfun
+                obj.task[1:52] = b"STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT"
+            end
+        else
+            return f, x
+        end
+    end
+end
+
 function typ_bnd(lb,ub)
         lb > ub && error("Inconsistent bounds")
 

test/wrapper.jl

@@ -24,6 +24,12 @@ function g!(z, x)
     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
@@ -72,6 +78,10 @@ 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
@@ -98,4 +108,8 @@ end
     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