ENH: rebased version of gh-14279 (#15157)

Addresses issue gh-2063 (currently closed, but unresolved). It adds the
callback function to the COBYLA optimizer in `minimize`.

I fully realize that in the COBYLA algorithm, each iteration is just a
single function evaluation, and that the user is free to do whatever
they want to in terms of callback behavior in that function. However,
currently the implementation of how the callback argument is handled is
just inconsistent between COBYLA and the other methods. Furthermore, a
user may want to separate their objective function from the callback.
This PR simply brings the COBYLA method in line with all the other
methods.


Co-authored-by: Daniël de Vries <danieldevries6@gmail.com>

scipy/optimize/_cobyla_py.py

@@ -36,7 +36,8 @@ def wrapper(*args, **kwargs):
 
 @synchronized
 def fmin_cobyla(func, x0, cons, args=(), consargs=None, rhobeg=1.0,
-                rhoend=1e-4, maxfun=1000, disp=None, catol=2e-4):
+                rhoend=1e-4, maxfun=1000, disp=None, catol=2e-4,
+                *, callback=None):
     """
     Minimize a function using the Constrained Optimization By Linear
     Approximation (COBYLA) method. This method wraps a FORTRAN
@@ -70,6 +71,9 @@ def fmin_cobyla(func, x0, cons, args=(), consargs=None, rhobeg=1.0,
         Maximum number of function evaluations.
     catol : float, optional
         Absolute tolerance for constraint violations.
+    callback : callable, optional
+        Called after each iteration, as ``callback(x)``, where ``x`` is the
+        current parameter vector.
 
     Returns
     -------
@@ -171,7 +175,8 @@ def fmin_cobyla(func, x0, cons, args=(), consargs=None, rhobeg=1.0,
             'tol': rhoend,
             'disp': disp,
             'maxiter': maxfun,
-            'catol': catol}
+            'catol': catol,
+            'callback': callback}
 
     sol = _minimize_cobyla(func, x0, args, constraints=con,
                            **opts)
@@ -182,7 +187,8 @@ def fmin_cobyla(func, x0, cons, args=(), consargs=None, rhobeg=1.0,
 @synchronized
 def _minimize_cobyla(fun, x0, args=(), constraints=(),
                      rhobeg=1.0, tol=1e-4, maxiter=1000,
-                     disp=False, catol=2e-4, **unknown_options):
+                     disp=False, catol=2e-4, callback=None,
+                     **unknown_options):
     """
     Minimize a scalar function of one or more variables using the
     Constrained Optimization BY Linear Approximation (COBYLA) algorithm.
