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_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()