From 56737102899db3a8cc893162caadda914222b897 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Thu, 2 Apr 2026 18:50:06 -0400 Subject: [PATCH 1/3] Add rotation lifter examples --- popcor/examples/rotation_lifter.py | 75 ++++++++++++++++++++++++++++++ tests/test_rotation_lifter.py | 65 ++++++++++++++++++++++++++ 2 files changed, 140 insertions(+) diff --git a/popcor/examples/rotation_lifter.py b/popcor/examples/rotation_lifter.py index 1421a6c..3f40f76 100644 --- a/popcor/examples/rotation_lifter.py +++ b/popcor/examples/rotation_lifter.py @@ -106,6 +106,7 @@ class RotationLifter(StateLifter): LEVELS: List[str] = ["no", "bm"] HOM: str = "h" VARIABLE_LIST: List[List[str]] = [["h", "c_0"], ["h", "c_0", "c_1"]] + SO2_EXAMPLE_TYPES: Tuple[str, str] = ("A", "B") ADD_DETERMINANT: bool = False NOISE: float = 1e-3 @@ -159,6 +160,80 @@ def sample_theta(self) -> np.ndarray: C[:, i * self.d : (i + 1) * self.d] = Rotation.random().as_matrix() return C + @staticmethod + def so2_theta(angle: float) -> np.ndarray: + """Return a 2x2 SO(2) rotation matrix for a scalar angle in radians.""" + c = np.cos(angle) + s = np.sin(angle) + return np.array([[c, -s], [s, c]]) + + @classmethod + def get_so2_example_coefficients( + cls, example_type: str = "A" + ) -> Tuple[np.ndarray, np.ndarray]: + """Return (quadratic, linear) coefficients for deterministic SO(2) examples. + + The objective is defined on the first column of R(theta), i.e. v=[cos(theta), sin(theta)], + as f(theta) = v.T @ Q2 @ v + q1.T @ v. + + Example types: + - "A": two global minima on the circle (symmetric antipodal minima) + - "B": one global minimum and one local minimum + """ + if example_type == "A": + q2 = np.array([[0.8, 0.05], [0.05, 0.2]]) + q1 = np.array([0.0, 0.0]) + elif example_type == "B": + q2 = np.array([[0.8, 0.05], [0.05, 0.2]]) + q1 = np.array([0.6, 0.0]) + else: + raise ValueError( + f"Unknown example_type {example_type}. Expected one of {cls.SO2_EXAMPLE_TYPES}." + ) + return q2, q1 + + def get_so2_example_Q( + self, example_type: str = "A", output_poly: bool = False + ) -> PolyMatrix | np.ndarray | sp.csr_matrix | sp.csc_matrix: + """Return a deterministic Q for 2D single-rotation examples. + + This helper builds anisotropic quadratic costs on the unit circle and is intended for + tutorial / visualization notebooks. + """ + if not (self.d == 2 and self.n_rot == 1 and self.level == "no"): + raise ValueError( + "SO(2) example Q is currently implemented only for d=2, n_rot=1, level='no'." + ) + + q2, q1 = self.get_so2_example_coefficients(example_type=example_type) + + Q = PolyMatrix(symmetric=True) + + # Vectorization is column-major: [R11, R21, R12, R22]. + q2_full = np.zeros((self.d**2, self.d**2)) + q2_full[:2, :2] = q2 + Q["c_0", "c_0"] = q2_full + + # Linear terms are represented through the homogeneous block. + q1_full = np.zeros((1, self.d**2)) + q1_full[0, :2] = 0.5 * q1 + Q[self.HOM, "c_0"] = q1_full + + if output_poly: + return Q + return Q.get_matrix(self.var_dict) + + def get_so2_example_cost(self, angle: float, example_type: str = "A") -> float: + """Evaluate deterministic SO(2) example cost at a given angle.""" + if not (self.d == 2 and self.n_rot == 1 and self.level == "no"): + raise ValueError( + "SO(2) example cost is currently implemented only for d=2, n_rot=1, level='no'." + ) + theta = self.so2_theta(angle) + x = self.get_x(theta=theta) + Q = self.get_so2_example_Q(example_type=example_type, output_poly=False) + return float(x.T @ Q @ x) + def get_x( self, theta: np.ndarray | None = None, diff --git a/tests/test_rotation_lifter.py b/tests/test_rotation_lifter.py index 931225d..2ff3278 100644 --- a/tests/test_rotation_lifter.py +++ b/tests/test_rotation_lifter.py @@ -135,6 +135,71 @@ def test_solve_local(): print("done") +@pytest.mark.parametrize( + "example_type, expect_equal_minima", + [ + ("A", True), + ("B", False), + ], +) +def test_so2_example_minima_structure(example_type, expect_equal_minima): + """Validate deterministic SO(2) tutorial examples have intended minima patterns.""" + lifter = RotationLifter(d=2, n_rot=1, n_abs=0, n_rel=0, level="no") + + angles = np.linspace(0.0, 2.0 * np.pi, 4001) + costs = np.array([lifter.get_so2_example_cost(t, example_type) for t in angles]) + + minima_idx = [] + for i in range(1, len(angles) - 1): + if costs[i] <= costs[i - 1] and costs[i] <= costs[i + 1]: + if minima_idx and (i - minima_idx[-1]) < 20: + if costs[i] < costs[minima_idx[-1]]: + minima_idx[-1] = i + else: + minima_idx.append(i) + + assert len(minima_idx) == 2 + c1, c2 = float(costs[minima_idx[0]]), float(costs[minima_idx[1]]) + + if expect_equal_minima: + assert abs(c1 - c2) < 1e-5 + else: + assert abs(c1 - c2) > 1e-3 + + +@pytest.mark.parametrize("example_type", ["A", "B"]) +def test_so2_example_sdp_recovery(example_type): + """Check SDP behavior for deterministic SO(2) examples. + + Type A is symmetric with two equal global minima, so the SDP optimum is generally + not rank-1 in the angle-moment block. Type B has a unique global minimum and should + be near rank-1. + """ + lifter = RotationLifter(d=2, n_rot=1, n_abs=0, n_rel=0, level="no") + + Q = lifter.get_so2_example_Q(example_type=example_type) + A_known, b_known = lifter.get_A_known(add_redundant=True) + A_0, b_0 = lifter.get_A0() + constraints = list(zip(A_0 + A_known, b_0 + b_known)) + + X, info_sdp = solve_sdp(Q, constraints, verbose=False) + _, info_rank = rank_project(X, p=1) + + angles = np.linspace(-np.pi, np.pi, 5001) + costs = np.array([lifter.get_so2_example_cost(t, example_type) for t in angles]) + c_min = float(np.min(costs)) + + # SDP lower bound should match the global minimum cost very tightly. + assert abs(float(info_sdp["cost"]) - c_min) < 1e-6 + + if example_type == "A": + # Symmetric two-minimum case should not force rank-1 recovery. + assert info_rank["EVR"] < 10.0 + else: + # Unique minimum case should be near rank-1. + assert info_rank["EVR"] > 1e6 + + @pytest.mark.parametrize("level", ["no", "bm"]) def test_solve_sdp(level): """Solve the SDP relaxation and compare the recovered rotations to ground truth.""" From 831788f36b5b8ac2d30977ef48f862c3da1ae9c1 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Thu, 2 Apr 2026 18:51:39 -0400 Subject: [PATCH 2/3] Update version --- popcor/__init__.py | 2 +- setup.cfg | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/popcor/__init__.py b/popcor/__init__.py index 79b07d0..2d17fcf 100644 --- a/popcor/__init__.py +++ b/popcor/__init__.py @@ -1,4 +1,4 @@ from .auto_template import AutoTemplate from .auto_tight import AutoTight -__version__ = "0.0.2" +__version__ = "0.0.3" diff --git a/setup.cfg b/setup.cfg index e185b01..fef3c87 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = popcor -version = 0.0.2 +version = 0.0.3 authors = [ {name = "Frederike Dümbgen", email = "frederike.duembgen@gmail.com" },] description = POPCOR -- Polynomial Optimization for Certifiably Optimal Robotics From 828d1ae603ab0c84b9f2027f2d4a4b0cbeb9e498 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Fri, 3 Apr 2026 18:12:24 -0400 Subject: [PATCH 3/3] Add changes to CHANGELOG --- CHANGELOG.md | 14 ++++++++++---- environment.yml | 5 +++-- setup.cfg | 1 - 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 14fd333..95bb102 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,10 +5,16 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](http://keepachangelog.com/) and this project adheres to [Semantic Versioning](http://semver.org/). -## [Unreleased] -- 2025-08-23 +## [Unreleased] - -(0.0.0a)= ### Added + +(0.0.3) +## [0.0.3] - 2026-04-03 + +(0.0.3a)= +### Added +- Add example cases to RotationLifter and corresponding tests. - RotationLifter: support to plot 2d frames, improved documentation for rotation conventions - RotationLifter: implemented rank-d ("bm") lifting - RangeOnlyLifter: new landmark sampling methods to fill available space @@ -17,13 +23,13 @@ and this project adheres to [Semantic Versioning](http://semver.org/). - StateLifter: started support for rank-d instead of rank-1 formulations - Added more documentation and type hints -(0.0.0c)= +(0.0.3c)= ### Changed - RangeOnlyLifter: default sampling for landmarks in RO is now the one filling space - RangeOnlyLifter: theta is sampled at least MIN_DIST away from landmarks - Always returning constraints matrices along with right-hand-side values -(0.0.0a)= +(0.0.3f)= ### Fixed - Link fixes in documentation diff --git a/environment.yml b/environment.yml index 15863eb..d2ea41f 100644 --- a/environment.yml +++ b/environment.yml @@ -12,11 +12,12 @@ dependencies: - pytest - black>=23.1 - - seaborn>=0.12.2 - - progressbar2>=4.2.0 + - jupyter>=1.0.0 - mosek + - progressbar2>=4.2.0 - pandas==1.5.3 + - seaborn>=0.12.2 - spatialmath-python=1.1.11 - suitesparse # used by sparseqr - pip: diff --git a/setup.cfg b/setup.cfg index fef3c87..8f55669 100644 --- a/setup.cfg +++ b/setup.cfg @@ -21,7 +21,6 @@ packages = install_requires = numpy>=1.23 scipy>=1.9 - cvxpy>=1.3 pymanopt>=2.1 chompack>=2.3 autograd