@@ -201,6 +207,9 @@ def _minimize_cobyla(fun, x0, args=(), constraints=(),
         Maximum number of function evaluations.
     catol : float
         Tolerance (absolute) for constraint violations
+    callback : callable, optional
+        Called after each iteration, as ``callback(x)``, where ``x`` is the
+        current parameter vector.
 
     """
     _check_unknown_options(unknown_options)
@@ -256,10 +265,14 @@ def calcfc(x, con):
             i += size
         return f
 
+    def wrapped_callback(x):
+        if callback is not None:
+            callback(np.copy(x))
+
     info = np.zeros(4, np.float64)
     xopt, info = cobyla.minimize(calcfc, m=m, x=np.copy(x0), rhobeg=rhobeg,
                                   rhoend=rhoend, iprint=iprint, maxfun=maxfun,
-                                  dinfo=info)
+                                  dinfo=info, callback=wrapped_callback)
 
     if info[3] > catol:
         # Check constraint violation

scipy/optimize/_minimize.py

@@ -563,9 +563,6 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
     if meth == 'cobyla' and bounds is not None:
         warn('Method %s cannot handle bounds.' % method,
              RuntimeWarning)
-    # - callback
-    if (meth in ('cobyla',) and callback is not None):
-        warn('Method %s does not support callback.' % method, RuntimeWarning)
     # - return_all
     if (meth in ('l-bfgs-b', 'tnc', 'cobyla', 'slsqp') and
             options.get('return_all', False)):
@@ -687,7 +684,8 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
         res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,
                             **options)
     elif meth == 'cobyla':
-        res = _minimize_cobyla(fun, x0, args, constraints, **options)
+        res = _minimize_cobyla(fun, x0, args, constraints, callback=callback,
+                                **options)
     elif meth == 'slsqp':
         res = _minimize_slsqp(fun, x0, args, jac, bounds,
                               constraints, callback=callback, **options)

scipy/optimize/cobyla/cobyla.pyf

@@ -8,11 +8,16 @@ python module _cobyla__user__routines
             double precision intent(out) :: f
             double precision intent(in), dimension(m), depend(m) :: con
         end subroutine calcfc
+        subroutine callback(n,m,x)
+            integer intent(in,hide) :: n
+            integer intent(in,hide) :: m
+            double precision dimension(n),depend(n),intent(in) :: x
+        end subroutine callback
     end interface _cobyla_user_interface
 end python module _cobyla__user__routines
 python module _cobyla ! in 
     interface  ! in :__cobyla
-        subroutine minimize(calcfc,n,m,x,rhobeg,rhoend,iprint,maxfun,w,iact,dinfo)
+        subroutine minimize(calcfc,n,m,x,rhobeg,rhoend,iprint,maxfun,w,iact,dinfo,callback)
             fortranname cobyla
             use _cobyla__user__routines
             external calcfc
@@ -26,6 +31,7 @@ python module _cobyla ! in
             double precision dimension(n*(3*n+2*m+11)+4*m+6), intent(cache,hide),depend(n,m) :: w
             integer dimension(m + 1),intent(cache,hide),depend(m) :: iact
             double precision dimension(4), intent(in,out) :: dinfo
+            external callback
         end subroutine minimize
     end interface 
 end python module _cobyla

scipy/optimize/cobyla/cobyla2.f

@@ -1,9 +1,9 @@
 C------------------------------------------------------------------------ 
 C
       SUBROUTINE COBYLA (CALCFC, N,M,X,RHOBEG,RHOEND,IPRINT,MAXFUN,
-     & W,IACT, DINFO)
+     & W,IACT, DINFO, CALLBACK)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      EXTERNAL CALCFC
+      EXTERNAL CALCFC,CALLBACK
       DIMENSION X(*),W(*),IACT(*), DINFO(*)
 C
 C     This subroutine minimizes an objective function F(X) subject to M
@@ -74,16 +74,16 @@ SUBROUTINE COBYLA (CALCFC, N,M,X,RHOBEG,RHOEND,IPRINT,MAXFUN,
       IWORK=IDX+N
       CALL COBYLB (CALCFC,N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(ICON),
      1  W(ISIM),W(ISIMI),W(IDATM),W(IA),W(IVSIG),W(IVETA),W(ISIGB),
-     2  W(IDX),W(IWORK),IACT,DINFO)
+     2  W(IDX),W(IWORK),IACT,DINFO,CALLBACK)
       RETURN
       END
 C------------------------------------------------------------------------------
       SUBROUTINE COBYLB (CALCFC,N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,
-     1  CON,SIM,SIMI,DATMAT,A,VSIG,VETA,SIGBAR,DX,W,IACT,DINFO)
+     1  CON,SIM,SIMI,DATMAT,A,VSIG,VETA,SIGBAR,DX,W,IACT,DINFO,CALLBACK)
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION X(*),CON(*),SIM(N,*),SIMI(N,*),DATMAT(MPP,*),
      1  A(N,*),VSIG(*),VETA(*),SIGBAR(*),DX(*),W(*),IACT(*),DINFO(*)
-      EXTERNAL CALCFC
+      EXTERNAL CALCFC,CALLBACK
 C
 C     Set the initial values of some parameters. The last column of SIM holds
 C     the optimal vertex of the current simplex, and the preceding N columns
@@ -134,7 +134,10 @@ SUBROUTINE COBYLB (CALCFC,N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,
 C     instructions are also used for calling CALCFC during the iterations of
 C     the algorithm.
 C
-   40 IF (NFVALS .GE. MAXFUN .AND. NFVALS .GT. 0) THEN
+   40 IF (NFVALS .GT. 0) THEN
+         CALL CALLBACK(N,M,X)
+      END IF
+      IF (NFVALS .GE. MAXFUN .AND. NFVALS .GT. 0) THEN
          IF (IPRINT .GE. 1) PRINT 50
    50      FORMAT (/3X,'Return from subroutine COBYLA because the ',
      1        'MAXFUN limit has been reached.')
@@ -561,6 +564,7 @@ SUBROUTINE COBYLB (CALCFC,N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,
           PRINT 70, NFVALS,F,RESMAX,(X(I),I=1,IPTEM)
           IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
       END IF
+      CALL CALLBACK(N,M,X)
       MAXFUN=NFVALS
       DINFO(2)=DBLE(NFVALS)
       DINFO(3)=DBLE(F)

scipy/optimize/tests/test_cobyla.py

@@ -1,7 +1,7 @@
 import math
 import numpy as np
 
-from numpy.testing import assert_allclose, assert_
+from numpy.testing import assert_allclose, assert_, assert_array_equal
 
 from scipy.optimize import fmin_cobyla, minimize
 
@@ -29,16 +29,31 @@ def test_simple(self):
         assert_allclose(x, self.solution, atol=1e-4)
 
     def test_minimize_simple(self):
+        class Callback:
+            def __init__(self):
+                self.n_calls = 0
+                self.last_x = None
+
+            def __call__(self, x):
+                self.n_calls += 1
+                self.last_x = x
+
+        callback = Callback()
+
         # Minimize with method='COBYLA'
         cons = ({'type': 'ineq', 'fun': self.con1},
                 {'type': 'ineq', 'fun': self.con2})
         sol = minimize(self.fun, self.x0, method='cobyla', constraints=cons,
-                       options=self.opts)
+                       callback=callback, options=self.opts)
         assert_allclose(sol.x, self.solution, atol=1e-4)
         assert_(sol.success, sol.message)
         assert_(sol.maxcv < 1e-5, sol)
         assert_(sol.nfev < 70, sol)
         assert_(sol.fun < self.fun(self.solution) + 1e-3, sol)
+        assert_(sol.nfev == callback.n_calls,
+                "Callback is not called exactly once for every function eval.")
+        assert_array_equal(sol.x, callback.last_x,
+                           "Last design vector sent to the callback is not equal to returned value.")
 
     def test_minimize_constraint_violation(self):
         np.random.seed(1234)

scipy/optimize/tests/test_optimize.py

@@ -1110,10 +1110,6 @@ def test_minimize_callback_copies_array(self, method):
         # Check that arrays passed to callbacks are not modified
         # inplace by the optimizer afterward
 
-        # cobyla doesn't have callback
-        if method == 'cobyla':
-            return
-
         if method in ('fmin_tnc', 'fmin_l_bfgs_b'):
             func = lambda x: (optimize.rosen(x), optimize.rosen_der(x))
         else: