Add optimization test with FD derivatives

pyoptsparse/pyCONMIN/pyCONMIN.py

@@ -212,7 +212,7 @@ def cnmngrad(n1, n2, x, f, g, ct, df, a, ic, nac):
             dabfun = self.getOption("DABFUN")
 
             itrm = self.getOption("ITRM")
-            nfeasct = self.getOption("ITRM")
+            nfeasct = self.getOption("NFEASCT")
             nfdg = 1  # User will supply all gradients
 
             # Counters for functions and gradients

pyoptsparse/pyOpt_gradient.py

@@ -44,6 +44,12 @@ def __init__(self, optProb: Optimization, sensType: str, sensStep: float = None,
                 self.sensStep = 1e-40j
         else:
             self.sensStep = sensStep
+
+        # Complex step divides by the imaginary part of the step, so a purely
+        # real step would silently yield NaN gradients.
+        if self.sensType == "cs" and np.imag(self.sensStep) == 0:
+            raise ValueError(f"The complex step size must have a nonzero imaginary part, got {self.sensStep}.")
+
         self.sensMode = sensMode
         self.comm = comm
 

pyoptsparse/testing/pyOpt_testing.py

@@ -356,46 +356,42 @@ def check_hist_file(self, tol):
             for varName in self.DVs:
                 assert_allclose(val[varName].flatten(), self.xStar[self.sol_index][varName], atol=tol, rtol=tol)
 
-    def optimize_with_hotstart(self, tol, optOptions=None, x0=None):
+    def optimize_with_hotstart(self, tol, x0=None, **kwargs):
         """
         This code will perform 4 optimizations, one real opt and three restarts.
         In this process, it will check various combinations of storeHistory and hotStart filenames.
         It will also call `check_hist_file` after the first optimization.
         """
         # we use a non-default starting point to test that the hotstart works
         # even if it does not match optProb initial values
-        sol = self.optimize(storeHistory=True, optOptions=optOptions, setDV=x0)
+        sol = self.optimize(storeHistory=True, setDV=x0, **kwargs)
         self.assert_solution_allclose(sol, tol)
         self.assertGreater(self.nf, 0)
         if self.optName in GRAD_BASED_OPTIMIZERS:
             self.assertGreater(self.ng, 0)
         self.check_hist_file(tol)
 
         # re-optimize with hotstart
-        sol = self.optimize(storeHistory=False, hotStart=True, optOptions=optOptions)
+        sol = self.optimize(storeHistory=False, hotStart=True, **kwargs)
         self.assert_solution_allclose(sol, tol)
         # we should have zero actual function/gradient evaluations
         self.assertEqual(self.nf, 0)
         self.assertEqual(self.ng, 0)
         # another test with hotstart, this time with storeHistory = hotStart
-        sol = self.optimize(storeHistory=True, hotStart=True, optOptions=optOptions)
+        sol = self.optimize(storeHistory=True, hotStart=True, **kwargs)
         self.assert_solution_allclose(sol, tol)
         # we should have zero actual function/gradient evaluations
         self.assertEqual(self.nf, 0)
         self.assertEqual(self.ng, 0)
         # another test with hotstart, this time with a non-existing history file
         # this will perform a cold start
-        self.optimize(storeHistory=True, hotStart="notexisting.hst", optOptions=optOptions)
+        self.optimize(storeHistory=True, hotStart="notexisting.hst", **kwargs)
         self.assertGreater(self.nf, 0)
         if self.optName in GRAD_BASED_OPTIMIZERS:
             self.assertGreater(self.ng, 0)
         self.check_hist_file(tol)
         # final test with hotstart, this time with a different storeHistory
-        sol = self.optimize(
-            storeHistory=f"{self.id()}_new_hotstart.hst",
-            hotStart=True,
-            optOptions=optOptions,
-        )
+        sol = self.optimize(storeHistory=f"{self.id()}_new_hotstart.hst", hotStart=True, **kwargs)
         self.assert_solution_allclose(sol, tol)
         # we should have zero actual function/gradient evaluations
         self.assertEqual(self.nf, 0)

tests/test_gradient.py

@@ -0,0 +1,88 @@
+# Standard Python modules
+import unittest
+
+# External modules
+import numpy as np
+from numpy.testing import assert_allclose
+from parameterized import parameterized
+
+# First party modules
+from pyoptsparse import Optimization
+from pyoptsparse.pyOpt_gradient import Gradient
+
+# Base point at which we evaluate the derivatives
+X0 = {"x": [1.5, -2.0], "y": 0.5}
+
+# Analytic Jacobian at X0
+ANALYTIC = {
+    "obj": {"x": [3.0, -8.0], "y": 3.0},
+    "c": {"x": [[-2.0, 1.5]], "y": 1.0},
+}
+
+
+def objfunc(xdict):
+    """
+    Obj = x0^2 + 2*x1^2 + 3*y^2
+    c   = x0*x1 + y
+    """
+    x, y = xdict["x"], xdict["y"]
+    funcs = {}
+    funcs["obj"] = x[0] ** 2 + 2 * x[1] ** 2 + 3 * y**2
+    funcs["c"] = x[0] * x[1] + y
+    return funcs, False
+
+
+def build_optProb(objfun=objfunc, xScale=1.0, conScale=1.0):
+    optProb = Optimization("grad-test", objfun)
+    optProb.addVarGroup("x", 2, lower=-10, upper=10, value=X0["x"], scale=xScale)
+    optProb.addVar("y", lower=-10, upper=10, value=X0["y"], scale=xScale)
+    optProb.addObj("obj")
+    optProb.addCon("c", lower=-100, upper=100, scale=conScale)
+    optProb.finalize()
+    return optProb
+
+
+def assert_sens_matches_analytic(funcsSens, atol):
+    for funcKey, perGroup in ANALYTIC.items():
+        for dvGroup, expected in perGroup.items():
+            assert_allclose(funcsSens[funcKey][dvGroup], expected, atol=atol)
+
+
+class TestGradient(unittest.TestCase):
+    @parameterized.expand(["fd", "fdr", "cd", "cdr", "cs"])
+    def test_mode_matches_analytic(self, sensType):
+        optProb = build_optProb()
+        funcs, _ = objfunc(X0)
+        grad = Gradient(optProb, sensType=sensType)
+        funcsSens, fail = grad(X0, funcs)
+        self.assertFalse(fail)
+        atol = 1e-12 if sensType == "cs" else 1e-5
+        assert_sens_matches_analytic(funcsSens, atol=atol)
+
+        # test that we get real derivs for cs
+        if sensType == "cs":
+            for funcKey in ANALYTIC:
+                for dvGroup in ANALYTIC[funcKey]:
+                    self.assertFalse(np.iscomplexobj(funcsSens[funcKey][dvGroup]))
+
+    def test_failed_eval(self):
+        def always_fail(xdict):
+            funcs, _ = objfunc(xdict)
+            return funcs, True
+
+        optProb = build_optProb(objfun=always_fail)
+        funcs, _ = objfunc(X0)
+        grad = Gradient(optProb, sensType="fd")
+        _, fail = grad(X0, funcs)
+        self.assertTrue(fail)
+
+    def test_scaling(self):
+        optProb = build_optProb(xScale=7.0, conScale=0.3)
+        funcs, _ = objfunc(X0)
+        grad = Gradient(optProb, sensType="cs")
+        funcsSens, _ = grad(X0, funcs)
+        assert_sens_matches_analytic(funcsSens, 1e-12)
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/test_optProb.py

@@ -18,6 +18,7 @@
 
 # First party modules
 from pyoptsparse import OPT, Optimization
+from pyoptsparse.pyOpt_utils import INFINITY, convertToCSR, convertToDense
 from pyoptsparse.testing.pyOpt_testing import assert_optProb_size
 
 
@@ -295,5 +296,123 @@ def test_parallel_add(self):
         self.assertEqual(allConNames[0], allConNames[1])
 
 
+class TestScaling(unittest.TestCase):
+    def setUp(self):
+        # Distinct, non-trivial per-element scales and offsets so that any
+        # mixed-up indexing or row/column confusion shows up.
+        self.xScale = {"x": [2.0, 0.5, 4.0], "y": [10.0, 0.1]}
+        self.xOffset = {"x": [1.0, -2.0, 0.5], "y": [0.0, 3.0]}
+        self.objScale = 3.0
+        self.conScaleVals = {"c1": [5.0, 0.2], "c2": [7.0]}
+
+        def objfunc(xdict):
+            # Never actually called in these tests, but required by the API.
+            return {"obj": 0.0, "c1": np.zeros(2), "c2": np.zeros(1)}, False
+
+        optProb = Optimization("scaling-test", objfunc)
+        optProb.addVarGroup("x", 3, lower=-10, upper=10, scale=self.xScale["x"], offset=self.xOffset["x"])
+        optProb.addVarGroup("y", 2, lower=-10, upper=10, scale=self.xScale["y"], offset=self.xOffset["y"])
+        optProb.addObj("obj", scale=self.objScale)
+        optProb.addConGroup("c1", 2, lower=-1, upper=1, scale=self.conScaleVals["c1"])
+        optProb.addConGroup("c2", 1, lower=-1, upper=1, scale=self.conScaleVals["c2"])
+        optProb.finalize()
+
+        self.optProb = optProb
+        self.ndvs = optProb.ndvs
+        self.nCon = optProb.nCon
+        # invXScale = 1/scale, in DV order (x then y)
+        self.invXScale = optProb.invXScale
+        # conScale in natural (un-reordered) order: c1, c2
+        self.conScale = optProb.conScale
+
+    def test_finalize_populated_scales(self):
+        assert_allclose(self.invXScale, 1.0 / np.array([2.0, 0.5, 4.0, 10.0, 0.1]))
+        assert_allclose(self.conScale, [5.0, 0.2, 7.0])
+        assert_allclose(self.optProb.xOffset, [1.0, -2.0, 0.5, 0.0, 3.0])
+
+    def test_mapX_roundtrip_and_formula(self):
+        rng = np.random.default_rng(0)
+        x_user = rng.uniform(-5, 5, self.ndvs)
+        x_opt = self.optProb._mapXtoOpt(x_user)
+        # x_opt = (x_user - offset) / invXScale
+        assert_allclose(x_opt, (x_user - self.optProb.xOffset) / self.invXScale)
+        # round trip
+        assert_allclose(self.optProb._mapXtoUser(x_opt), x_user)
+
+    def test_mapObjGrad(self):
+        # Objective gradient mapping: g_opt = g_user * s_f * invXScale (column/chain-rule scaling).
+        rng = np.random.default_rng(1)
+        gobj = rng.uniform(-3, 3, (self.optProb.nObj, self.ndvs))
+        gobj_orig = gobj.copy()
+        gobj_opt = self.optProb._mapObjGradtoOpt(gobj)
+        assert_allclose(gobj_opt, gobj * self.objScale * self.invXScale)
+        # the method must not mutate its input
+        assert_allclose(gobj, gobj_orig)
+
+    def test_mapConJac_formula_and_roundtrip(self):
+        # Build an arbitrary dense Jacobian of the right shape and convert to CSR.
+        rng = np.random.default_rng(2)
+        dense = rng.uniform(-2, 2, (self.nCon, self.ndvs))
+        jac = convertToCSR(dense)
+
+        # _mapConJactoOpt works in place: J_opt = diag(conScale) . J . diag(invXScale)
+        self.optProb._mapConJactoOpt(jac)
+        expected = np.diag(self.conScale) @ dense @ np.diag(self.invXScale)
+        assert_allclose(convertToDense(jac), expected)
+
+        # _mapConJactoUser must invert it back to the original.
+        self.optProb._mapConJactoUser(jac)
+        assert_allclose(convertToDense(jac), dense)
+
+    def test_mapObj_value_roundtrip(self):
+        f_user = 2.5
+        f_opt = self.optProb._mapObjtoOpt(f_user)
+        assert_allclose(f_opt, f_user * self.objScale)
+        assert_allclose(self.optProb._mapObjtoUser(f_opt), f_user)
+
+    def test_mapCon_value_roundtrip(self):
+        c_user = [1.0, -2.0, 3.0]
+        c_opt = self.optProb._mapContoOpt(c_user)
+        assert_allclose(c_opt, c_user * self.conScale)
+        assert_allclose(self.optProb._mapContoUser(c_opt), c_user)
+
+    def test_combined_scale_and_offset(self):
+        """A DV group with both a non-unit scale and a non-zero offset is the
+        classic place to get the order of operations wrong.
+        """
+
+        def objfunc(xdict):
+            return {"obj": 0.0}, False
+
+        optProb = Optimization("edge", objfunc)
+        optProb.addVarGroup("x", 2, lower=-10, upper=10, scale=4.0, offset=3.0)
+        optProb.addObj("obj")
+        optProb.finalize()
+
+        x_user = [3.0, 7.0]  # note x_user[0] == offset
+        x_opt = optProb._mapXtoOpt(x_user)
+        # (x - 3) * 4
+        assert_allclose(x_opt, [0.0, 16.0])
+        assert_allclose(optProb._mapXtoUser(x_opt), x_user)
+
+    def test_infinite_bounds_not_scaled(self):
+        """INFINITY bounds must remain unbounded; scale/offset must not turn
+        them into finite numbers in the assembled bounds.
+        """
+
+        def objfunc(xdict):
+            return {"obj": 0.0}, False
+
+        optProb = Optimization("inf", objfunc)
+        optProb.addVarGroup("x", 1, lower=None, upper=None, scale=10.0, offset=5.0)
+        optProb.addObj("obj")
+        optProb.finalize()
+
+        var = optProb.variables["x"][0]
+        # Variable stores scaled bounds; unbounded sides stay at exactly +/- INFINITY.
+        self.assertEqual(var.lower, -INFINITY)
+        self.assertEqual(var.upper, INFINITY)
+
+
 if __name__ == "__main__":
     unittest.main()

tests/test_sphere.py

@@ -1,6 +1,7 @@
 """Test solution of Sphere problem"""
 
 # Standard Python modules
+from itertools import product
 import unittest
 
 # External modules
@@ -57,6 +58,9 @@ class TestSphere(OptTest):
             "maxGen": 100,
             "seed": 123,
         },
+        "CONMIN": {  # CONMIN diverges when gradient is near zero, here we stop on first optimal iterate
+            "ITRM": 1,
+        },
         "SNOPT": {
             "Major iterations limit": 10,
         },
@@ -101,6 +105,17 @@ def test_optimization(self, optName):
         optOptions = self.optOptions.get(optName, {})
         self.optimize_with_hotstart(self.tol[optName], optOptions=optOptions)
 
+    @parameterized.expand(
+        product(ALL_OPTIMIZERS, ["fd", "fdr", "cd", "cdr", "cs"]),
+        name_func=lambda f, n, p: f"{f.__name__}_{p.args[0]}_{p.args[1]}",
+    )
+    def test_optimization_approx_deriv(self, optName, sens):
+        self.optName = optName
+        self.setup_optProb()
+        optOptions = self.optOptions.get(optName, {})
+        sol = self.optimize(optOptions=optOptions, sens=sens)
+        self.assert_solution_allclose(sol, self.tol[optName])
+
     @parameterized.expand(["filtersqp", "funnelsqp"])
     def test_uno_presets(self, preset):
         self.optName = "Uno"