diff --git a/.gitignore b/.gitignore index 26a192e..2763e0d 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,4 @@ mosek_output.tmp bin/act docs/build/doctest docs/build/ +.mypy_cache diff --git a/CHANGELOG.md b/CHANGELOG.md index 95bb102..66985b6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,14 +5,27 @@ 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] - +## [Unreleased] ### Added -(0.0.3) +### Changed + +### Fixed + +## [0.0.4] - 2026-04-08 + +### Changed +- Fix the ground truth value for Poly4 type B +- Creating fixed examples is now handled by staticmethod create_example() + - Added this for: Poly4Lifter, Poly6Lifter, RotationLifter, RangeOnlySqLifter, RangeOnlyNsqLifter + - Pending: other lifters: create this functionality +- Created two plotting types: plot_cost and plot_setup + - Added this for: Poly4Lifter, Poly6Lifter, RotationLifter, RangeOnlySqLifter, RangeOnlyNsqLifter + - Pending: make sure all others also conform to this new structure + ## [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 @@ -23,19 +36,16 @@ 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.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.3f)= ### Fixed - Link fixes in documentation ## [0.0.2] - 2025-06-04 -(0.0.2a)= ### Added - RangeOnlyLifter base class for 2D/3D range-only localization (base_lifters/range_only_lifter.py) @@ -43,12 +53,10 @@ and this project adheres to [Semantic Versioning](http://semver.org/). - RangeOnlySqLifter: specialized lifter for squared-distance-based range-only localization - GitHub Pages deployment step to documentation.yml using peaceiris/actions-gh-pages for automatic documentation publishing -(0.0.2r)= ### Removed - docs/build/ files -(0.0.2f)= ### Fixed - Corrected return type annotations in get_grad and get_hess methods in state_lifter.py from float to np.ndarray diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 304560a..a938af1 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -13,7 +13,7 @@ Please try to the best of your abilities to: ## Adding a new lifter class -You can start with the [ExampleLifter](../popcor/examples/ExampleLifter.py) skeleton, and feel free to add more functionalities depending on the nature of the problem. You can also consider adding a new base class similar to [RobustPoseLifter](../popcor/base_lifters/RobustPoseLifter.py) or [StereoLifter](../popcor/base_lifters/StereoLifter.py) if you want to create multiple new lifters that all share similar functionalities. +You can start with the `example_lifter.py` skeleton in `popcor/examples/`, and feel free to add more functionalities depending on the nature of the problem. You can also consider adding a new base class similar to `robust_pose_lifter.py` or `stereo_lifter.py` (both in `popcor/base_lifters/`) if you want to create multiple new lifters that all share similar functionalities. ## Adding new functionalities diff --git a/docs/source/_static/ro_nsq_cost_A.png b/docs/source/_static/ro_nsq_cost_A.png new file mode 100644 index 0000000..1ab0f8e Binary files /dev/null and b/docs/source/_static/ro_nsq_cost_A.png differ diff --git a/docs/source/_static/ro_nsq_cost_B.png b/docs/source/_static/ro_nsq_cost_B.png new file mode 100644 index 0000000..ede0a13 Binary files /dev/null and b/docs/source/_static/ro_nsq_cost_B.png differ diff --git a/docs/source/_static/ro_nsq_cost_C.png b/docs/source/_static/ro_nsq_cost_C.png new file mode 100644 index 0000000..2aa5499 Binary files /dev/null and b/docs/source/_static/ro_nsq_cost_C.png differ diff --git a/docs/source/_static/ro_nsq_setup_A.png b/docs/source/_static/ro_nsq_setup_A.png new file mode 100644 index 0000000..ada44c6 Binary files /dev/null and b/docs/source/_static/ro_nsq_setup_A.png differ diff --git a/docs/source/_static/ro_nsq_setup_B.png b/docs/source/_static/ro_nsq_setup_B.png new file mode 100644 index 0000000..86b0626 Binary files /dev/null and b/docs/source/_static/ro_nsq_setup_B.png differ diff --git a/docs/source/_static/ro_nsq_setup_C.png b/docs/source/_static/ro_nsq_setup_C.png new file mode 100644 index 0000000..b7f8bed Binary files /dev/null and b/docs/source/_static/ro_nsq_setup_C.png differ diff --git a/docs/source/_static/ro_sq_cost_A.png b/docs/source/_static/ro_sq_cost_A.png new file mode 100644 index 0000000..a30f14a Binary files /dev/null and b/docs/source/_static/ro_sq_cost_A.png differ diff --git a/docs/source/_static/ro_sq_cost_B.png b/docs/source/_static/ro_sq_cost_B.png new file mode 100644 index 0000000..6203c61 Binary files /dev/null and b/docs/source/_static/ro_sq_cost_B.png differ diff --git a/docs/source/_static/ro_sq_cost_C.png b/docs/source/_static/ro_sq_cost_C.png new file mode 100644 index 0000000..6d19059 Binary files /dev/null and b/docs/source/_static/ro_sq_cost_C.png differ diff --git a/docs/source/_static/ro_sq_setup_A.png b/docs/source/_static/ro_sq_setup_A.png new file mode 100644 index 0000000..db3d251 Binary files /dev/null and b/docs/source/_static/ro_sq_setup_A.png differ diff --git a/docs/source/_static/ro_sq_setup_B.png b/docs/source/_static/ro_sq_setup_B.png new file mode 100644 index 0000000..d6d2919 Binary files /dev/null and b/docs/source/_static/ro_sq_setup_B.png differ diff --git a/docs/source/_static/ro_sq_setup_C.png b/docs/source/_static/ro_sq_setup_C.png new file mode 100644 index 0000000..0bd831d Binary files /dev/null and b/docs/source/_static/ro_sq_setup_C.png differ diff --git a/docs/source/conf.py b/docs/source/conf.py index 01b37b8..5e4a59f 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -14,11 +14,14 @@ import os import sys -import popcor - # The module you're documenting (assumes you've added the project root dir to sys.path) sys.path.insert(0, os.path.abspath("../..")) +try: + import popcor +except (ImportError, ModuleNotFoundError): + popcor = None + # -- Project information ----------------------------------------------------- project = "POPCOR" diff --git a/docs/source/examples/standard.rst b/docs/source/examples/standard.rst index e2bf32c..1458aec 100644 --- a/docs/source/examples/standard.rst +++ b/docs/source/examples/standard.rst @@ -12,6 +12,94 @@ Range-Only Localization :undoc-members: :show-inheritance: +RangeOnlySqLifter Example Plots +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The below plots are generated by calling the following lines of code. +You can find the full code in the ``__main__`` section of the corresponding module. + +.. code-block:: python + + lifter = RangeOnlySqLifter.create_example(example_type="A") # choose A, B or C + lifter.plot_cost() + + +.. figure:: /_static/ro_sq_cost_A.png + :alt: RangeOnlySqLifter cost heatmap type A + :align: center + :figclass: align-center + + RangeOnlySqLifter cost heatmap (Type A) + +.. figure:: /_static/ro_sq_cost_B.png + :alt: RangeOnlySqLifter cost heatmap type B + :align: center + :figclass: align-center + + RangeOnlySqLifter cost heatmap (Type B) + +.. figure:: /_static/ro_sq_cost_C.png + :alt: RangeOnlySqLifter cost heatmap type C + :align: center + :figclass: align-center + + RangeOnlySqLifter cost heatmap (Type C) + + +RangeOnlyNsqLifter Example Plots +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + +The below plots are generated by calling the following lines of code. +You can find the full code in the ``__main__`` section of the corresponding module. + +.. code-block:: python + + lifter = RangeOnlyNsqLifter.create_example(example_type="A") # choose A, B or C + lifter.plot_cost() + +.. figure:: /_static/ro_nsq_cost_A.png + :alt: RangeOnlyNsqLifter cost heatmap type A + :align: center + :figclass: align-center + + RangeOnlyNsqLifter cost heatmap (Type A) + +.. figure:: /_static/ro_nsq_setup_A.png + :alt: RangeOnlyNsqLifter setup plot type A + :align: center + :figclass: align-center + + RangeOnlyNsqLifter setup (Type A) + +.. figure:: /_static/ro_nsq_cost_B.png + :alt: RangeOnlyNsqLifter cost heatmap type B + :align: center + :figclass: align-center + + RangeOnlyNsqLifter cost heatmap (Type B) + +.. figure:: /_static/ro_nsq_setup_B.png + :alt: RangeOnlyNsqLifter setup plot type B + :align: center + :figclass: align-center + + RangeOnlyNsqLifter setup (Type B) + +.. figure:: /_static/ro_nsq_cost_C.png + :alt: RangeOnlyNsqLifter cost heatmap type C + :align: center + :figclass: align-center + + RangeOnlyNsqLifter cost heatmap (Type C) + +.. figure:: /_static/ro_nsq_setup_C.png + :alt: RangeOnlyNsqLifter setup plot type C + :align: center + :figclass: align-center + + RangeOnlyNsqLifter setup (Type C) + Stereo-Camera Localization -------------------------- diff --git a/docs/source/examples/toy.rst b/docs/source/examples/toy.rst index 69666d1..e785ae0 100644 --- a/docs/source/examples/toy.rst +++ b/docs/source/examples/toy.rst @@ -4,6 +4,14 @@ Toy Examples Univariate Polynomials ---------------------- +The below plots are generated by calling the following lines of code. +You can find the full code in the ``__main__`` section of the corresponding module. + +.. code-block:: python + + lifter = Poly4Lifter.create_example(example_type="A") # choose A or B + lifter.plot_cost() + .. autoclass:: popcor.examples.Poly4Lifter :undoc-members: :show-inheritance: diff --git a/docs/source/quickstart.rst b/docs/source/quickstart.rst index e78a249..7be811d 100644 --- a/docs/source/quickstart.rst +++ b/docs/source/quickstart.rst @@ -139,7 +139,7 @@ More information on how to use AutoTemplate can be found :ref:`here =1.6.0 + - -r docs/requirements.txt - -e . diff --git a/environment_small.yml b/environment_small.yml index 3e5ae60..485b29c 100644 --- a/environment_small.yml +++ b/environment_small.yml @@ -11,4 +11,5 @@ dependencies: - suitesparse # used by sparseqr - pip: + - sparseqr>=1.6.0 - -e . diff --git a/popcor/__init__.py b/popcor/__init__.py index 2d17fcf..4f24f1a 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.3" +__version__ = "0.0.4" diff --git a/popcor/auto_tight.py b/popcor/auto_tight.py index 4f1dec5..ca90a97 100644 --- a/popcor/auto_tight.py +++ b/popcor/auto_tight.py @@ -303,7 +303,7 @@ def generate_matrices_simple( elif normalize and sparse: Ai /= splinalg.norm(Ai, ord="fro") if sparse: - assert isinstance(Ai, (sp.csr_matrix, sp.csc_matrix)) + assert isinstance(Ai, (sp.csr_array, sp.csc_array)) Ai.eliminate_zeros() else: assert isinstance(Ai, np.ndarray) @@ -338,7 +338,7 @@ def generate_matrices( ai = lifter.get_reduced_a(bi=basis[i], var_subset=var_dict, sparse=True) basis_reduced.append(ai) basis_reduced = sp.vstack(basis_reduced) - assert isinstance(basis_reduced, sp.csr_matrix) + assert isinstance(basis_reduced, sp.csr_array) if AutoTight.REDUCE_DEPENDENT: bad_idx = find_dependent_columns(basis_reduced.T, tolerance=1e-6) else: @@ -355,7 +355,7 @@ def generate_matrices( elif normalize and sparse: Ai /= splinalg.norm(Ai, ord="fro") if sparse: - assert isinstance(Ai, (sp.csr_matrix, sp.csc_matrix)) + assert isinstance(Ai, (sp.csr_array, sp.csc_array)) Ai.eliminate_zeros() else: assert isinstance(Ai, np.ndarray) diff --git a/popcor/base_lifters/poly_lifter.py b/popcor/base_lifters/poly_lifter.py index 7ee7836..d644290 100644 --- a/popcor/base_lifters/poly_lifter.py +++ b/popcor/base_lifters/poly_lifter.py @@ -49,14 +49,34 @@ def local_solver(self, t0, y=None, *args, **kwargs): info = {"success": sol.success, "cost": sol.fun} return sol.x, info, sol.fun - def plot(self, thetas, label=None): + def plot(self, thetas=None, label=None, estimates=None): from popcor.utils.plotting_tools_poly import coordinate_system fig, ax = coordinate_system() - ys = [self.get_cost(t) for t in thetas] - ax.plot(thetas, ys, label=label) - ymin = min(-max(ys) / 3, min(ys)) - ax.set_ylim(ymin, max(ys)) + + # Handle the case where estimates is provided but thetas is not + if thetas is None and estimates is None: + raise ValueError("Either thetas or estimates must be provided") + + if thetas is not None: + ys = [self.get_cost(t) for t in thetas] + ax.plot(thetas, ys, label=label) + ymin = min(-max(ys) / 3, min(ys)) + ax.set_ylim(ymin, max(ys)) + + # Plot estimates if provided + if estimates is not None: + for est_label, theta_est in estimates.items(): + cost_est = self.get_cost(theta_est) + ax.scatter( + [theta_est], + [cost_est], + label=est_label, + s=100, + marker="x", + linewidth=2, + ) + return fig, ax def __repr__(self): diff --git a/popcor/base_lifters/range_only_lifter.py b/popcor/base_lifters/range_only_lifter.py index 258606f..f637e34 100644 --- a/popcor/base_lifters/range_only_lifter.py +++ b/popcor/base_lifters/range_only_lifter.py @@ -1,5 +1,6 @@ import warnings from abc import ABC, abstractmethod +from typing import Any import autograd.numpy as anp import matplotlib.pylab as plt @@ -37,6 +38,8 @@ class RangeOnlyLifter(StateLifter, ABC): NOISE = 1e-2 # std deviation of distance noise SCALE = 2.0 # size of the region of intereist: [0, SCALE]^d MIN_DIST = 1e-2 # minimum distance between landmarks and positions + EXAMPLE_TYPES: tuple[str, str, str] = ("A", "B", "C") + example_type: str | None = None def __init__( self, @@ -101,6 +104,40 @@ def create_good(n_positions, n_landmarks, d=2): ) return landmarks, theta + @classmethod + def create_example( + cls, + example_type: str = "A", + n_positions: int = 3, + n_landmarks: int = 4, + d: int = 2, + ) -> "RangeOnlyLifter": + """Create a deterministic tutorial setup selected by example type. + + - "A": good initialization. + - "B": bad initialization. + - "C": fixed bad initialization. + """ + if example_type not in cls.EXAMPLE_TYPES: + raise ValueError( + f"Unknown example_type {example_type}. Expected one of {cls.EXAMPLE_TYPES}." + ) + + if example_type == "A": + landmarks, theta = RangeOnlyLifter.create_good(n_positions, n_landmarks, d) + elif example_type == "B": + landmarks, theta = RangeOnlyLifter.create_bad(n_positions, n_landmarks, d) + else: + landmarks, theta = RangeOnlyLifter.create_bad_fixed( + n_positions, n_landmarks, d + ) + + lifter = cls(n_positions, n_landmarks, d) + lifter.overwrite_theta(theta) + lifter.landmarks = landmarks + lifter.example_type = example_type + return lifter + @staticmethod def sample_landmarks_filling_space(n_landmarks, d=2): landmarks = np.random.rand(n_landmarks, d) @@ -344,7 +381,7 @@ def fun(x): cost = sol.fun info["max res"] = np.max(np.abs(residuals)) hess = self.get_hess(that, y) - assert isinstance(hess, sp.csc_matrix) + assert isinstance(hess, sp.csc_array) eigs = np.linalg.eigvalsh(hess.toarray()) info["cond Hess"] = eigs[-1] / eigs[0] else: @@ -392,6 +429,106 @@ def plot(self, y=None, xlims=[0, 2], ylims=[0, 2], ax=None, estimates={}): ax.set_aspect("equal") return fig, ax, im + def plot_cost( + self, + y: np.ndarray | None = None, + xlims: tuple[float, float] | None = None, + ylims: tuple[float, float] | None = None, + grid_size: int = 120, + ) -> tuple[Any, Any, Any]: + """Plot a 2D heatmap of the cost by varying the first position only.""" + if self.d != 2: + raise ValueError("plot_cost is implemented only for d=2.") + + from matplotlib.colors import LogNorm + + if y is None: + if self.y_ is None: + self.y_ = self.simulate_y(noise=0.0) + y = self.y_ + + theta_ref = self.theta.reshape(self.n_positions, self.d).copy() + anchors = self.landmarks + ref_xy = theta_ref[0] + + if xlims is None: + xmin = min(np.min(anchors[:, 0]), ref_xy[0]) - 0.3 + xmax = max(np.max(anchors[:, 0]), ref_xy[0]) + 0.3 + xlims = (xmin, xmax) + if ylims is None: + ymin = min(np.min(anchors[:, 1]), ref_xy[1]) - 0.3 + ymax = max(np.max(anchors[:, 1]), ref_xy[1]) + 0.3 + ylims = (ymin, ymax) + + xs = np.linspace(xlims[0], xlims[1], grid_size) + ys = np.linspace(ylims[0], ylims[1], grid_size) + xx, yy = np.meshgrid(xs, ys) + + zz = np.empty_like(xx) + for i in range(xx.shape[0]): + for j in range(xx.shape[1]): + theta_test = theta_ref.copy() + theta_test[0, :] = np.array([xx[i, j], yy[i, j]]) + zz[i, j] = self.get_cost(theta=theta_test, y=y) + + fig, ax = plt.subplots() + im = ax.pcolormesh( + xx, + yy, + zz + 1e-12, + shading="auto", + norm=LogNorm(vmin=max(1e-12, float(np.min(zz + 1e-12)))), + cmap="viridis", + alpha=0.85, + ) + ax.scatter(*anchors[:, :2].T, color="k", marker="+", label="anchors") + ax.scatter(*theta_ref[:, :2].T, color="C0", marker="o", label="gt") + ax.set_aspect("equal") + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.legend() + fig.colorbar(im, ax=ax, label="cost") + return fig, ax, im + + def plot_setup(self, estimates: dict[str, np.ndarray] | None = None) -> tuple: + """Plot anchors and ground-truth positions, with optional estimates, in 2D or 3D.""" + if estimates is None: + estimates = {} + + theta_gt = self.theta + if np.ndim(theta_gt) < 2: + theta_gt = theta_gt.reshape((self.n_positions, self.d)) + + if self.d == 2: + fig, ax = plt.subplots() + ax.scatter(*self.landmarks[:, :2].T, color="k", marker="+", label="anchors") + ax.scatter(*theta_gt[:, :2].T, color="C0", marker="o", label="gt") + for label, theta_est in estimates.items(): + if np.ndim(theta_est) < 2: + theta_est = theta_est.reshape((self.n_positions, self.d)) + ax.scatter(*theta_est[:, :2].T, marker="x", label=label) + ax.set_aspect("equal") + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.legend() + return fig, ax + elif self.d == 3: + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + ax.scatter(*self.landmarks.T, color="k", marker="+", label="anchors") + ax.scatter(*theta_gt.T, color="C0", marker="o", label="gt") + for label, theta_est in estimates.items(): + if np.ndim(theta_est) < 2: + theta_est = theta_est.reshape((self.n_positions, self.d)) + ax.scatter(*theta_est.T, marker="x", label=label) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + ax.legend() + return fig, ax + else: + raise ValueError("plot_setup is implemented only for d=2 or d=3.") + @property @abstractmethod def size_z(self) -> int: diff --git a/popcor/examples/poly4_lifter.py b/popcor/examples/poly4_lifter.py index ae4edeb..6300074 100644 --- a/popcor/examples/poly4_lifter.py +++ b/popcor/examples/poly4_lifter.py @@ -1,6 +1,7 @@ """Poly4Lifter class for fourth-degree polynomial lifter examples.""" import numpy as np +from poly_matrix import PolyMatrix from popcor.base_lifters import PolyLifter @@ -8,38 +9,52 @@ class Poly4Lifter(PolyLifter): """Fourth-degree polynomial lifter. - Two types are provided: + Two example types are provided and selected with create_example(example_type=...). - - poly_type="A": one global minimum, one local minimum - - poly_type="B": two global minima + - example_type="A": one global minimum at -1, one local minimum at 2. + - example_type="B": two global minima at 1 and 3. """ + EXAMPLE_TYPES: tuple[str, str] = ("A", "B") + example_type: str = "A" + @property def VARIABLE_LIST(self) -> list[list[str]]: return [[self.HOM, "t", "z0"]] - def __init__(self, poly_type: str = "A") -> None: - # actual minimum - assert poly_type in ["A", "B"] - self.poly_type: str = poly_type + def __init__(self) -> None: + self.example_type = "A" super().__init__(degree=4) + @staticmethod + def create_example(example_type: str = "A") -> "Poly4Lifter": + """Create a fourth-degree polynomial example of the requested type.""" + if example_type not in Poly4Lifter.EXAMPLE_TYPES: + raise ValueError( + f"Unknown example_type: {example_type}. Expected one of {Poly4Lifter.EXAMPLE_TYPES}." + ) + + lifter = Poly4Lifter() + lifter.example_type = example_type + lifter.generate_random_setup() + return lifter + def get_Q( self, output_poly: bool = False, noise: float | None = None ) -> np.ndarray: - """Returns the Q matrix for the selected polynomial type.""" + """Return the Q matrix for the selected example type.""" if output_poly: raise ValueError("output_poly not implemented for Poly4Lifter.") - if self.poly_type == "A": + if self.example_type == "A": # Q matrix for type A return np.r_[ np.c_[2, 1, 0], np.c_[1, -1 / 2, -1 / 3], np.c_[0, -1 / 3, 1 / 4] ] - elif self.poly_type == "B": + elif self.example_type == "B": # Q matrix for type B, constructed such that f'(t) = (t-1)*(t-2)*(t-3) return np.r_[np.c_[3, -3, 0], np.c_[-3, 11 / 2, -1], np.c_[0, -1, 1 / 4]] else: - raise ValueError(f"Unknown poly_type: {self.poly_type}") + raise ValueError(f"Unknown example_type: {self.example_type}") def get_A_known( self, @@ -47,8 +62,6 @@ def get_A_known( add_redundant: bool = False, var_dict: dict | None = None, ) -> tuple[list[np.ndarray] | list, list[int]]: - from poly_matrix import PolyMatrix - if add_redundant: print("No redundant constraints for 4-degree polynomial.") @@ -62,7 +75,32 @@ def get_A_known( return [A_1.get_matrix(self.var_dict)], [0] def generate_random_setup(self) -> None: - self.theta_: np.ndarray = np.array([-1]) + if self.example_type == "A": + self.theta_: np.ndarray = np.array([-1]) + elif self.example_type == "B": + self.theta_: np.ndarray = np.array([1]) + else: + raise ValueError(f"Unknown example_type: {self.example_type}") + + def plot_cost( + self, + y: np.ndarray | None = None, + xlims: tuple[float, float] | None = None, + ylims: tuple[float, float] | None = None, + grid_size: int = 120, + ) -> tuple[object, object, object]: + """Plot cost using the same signature as range-only lifters.""" + theta_ref = float(np.asarray(self.theta).reshape(-1)[0]) + if xlims is None: + xlims = (theta_ref - 2.0, theta_ref + 3.0) + + thetas = np.linspace(xlims[0], xlims[1], grid_size) + fig, ax = self.plot(thetas, label="cost", estimates={"theta_ref": theta_ref}) + line = ax.lines[-1] if len(ax.lines) else None + if ylims is not None: + ax.set_ylim(*ylims) + ax.legend() + return fig, ax, line def get_D(self, that: float) -> np.ndarray: D = np.array( @@ -84,14 +122,18 @@ def get_D(self, that: float) -> np.ndarray: base_dir: str = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) thetas: np.ndarray = np.linspace(-2, 3, 100) - poly_lifter: Poly4Lifter = Poly4Lifter(poly_type="A") + poly_lifter: Poly4Lifter = Poly4Lifter.create_example(example_type="A") fig, ax = poly_lifter.plot(thetas) - fig.savefig(os.path.join(base_dir, "docs", "figures", "poly4_lifter_A.png")) + fig.savefig( + os.path.join(base_dir, "docs", "source", "_static", "poly4_lifter_A.png") + ) thetas = np.linspace(0, 4, 100) - poly_lifter = Poly4Lifter(poly_type="B") + poly_lifter = Poly4Lifter.create_example(example_type="B") fig, ax = poly_lifter.plot(thetas) - fig.savefig(os.path.join(base_dir, "docs", "figures", "poly4_lifter_B.png")) + fig.savefig( + os.path.join(base_dir, "docs", "source", "_static", "poly4_lifter_B.png") + ) plt.show(block=False) print("done") diff --git a/popcor/examples/poly6_lifter.py b/popcor/examples/poly6_lifter.py index 265b47c..1c98660 100644 --- a/popcor/examples/poly6_lifter.py +++ b/popcor/examples/poly6_lifter.py @@ -1,7 +1,6 @@ """Poly6Lifter class for sixth-degree polynomial examples.""" import numpy as np -import scipy.sparse as sp from poly_matrix import PolyMatrix from popcor.base_lifters import PolyLifter @@ -10,28 +9,43 @@ class Poly6Lifter(PolyLifter): """Sixth-degree polynomial examples. - Two types are provided: + Two example types are provided and selected with create_example(example_type=...). - - poly_type="A": one global minimum, two local minima, two local maxima - - poly_type="B": one global minimum, one local minimum, one local maximum + - example_type="A": one global minimum, two local minima, two local maxima + - example_type="B": one global minimum, one local minimum, one local maximum """ + EXAMPLE_TYPES: tuple[str, str] = ("A", "B") + example_type: str = "A" + @property def VARIABLE_LIST(self) -> list[list[str]]: return [[self.HOM, "t", "z0", "z1"]] - def __init__(self, poly_type: str = "A") -> None: - assert poly_type in ["A", "B"] - self.poly_type: str = poly_type + def __init__(self) -> None: + self.example_type = "A" super().__init__(degree=6) + @staticmethod + def create_example(example_type: str = "A") -> "Poly6Lifter": + """Create a sixth-degree polynomial example of the requested type.""" + if example_type not in Poly6Lifter.EXAMPLE_TYPES: + raise ValueError( + f"Unknown example_type: {example_type}. Expected one of {Poly6Lifter.EXAMPLE_TYPES}." + ) + + lifter = Poly6Lifter() + lifter.example_type = example_type + lifter.generate_random_setup() + return lifter + def get_Q( self, output_poly: bool = False, noise: float | None = None ) -> np.ndarray: - """Returns the Q matrix for the selected polynomial type.""" + """Return the Q matrix for the selected example type.""" if output_poly: raise ValueError("output_poly not implemented for Poly6Lifter.") - if self.poly_type == "A": + if self.example_type == "A": return 0.1 * np.array( [ [25, 12, 0, 0], @@ -40,7 +54,7 @@ def get_Q( [0, 0, -0.9, 1 / 6], ] ) - elif self.poly_type == "B": + elif self.example_type == "B": return np.array( [ [5.0000, 1.3167, -1.4481, 0], @@ -50,7 +64,7 @@ def get_Q( ] ) else: - raise ValueError(f"Unknown poly_type: {self.poly_type}") + raise ValueError(f"Unknown example_type: {self.example_type}") def get_A_known( self, @@ -59,8 +73,6 @@ def get_A_known( var_dict: dict | None = None, ) -> tuple[list, list[float]]: """Returns the list of known A matrices and their corresponding values.""" - from poly_matrix import PolyMatrix - A_list: list = [] # z_0 = t^2 @@ -101,6 +113,26 @@ def get_D(self, that: float) -> np.ndarray: ) return D + def plot_cost( + self, + y: np.ndarray | None = None, + xlims: tuple[float, float] | None = None, + ylims: tuple[float, float] | None = None, + grid_size: int = 120, + ) -> tuple[object, object, object]: + """Plot cost using the same signature as range-only lifters.""" + theta_ref = float(np.asarray(self.theta).reshape(-1)[0]) + if xlims is None: + xlims = (theta_ref - 3.0, theta_ref + 4.0) + + thetas = np.linspace(xlims[0], xlims[1], grid_size) + fig, ax = self.plot(thetas, label="cost", estimates={"theta_ref": theta_ref}) + line = ax.lines[-1] if len(ax.lines) else None + if ylims is not None: + ax.set_ylim(*ylims) + ax.legend() + return fig, ax, line + def generate_random_setup(self) -> None: """Initializes a random setup for theta_.""" self.theta_ = np.array([-1]) @@ -115,14 +147,18 @@ def generate_random_setup(self) -> None: base_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) thetas = np.linspace(-1.5, 4.5, 100) - poly_lifter = Poly6Lifter(poly_type="A") + poly_lifter = Poly6Lifter.create_example(example_type="A") fig, ax = poly_lifter.plot(thetas) - fig.savefig(os.path.join(base_dir, "docs", "figures", "poly6_lifter_A.png")) + fig.savefig( + os.path.join(base_dir, "docs", "source", "_static", "poly6_lifter_A.png") + ) thetas = np.linspace(-3, 3, 100) - poly_lifter = Poly6Lifter(poly_type="B") + poly_lifter = Poly6Lifter.create_example(example_type="B") fig, ax = poly_lifter.plot(thetas) - fig.savefig(os.path.join(base_dir, "docs", "figures", "poly6_lifter_B.png")) + fig.savefig( + os.path.join(base_dir, "docs", "source", "_static", "poly6_lifter_B.png") + ) plt.show(block=False) print("done") diff --git a/popcor/examples/ro_nsq_lifter.py b/popcor/examples/ro_nsq_lifter.py index 7802020..5bfa69c 100644 --- a/popcor/examples/ro_nsq_lifter.py +++ b/popcor/examples/ro_nsq_lifter.py @@ -43,6 +43,8 @@ class RangeOnlyNsqLifter(RangeOnlyLifter): "normals": "$\\boldsymbol{y}_n$", "simple": "$z_n$", } + EXAMPLE_TYPES: tuple[str, str, str] = ("A", "B", "C") + example_type: str | None = None MONOMIAL_DEGREE: int = 1 SCALE: float = 1.0 @@ -51,33 +53,21 @@ def create_good( n_positions: int, n_landmarks: int, d: int = 2 ) -> "RangeOnlyNsqLifter": """Create a lifter with good initial positions.""" - landmarks, theta = RangeOnlyLifter.create_good(n_positions, n_landmarks, d) - lifter = RangeOnlyNsqLifter(n_positions, n_landmarks, d) - lifter.overwrite_theta(theta) - lifter.landmarks = landmarks - return lifter + return RangeOnlyNsqLifter.create_example("A", n_positions, n_landmarks, d) @staticmethod def create_bad( n_positions: int, n_landmarks: int, d: int = 2 ) -> "RangeOnlyNsqLifter": """Create a lifter with bad initial positions.""" - landmarks, theta = RangeOnlyLifter.create_bad(n_positions, n_landmarks, d) - lifter = RangeOnlyNsqLifter(n_positions, n_landmarks, d) - lifter.overwrite_theta(theta) - lifter.landmarks = landmarks - return lifter + return RangeOnlyNsqLifter.create_example("B", n_positions, n_landmarks, d) @staticmethod def create_bad_fixed( n_positions: int, n_landmarks: int, d: int = 2 ) -> "RangeOnlyNsqLifter": """Create a lifter with fixed bad initial positions.""" - landmarks, theta = RangeOnlyLifter.create_bad_fixed(n_positions, n_landmarks, d) - lifter = RangeOnlyNsqLifter(n_positions, n_landmarks, d) - lifter.overwrite_theta(theta) - lifter.landmarks = landmarks - return lifter + return RangeOnlyNsqLifter.create_example("C", n_positions, n_landmarks, d) def __init__( self, @@ -124,7 +114,7 @@ def get_A_known( var_dict: dict | None = None, output_poly: bool = False, add_redundant: bool = False, - ) -> list[np.ndarray | PolyMatrix]: + ) -> tuple[list[np.ndarray | PolyMatrix], list[float]]: """Return known quadratic constraints for the lifted variables.""" if var_dict is None: var_dict = self.var_dict @@ -147,7 +137,7 @@ def get_A_known( raise NotImplementedError( "get_A_known not implemented yet for simple level" ) - return A_list + return A_list, [0.0] * len(A_list) def get_residuals( self, t: np.ndarray, y: np.ndarray, ad: bool = False @@ -296,4 +286,54 @@ def __repr__(self) -> str: if __name__ == "__main__": - lifter = RangeOnlyNsqLifter(n_positions=3, n_landmarks=4, d=2) + import os + + import matplotlib.pyplot as plt + from cert_tools.linalg_tools import rank_project + from cert_tools.sdp_solvers import solve_sdp + + np.random.seed(0) + base_dir: str = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) + + for example_type in RangeOnlyNsqLifter.EXAMPLE_TYPES: + lifter = RangeOnlyNsqLifter.create_example( + example_type=example_type, + n_positions=1, + n_landmarks=4, + d=2, + ) + y = lifter.simulate_y(noise=0.0) + + Q = lifter.get_Q_from_y(y=y, output_poly=False) + A_known, b_known = lifter.get_A_known(output_poly=False, add_redundant=True) + A_0, b_0 = lifter.get_A0() + constraints = list(zip(A_known + A_0, b_known + b_0)) + + X, _ = solve_sdp(Q, constraints, verbose=False) + x_round, _ = rank_project(X, p=1) + x_round = x_round.flatten() + x_round = x_round / x_round[0] + theta_est = x_round[1 : 1 + lifter.N] + + fig, ax, _ = lifter.plot_cost(y=y) + ax.set_title(f"RangeOnlyNsqLifter cost ({example_type})") + fig.savefig( + os.path.join( + base_dir, "docs", "source", "_static", f"ro_nsq_cost_{example_type}.png" + ) + ) + + fig, ax = lifter.plot_setup(estimates={"estimate": theta_est}) + ax.set_title(f"RangeOnlyNsqLifter setup ({example_type})") + fig.savefig( + os.path.join( + base_dir, + "docs", + "source", + "_static", + f"ro_nsq_setup_{example_type}.png", + ) + ) + + plt.show(block=False) + print("done") diff --git a/popcor/examples/ro_sq_lifter.py b/popcor/examples/ro_sq_lifter.py index e101c1b..dcc0540 100644 --- a/popcor/examples/ro_sq_lifter.py +++ b/popcor/examples/ro_sq_lifter.py @@ -44,37 +44,27 @@ class RangeOnlySqLifter(RangeOnlyLifter): "no": "$z_n$", "quad": "$\\boldsymbol{y}_n$", } + EXAMPLE_TYPES: tuple[str, str, str] = ("A", "B", "C") + example_type: str | None = None MONOMIAL_DEGREE = 2 @staticmethod def create_good( n_positions: int, n_landmarks: int, d: int = 2 ) -> "RangeOnlySqLifter": - landmarks, theta = RangeOnlyLifter.create_good(n_positions, n_landmarks, d) - lifter = RangeOnlySqLifter(n_positions, n_landmarks, d) - lifter.overwrite_theta(theta) - lifter.landmarks = landmarks - return lifter + return RangeOnlySqLifter.create_example("A", n_positions, n_landmarks, d) @staticmethod def create_bad( n_positions: int, n_landmarks: int, d: int = 2 ) -> "RangeOnlySqLifter": - landmarks, theta = RangeOnlyLifter.create_bad(n_positions, n_landmarks, d) - lifter = RangeOnlySqLifter(n_positions, n_landmarks, d) - lifter.overwrite_theta(theta) - lifter.landmarks = landmarks - return lifter + return RangeOnlySqLifter.create_example("B", n_positions, n_landmarks, d) @staticmethod def create_bad_fixed( n_positions: int, n_landmarks: int, d: int = 2 ) -> "RangeOnlySqLifter": - landmarks, theta = RangeOnlyLifter.create_bad_fixed(n_positions, n_landmarks, d) - lifter = RangeOnlySqLifter(n_positions, n_landmarks, d) - lifter.overwrite_theta(theta) - lifter.landmarks = landmarks - return lifter + return RangeOnlySqLifter.create_example("C", n_positions, n_landmarks, d) def __init__( self, @@ -353,14 +343,14 @@ def get_grad( sub_idx_x = self.get_sub_idx_x(sub_idx) return 2 * J.T[:, sub_idx_x] @ Q[sub_idx_x, :][:, sub_idx_x] @ x[sub_idx_x] # type: ignore - def get_J(self, t: np.ndarray, y: np.ndarray) -> sp.csr_matrix: + def get_J(self, t: np.ndarray, y: np.ndarray) -> sp.csr_array: J = sp.csr_matrix( (np.ones(self.N), (range(1, self.N + 1), range(self.N))), shape=(self.N + 1, self.N), ) J_lift = self.get_J_lifting(t) J = sp.vstack([J, J_lift]) - assert isinstance(J, sp.csr_matrix) + assert isinstance(J, sp.csr_array) return J def get_hess(self, t: np.ndarray, y: np.ndarray) -> np.ndarray: @@ -391,4 +381,50 @@ def get_D(self, that: np.ndarray) -> sp.csc_array: if __name__ == "__main__": - lifter = RangeOnlySqLifter(n_positions=3, n_landmarks=4, d=2) + import os + + import matplotlib.pyplot as plt + from cert_tools.linalg_tools import rank_project + from cert_tools.sdp_solvers import solve_sdp + + np.random.seed(0) + base_dir: str = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) + + for example_type in RangeOnlySqLifter.EXAMPLE_TYPES: + lifter = RangeOnlySqLifter.create_example( + example_type=example_type, + n_positions=1, + n_landmarks=4, + d=2, + ) + y = lifter.simulate_y(noise=0.0) + + Q = lifter.get_Q_from_y(y=y, output_poly=False) + A_known, b_known = lifter.get_A_known(output_poly=False, add_redundant=True) + A_0, b_0 = lifter.get_A0() + constraints = list(zip(A_known + A_0, b_known + b_0)) + + X, _ = solve_sdp(Q, constraints, verbose=False) + x_round, _ = rank_project(X, p=1) + x_round = x_round.flatten() + x_round = x_round / x_round[0] + theta_est = x_round[1 : 1 + lifter.N] + + fig, ax, _ = lifter.plot_cost(y=y) + ax.set_title(f"RangeOnlySqLifter cost ({example_type})") + fig.savefig( + os.path.join( + base_dir, "docs", "source", "_static", f"ro_sq_cost_{example_type}.png" + ) + ) + + fig, ax = lifter.plot_setup(estimates={"estimate": theta_est}) + ax.set_title(f"RangeOnlySqLifter setup ({example_type})") + fig.savefig( + os.path.join( + base_dir, "docs", "source", "_static", f"ro_sq_setup_{example_type}.png" + ) + ) + + plt.show(block=False) + print("done") diff --git a/popcor/examples/rotation_lifter.py b/popcor/examples/rotation_lifter.py index 3f40f76..b2ede02 100644 --- a/popcor/examples/rotation_lifter.py +++ b/popcor/examples/rotation_lifter.py @@ -5,7 +5,7 @@ import numpy as np import scipy.sparse as sp -from poly_matrix.poly_matrix import PolyMatrix +from poly_matrix import PolyMatrix from scipy.spatial.transform import Rotation from popcor.base_lifters import StateLifter @@ -106,7 +106,8 @@ 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") + EXAMPLE_TYPES: Tuple[str, str] = ("A", "B") + example_type: str | None = None ADD_DETERMINANT: bool = False NOISE: float = 1e-3 @@ -136,6 +137,18 @@ def __init__( d=d, ) + @staticmethod + def create_example(example_type: str = "A") -> "RotationLifter": + """Create a deterministic SO(2) tutorial example of the requested type.""" + if example_type not in RotationLifter.EXAMPLE_TYPES: + raise ValueError( + f"Unknown example_type {example_type}. Expected one of {RotationLifter.EXAMPLE_TYPES}." + ) + + lifter = RotationLifter(d=2, n_rot=1, n_abs=0, n_rel=0, level="no") + lifter.example_type = example_type + return lifter + @property def var_dict(self) -> Dict[str, int]: """Return dictionary mapping variable names to their dimensions in the lifted representation.""" @@ -167,73 +180,6 @@ def so2_theta(angle: float) -> np.ndarray: 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, @@ -360,6 +306,39 @@ def get_Q( self, noise: float | None = None, output_poly: bool = False ) -> PolyMatrix | np.ndarray | sp.csr_matrix | sp.csc_matrix: """Return the cost matrix Q (poly or ndarray) constructed from simulated measurements.""" + if getattr(self, "example_type", None) is not None: + 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'." + ) + + if self.example_type == "A": + q2 = np.array([[0.8, 0.05], [0.05, 0.2]]) + q1 = np.array([0.0, 0.0]) + elif self.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 {self.example_type}. Expected one of {self.EXAMPLE_TYPES}." + ) + + 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) + if noise is None: noise = self.NOISE if self.y_ is None: @@ -544,27 +523,49 @@ def get_A_known( self.test_and_add(A_list, Ai, output_poly, b_list, 0.0) return A_list, b_list - def plot(self, estimates: Dict[str, np.ndarray] = {}) -> Tuple[Any, Any]: - """Plot ground-truth frames and optional estimated frames; returns (fig, ax).""" + def plot_cost( + self, thetas: np.ndarray, label: str | None = None + ) -> Tuple[Any, Any]: + """Plot the cost profile as a polar plot for d=2 single-rotation examples.""" + import matplotlib.pyplot as plt + + if not (self.d == 2 and self.n_rot == 1 and self.level == "no"): + raise ValueError( + "Polar cost plotting is implemented only for d=2, n_rot=1, level='no'." + ) + costs = np.array([self.get_cost(self.so2_theta(angle)) for angle in thetas]) + fig, ax = plt.subplots(subplot_kw={"projection": "polar"}) + ax.plot(thetas, costs, label=label) + if label is not None: + ax.legend() + return fig, ax + + def plot_setup( + self, estimates: Dict[str, np.ndarray] | None = None + ) -> Tuple[Any, Any]: + """Plot ground-truth frames and optional estimated frames.""" import itertools import matplotlib.pyplot as plt from popcor.utils.plotting_tools import plot_frame + if estimates is None: + estimates = {} + fig, ax = plt.subplots() - label = "gt" + gt_label = "gt" for i in range(self.n_rot): plot_frame( ax=ax, theta=self.theta[:, i * self.d : (i + 1) * self.d], - label=label, + label=gt_label, ls="-", scale=0.5, marker="", r_wc_w=np.hstack([i * 2.0] + [0.0] * (self.d - 1)), # type: ignore ) - label = None + gt_label = None linestyles = itertools.cycle(LINESTYLES) for label, theta in estimates.items(): @@ -585,18 +586,48 @@ def plot(self, estimates: Dict[str, np.ndarray] = {}) -> Tuple[Any, Any]: ax.legend() return fig, ax + def plot( + self, + thetas: np.ndarray | None = None, + label: str | None = None, + estimates: Dict[str, np.ndarray] | None = None, + ) -> Tuple[Any, Any]: + """Compatibility wrapper for plot_cost/plot_setup.""" + if thetas is not None: + return self.plot_cost(thetas=thetas, label=label) + return self.plot_setup(estimates=estimates) + def __repr__(self) -> str: return f"rotation_lifter{self.d}d_{self.level}" if __name__ == "__main__": + import os + import matplotlib.pyplot as plt - import numpy as np from cert_tools.linalg_tools import rank_project from cert_tools.sdp_solvers import solve_sdp from popcor.utils.plotting_tools import plot_matrix + base_dir: str = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) + + angles: np.ndarray = np.linspace(-np.pi, np.pi, 500) + for example_type in RotationLifter.EXAMPLE_TYPES: + lifter = RotationLifter.create_example(example_type=example_type) + fig, ax = lifter.plot_cost(angles, label=f"example {example_type}") + ax.set_title(f"SO(2) example {example_type}") + fig.savefig( + os.path.join( + base_dir, + "docs", + "source", + "_static", + f"rotation_lifter_{example_type}.png", + ) + ) + + # Create estimation pipeline level: str = "no" np.random.seed(0) lifter = RotationLifter( @@ -615,7 +646,7 @@ def __repr__(self) -> str: theta_i, *_ = lifter.local_solver(theta_init, y, verbose=False) estimates[f"init random {i}"] = theta_i - fig, ax = lifter.plot(estimates=estimates) + fig, ax = lifter.plot_setup(estimates=estimates) ax.legend() plt.show(block=False) @@ -644,6 +675,6 @@ def __repr__(self) -> str: theta_opt = lifter.get_theta(x[:, :rank]) estimates = {"init gt": theta_gt, "SDP": theta_opt} - fig, ax = lifter.plot(estimates=estimates) + fig, ax = lifter.plot_setup(estimates=estimates) print("done") diff --git a/setup.cfg b/setup.cfg index 8f55669..b2a2b01 100644 --- a/setup.cfg +++ b/setup.cfg @@ -25,9 +25,8 @@ install_requires = chompack>=2.3 autograd asrl-pylgmath>=1.0.3 - sparseqr poly_matrix @ git+https://github.com/utiasASRL/poly_matrix.git@v0.3.1#egg=poly_matrix - cert_tools @ git+https://github.com/utiasASRL/certifiable-tools.git@v0.0.5#egg=cert_tools + cert_tools @ git+https://github.com/utiasASRL/certifiable-tools.git@v0.0.7#egg=cert_tools [flake8] ignore = W292, W391, F541, F841, E203, E501, W503, E741 diff --git a/tests/test_quickstart.py b/tests/test_quickstart.py index 74aa3b9..847484e 100644 --- a/tests/test_quickstart.py +++ b/tests/test_quickstart.py @@ -8,7 +8,7 @@ def test_setup_problem(): """Test the problem setup example from the documentation.""" from popcor.examples import Poly4Lifter - lifter = Poly4Lifter() + lifter = Poly4Lifter.create_example(example_type="A") Q = lifter.get_Q() @@ -30,7 +30,7 @@ def test_solve_sdp(): from popcor.examples import Poly4Lifter - lifter = Poly4Lifter() + lifter = Poly4Lifter.create_example(example_type="A") # the cost matrix Q = lifter.get_Q() diff --git a/tests/test_rotation_lifter.py b/tests/test_rotation_lifter.py index 2ff3278..285483e 100644 --- a/tests/test_rotation_lifter.py +++ b/tests/test_rotation_lifter.py @@ -144,10 +144,10 @@ def test_solve_local(): ) 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") + lifter = RotationLifter.create_example(example_type=example_type) 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]) + costs = np.array([lifter.get_cost(lifter.so2_theta(t)) for t in angles]) minima_idx = [] for i in range(1, len(angles) - 1): @@ -175,9 +175,9 @@ def test_so2_example_sdp_recovery(example_type): 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") + lifter = RotationLifter.create_example(example_type=example_type) - Q = lifter.get_so2_example_Q(example_type=example_type) + Q = lifter.get_Q() 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)) @@ -186,7 +186,7 @@ def test_so2_example_sdp_recovery(example_type): _, 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]) + costs = np.array([lifter.get_cost(lifter.so2_theta(t)) for t in angles]) c_min = float(np.min(costs)) # SDP lower bound should match the global minimum cost very tightly.