diff --git a/pyoptsparse/pyCONMIN/pyCONMIN.py b/pyoptsparse/pyCONMIN/pyCONMIN.py index a2b2eb4a..82cfedff 100644 --- a/pyoptsparse/pyCONMIN/pyCONMIN.py +++ b/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 diff --git a/pyoptsparse/pyOpt_gradient.py b/pyoptsparse/pyOpt_gradient.py index 42b9b6aa..ef0d74c0 100644 --- a/pyoptsparse/pyOpt_gradient.py +++ b/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 diff --git a/pyoptsparse/testing/pyOpt_testing.py b/pyoptsparse/testing/pyOpt_testing.py index 6c1ea6c5..3f8738f4 100644 --- a/pyoptsparse/testing/pyOpt_testing.py +++ b/pyoptsparse/testing/pyOpt_testing.py @@ -356,7 +356,7 @@ 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. @@ -364,7 +364,7 @@ def optimize_with_hotstart(self, tol, optOptions=None, x0=None): """ # 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: @@ -372,30 +372,26 @@ def optimize_with_hotstart(self, tol, optOptions=None, x0=None): 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) diff --git a/tests/test_gradient.py b/tests/test_gradient.py new file mode 100644 index 00000000..de277a8c --- /dev/null +++ b/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() diff --git a/tests/test_optProb.py b/tests/test_optProb.py index 6e41fb71..d2743259 100644 --- a/tests/test_optProb.py +++ b/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() diff --git a/tests/test_sphere.py b/tests/test_sphere.py index 2dc77014..2661454f 100644 --- a/tests/test_sphere.py +++ b/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" diff --git a/tests/test_utils.py b/tests/test_utils.py new file mode 100644 index 00000000..aa5607d6 --- /dev/null +++ b/tests/test_utils.py @@ -0,0 +1,254 @@ +""" +Unit tests for the sparse-matrix utilities in pyOpt_utils. +""" + +# Standard Python modules +import unittest + +# External modules +import numpy as np +from numpy.testing import assert_allclose, assert_array_equal + +# First party modules +from pyoptsparse.pyOpt_utils import ( + _broadcast_to_array, + convertToCOO, + convertToCSC, + convertToCSR, + convertToDense, + extractRows, + mapToCSC, + mapToCSR, + scaleColumns, + scaleRows, +) + +# All three sparse fixtures represent the same 3x3 matrix: +# [[1, 0, 2], +# [0, 3, 0], +# [4, 0, 5]] +# +# _COO is intentionally NOT in row-major order to verify that conversion +# routines handle unsorted input correctly. +# +# _CSR is the expected output of convertToCSR(_COO): elements within each +# row appear in COO-arrival order, not sorted by column. +# row 0: (col 0, 1.0), (col 2, 2.0) +# row 1: (col 1, 3.0) +# row 2: (col 2, 5.0), (col 0, 4.0) <- arrival order from _COO +# +# _CSC is the expected output of convertToCSC(_CSR): elements within each +# column appear in row-scan order from the CSR pass. +# col 0: (row 0, 1.0), (row 2, 4.0) +# col 1: (row 1, 3.0) +# col 2: (row 0, 2.0), (row 2, 5.0) +_DENSE = np.array( + [ + [1.0, 0.0, 2.0], + [0.0, 3.0, 0.0], + [4.0, 0.0, 5.0], + ] +) +_COO = { + "coo": [ + [2, 0, 1, 0, 2], + [2, 0, 1, 2, 0], + [5.0, 1.0, 3.0, 2.0, 4.0], + ], + "shape": [3, 3], +} +_CSR = { + "csr": [ + [0, 2, 3, 5], + [0, 2, 1, 2, 0], + [1.0, 2.0, 3.0, 5.0, 4.0], + ], + "shape": [3, 3], +} +_CSC = { + "csc": [ + [0, 2, 3, 5], + [0, 2, 1, 0, 2], + [1.0, 4.0, 3.0, 2.0, 5.0], + ], + "shape": [3, 3], +} + + +class TestConvertToDense(unittest.TestCase): + def test_from_coo(self): + assert_allclose(convertToDense(_COO), _DENSE) + + def test_from_csr(self): + assert_allclose(convertToDense(_CSR), _DENSE) + + def test_from_csc(self): + assert_allclose(convertToDense(_CSC), _DENSE) + + def test_from_dense_array(self): + assert_allclose(convertToDense(_DENSE), _DENSE) + + +class TestConvertToCOO(unittest.TestCase): + def test_from_csr(self): + coo = convertToCOO(_CSR) + rows, cols, data = coo["coo"] + assert_array_equal(rows, [0, 0, 1, 2, 2]) + assert_array_equal(cols, [0, 2, 1, 2, 0]) + assert_allclose(data, [1.0, 2.0, 3.0, 5.0, 4.0]) + self.assertEqual(coo["shape"], [3, 3]) + + def test_from_csc(self): + coo = convertToCOO(_CSC) + rows, cols, data = coo["coo"] + assert_array_equal(rows, [0, 2, 1, 0, 2]) + assert_array_equal(cols, [0, 0, 1, 2, 2]) + assert_allclose(data, [1.0, 4.0, 3.0, 2.0, 5.0]) + self.assertEqual(coo["shape"], [3, 3]) + + def test_from_dense_array(self): + coo = convertToCOO(_DENSE) + rows, cols, data = coo["coo"] + reconstructed = np.zeros(_DENSE.shape) + for r, c, v in zip(rows, cols, data, strict=True): + reconstructed[r, c] = v + assert_allclose(reconstructed, _DENSE) + + def test_idempotent(self): + self.assertIs(convertToCOO(_COO), _COO) + + +class TestConvertToCSR(unittest.TestCase): + def test_from_coo(self): + # COO arrives unordered; elements land in row buckets in COO-arrival order. + # row 0: (col 0, 1.0) then (col 2, 2.0); row 2: (col 2, 5.0) then (col 0, 4.0) + csr = convertToCSR(_COO) + rowp, col_idx, data = csr["csr"] + assert_array_equal(rowp, [0, 2, 3, 5]) + assert_array_equal(col_idx, [0, 2, 1, 2, 0]) + assert_allclose(data, [1.0, 2.0, 3.0, 5.0, 4.0]) + self.assertEqual(csr["shape"], [3, 3]) + + def test_from_csc(self): + # _CSC expands to COO in column-scan order; that COO then feeds the CSR builder. + # row 0: (col 0, 1.0) then (col 2, 2.0); row 2: (col 0, 4.0) then (col 2, 5.0) + csr = convertToCSR(_CSC) + rowp, col_idx, data = csr["csr"] + assert_array_equal(rowp, [0, 2, 3, 5]) + assert_array_equal(col_idx, [0, 2, 1, 0, 2]) + assert_allclose(data, [1.0, 2.0, 3.0, 4.0, 5.0]) + self.assertEqual(csr["shape"], [3, 3]) + + def test_from_dense_array(self): + csr = convertToCSR(_DENSE) + self.assertIn("csr", csr) + assert_allclose(convertToDense(csr), _DENSE) + + def test_idempotent(self): + self.assertIs(convertToCSR(_CSR), _CSR) + + +class TestConvertToCSC(unittest.TestCase): + def test_from_coo(self): + # Converts COO -> CSR -> CSC. Column-scan order from _CSR: + # col 0: (row 0, 1.0),(row 2, 4.0); col 1: (row 1, 3.0); col 2: (row 0, 2.0),(row 2, 5.0) + csc = convertToCSC(_COO) + colp, row_idx, data = csc["csc"] + assert_array_equal(colp, [0, 2, 3, 5]) + assert_array_equal(row_idx, [0, 2, 1, 0, 2]) + assert_allclose(data, [1.0, 4.0, 3.0, 2.0, 5.0]) + self.assertEqual(csc["shape"], [3, 3]) + + def test_from_csr(self): + # _CSR -> CSC: same column-scan order, same result as test_from_coo. + csc = convertToCSC(_CSR) + colp, row_idx, data = csc["csc"] + assert_array_equal(colp, [0, 2, 3, 5]) + assert_array_equal(row_idx, [0, 2, 1, 0, 2]) + assert_allclose(data, [1.0, 4.0, 3.0, 2.0, 5.0]) + self.assertEqual(csc["shape"], [3, 3]) + + def test_from_dense_array(self): + csc = convertToCSC(_DENSE) + self.assertIn("csc", csc) + assert_allclose(convertToDense(csc), _DENSE) + + def test_idempotent(self): + self.assertIs(convertToCSC(_CSC), _CSC) + + +class TestIndexMaps(unittest.TestCase): + """mapToCSR/mapToCSC return index arrays into the original data; the subtle + part is that the permutation must reproduce the matrix exactly. + """ + + def test_mapToCSR(self): + row_p, col_idx, idx_data = mapToCSR(_COO) + data = np.asarray(_COO["coo"][2]) + csr = {"csr": [row_p, col_idx, data[idx_data]], "shape": [3, 3]} + assert_allclose(convertToDense(csr), _DENSE) + # last entry of the row pointer is the nnz + self.assertEqual(row_p[-1], 5) + + def test_mapToCSC(self): + row_idx, col_p, idx_data = mapToCSC(_COO) + data = np.asarray(_COO["coo"][2]) + csc = {"csc": [col_p, row_idx, data[idx_data]], "shape": [3, 3]} + assert_allclose(convertToDense(csc), _DENSE) + self.assertEqual(col_p[-1], 5) + + +class TestRowColScaling(unittest.TestCase): + def test_scaleRows(self): + csr = convertToCSR(_DENSE) + factor = np.array([10.0, 100.0, 1000.0]) + scaleRows(csr, factor) + assert_allclose(convertToDense(csr), np.diag(factor) @ _DENSE) + + def test_scaleColumns(self): + csr = convertToCSR(_DENSE) + factor = np.array([2.0, 3.0, 4.0]) + scaleColumns(csr, factor) + assert_allclose(convertToDense(csr), _DENSE @ np.diag(factor)) + + def test_scale_wrong_length_raises(self): + csr = convertToCSR(_DENSE) + with self.assertRaises(ValueError): + scaleRows(csr, np.array([1.0, 2.0])) + with self.assertRaises(ValueError): + scaleColumns(csr, np.array([1.0, 2.0])) + + +class TestExtractRows(unittest.TestCase): + def test_extractRows(self): + csr = convertToCSR(_DENSE) + sub = extractRows(csr, [0, 2]) + assert_allclose(convertToDense(sub), _DENSE[[0, 2], :]) + self.assertEqual(sub["shape"], [2, 3]) + + +class TestBroadcastToArray(unittest.TestCase): + def test_scalar_broadcast(self): + out = _broadcast_to_array("scale", 2.0, 4) + assert_array_equal(out, np.full(4, 2.0)) + + def test_array_passthrough(self): + out = _broadcast_to_array("scale", [1.0, 2.0, 3.0], 3) + assert_array_equal(out, np.array([1.0, 2.0, 3.0])) + + def test_wrong_length_raises(self): + with self.assertRaises(ValueError): + _broadcast_to_array("scale", [1.0, 2.0], 3) + + def test_none_disallowed_by_default(self): + with self.assertRaises(ValueError): + _broadcast_to_array("lower", None, 3) + + def test_none_allowed_when_requested(self): + out = _broadcast_to_array("lower", None, 3, allow_none=True) + self.assertEqual(len(out), 3) + self.assertTrue(all(v is None for v in out)) + + +if __name__ == "__main__": + unittest.main()