Implement object-oriented uncertainty blocks in tests and examples

timeverettadams · timeverettadams · commit dd48ded98124 · 2025-05-14T14:05:51.000-04:00
diff --git a/examples/1_example_dk_iter_fixed_order.py b/examples/1_example_dk_iter_fixed_order.py
@@ -19,7 +19,12 @@ def example_dk_iter_fixed_order():
     )
 
     omega = np.logspace(-3, 3, 61)
-    block_structure = np.array([[1, 1], [1, 1], [2, 2]])
+    # block_structure = np.array([[1, 1], [1, 1], [2, 2]])
+    block_structure = [
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(2, 2),
+    ]
     K, N, mu, iter_results, info = dk_iter.synthesize(
         eg["P"],
         eg["n_y"],
diff --git a/examples/2_example_dk_iter_list_order.py b/examples/2_example_dk_iter_list_order.py
@@ -23,7 +23,11 @@ def example_dk_iter_list_order():
     )
 
     omega = np.logspace(-3, 3, 61)
-    block_structure = np.array([[1, 1], [1, 1], [2, 2]])
+    block_structure = [
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(2, 2),
+    ]
     K, N, mu, iter_results, info = dk_iter.synthesize(
         eg["P"],
         eg["n_y"],
diff --git a/examples/3_example_dk_iter_auto_order.py b/examples/3_example_dk_iter_auto_order.py
@@ -45,6 +45,11 @@ def example_dk_iter_auto_order():
 
     omega = np.logspace(-3, 3, 61)
     block_structure = np.array([[1, 1], [1, 1], [2, 2]])
+    block_structure = [
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(2, 2),
+    ]
     K, N, mu, iter_results, info = dk_iter.synthesize(
         eg["P"],
         eg["n_y"],
diff --git a/examples/4_example_dk_iter_interactive.py b/examples/4_example_dk_iter_interactive.py
@@ -37,7 +37,11 @@ def example_dk_iter_interactive():
     )
 
     omega = np.logspace(-3, 3, 61)
-    block_structure = np.array([[1, 1], [1, 1], [2, 2]])
+    block_structure = [
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(2, 2),
+    ]
     K, N, mu, iter_results, info = dk_iter.synthesize(
         eg["P"],
         eg["n_y"],
diff --git a/examples/5_example_dk_iteration_custom.py b/examples/5_example_dk_iteration_custom.py
@@ -70,7 +70,11 @@ def example_dk_iter_custom():
     )
 
     omega = np.logspace(-3, 3, 61)
-    block_structure = np.array([[1, 1], [1, 1], [2, 2]])
+    block_structure = [
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(1, 1),
+        dkpy.ComplexFullBlock(2, 2),
+    ]
     K, N, mu, iter_results, info = dk_iter.synthesize(
         eg["P"],
         eg["n_y"],
diff --git a/src/dkpy/d_scale_fit.py b/src/dkpy/d_scale_fit.py
@@ -6,7 +6,7 @@
 ]
 
 import abc
-from typing import Optional, Tuple, Union, List
+from typing import Optional, Tuple, Union, List, Sequence
 import warnings
 
 import control
@@ -31,7 +31,9 @@ def fit(
         omega: np.ndarray,
         D_omega: np.ndarray,
         order: Union[int, np.ndarray] = 0,
-        block_structure: Optional[np.ndarray] = None,
+        block_structure: Optional[
+            Sequence[Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]]
+        ] = None,
     ) -> Tuple[control.StateSpace, control.StateSpace]:
         """Fit D-scale magnitudes.
 
@@ -44,7 +46,7 @@ def fit(
             dimension.
         order : Union[int, np.ndarray]
             Transfer function order to fit. Can be specified per-entry.
-        block_structure : np.ndarray
+        block_structure : Sequence[RealDiagonalBlock | ComplexDiagonalBlock | ComplexFullBlock]
             2D array with 2 columns and as many rows as uncertainty blocks
             in Delta. The columns represent the number of rows and columns in
             each uncertainty block. See [#mussv]_.
@@ -112,7 +114,7 @@ def fit(
         D_omega: np.ndarray,
         order: Union[int, np.ndarray] = 0,
         block_structure: Optional[
-            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+            List[Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]]
         ] = None,
     ) -> Tuple[control.StateSpace, control.StateSpace]:
         # Get mask
@@ -192,12 +194,15 @@ def _mask_from_block_structure(
     num_blocks = len(block_structure)
     X_lst = []
     for i in range(num_blocks):
+        # Uncertainty block
         block = block_structure[i]
+        # Square uncertainty block condition
+        is_block_square = block.num_inputs == block.num_outputs
         if isinstance(block, RealDiagonalBlock):
             raise NotImplementedError("Real perturbations are not yet supported.")
         if isinstance(block, ComplexDiagonalBlock):
             raise NotImplementedError("Diagonal perturbations are not yet supported.")
-        if isinstance(block, ComplexFullBlock) and (not block.is_square):
+        if isinstance(block, ComplexFullBlock) and (not is_block_square):
             raise NotImplementedError("Nonsquare perturbations are not yet supported.")
         # Set last scaling to identity
         if i == num_blocks - 1:
diff --git a/src/dkpy/dk_iteration.py b/src/dkpy/dk_iteration.py
@@ -13,7 +13,7 @@
 
 import abc
 import logging
-from typing import Any, Dict, List, Optional, Tuple, Union
+from typing import Any, Dict, List, Optional, Tuple, Union, Sequence
 
 import control
 import numpy as np
@@ -26,6 +26,11 @@
     structured_singular_value,
     utilities,
 )
+from .uncertainty_structure import (
+    RealDiagonalBlock,
+    ComplexDiagonalBlock,
+    ComplexFullBlock,
+)
 
 
 class IterResult:
@@ -46,7 +51,9 @@ def __init__(
         mu_fit_omega: np.ndarray,
         D_fit_omega: np.ndarray,
         D_fit: control.StateSpace,
-        block_structure: np.ndarray,
+        block_structure: List[
+            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+        ],
     ):
         """Instantiate :class:`IterResult`.
 
@@ -91,7 +98,9 @@ def create_from_fit(
         K: control.StateSpace,
         D_fit: control.StateSpace,
         D_fit_inv: control.StateSpace,
-        block_structure: np.ndarray,
+        block_structure: List[
+            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+        ],
     ) -> "IterResult":
         """Instantiate :class:`IterResult` from fit D-scales.
 
@@ -204,7 +213,9 @@ def synthesize(
         n_y: int,
         n_u: int,
         omega: np.ndarray,
-        block_structure: np.ndarray,
+        block_structure: Sequence[
+            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+        ],
     ) -> Tuple[
         control.StateSpace,
         control.StateSpace,
@@ -341,7 +352,9 @@ def _get_fit_order(
         D_omega: np.ndarray,
         P: control.StateSpace,
         K: control.StateSpace,
-        block_structure: np.ndarray,
+        block_structure: List[
+            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+        ],
     ) -> Optional[Union[int, np.ndarray]]:
         """Get D-scale fit order.
 
@@ -440,7 +453,9 @@ def _get_fit_order(
         D_omega: np.ndarray,
         P: control.StateSpace,
         K: control.StateSpace,
-        block_structure: np.ndarray,
+        block_structure: List[
+            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+        ],
     ) -> Optional[Union[int, np.ndarray]]:
         if iteration < self.n_iterations:
             return self.fit_order
@@ -507,7 +522,9 @@ def _get_fit_order(
         D_omega: np.ndarray,
         P: control.StateSpace,
         K: control.StateSpace,
-        block_structure: np.ndarray,
+        block_structure: List[
+            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+        ],
     ) -> Optional[Union[int, np.ndarray]]:
         if iteration < len(self.fit_orders):
             return self.fit_orders[iteration]
@@ -589,7 +606,9 @@ def _get_fit_order(
         D_omega: np.ndarray,
         P: control.StateSpace,
         K: control.StateSpace,
-        block_structure: np.ndarray,
+        block_structure: List[
+            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+        ],
     ) -> Optional[Union[int, np.ndarray]]:
         # Check termination conditions
         if (self.max_iterations is not None) and (iteration >= self.max_iterations):
@@ -686,14 +705,16 @@ def __init__(
 
     def _get_fit_order(
         self,
-        iteration,
-        omega,
-        mu_omega,
-        D_omega,
-        P,
-        K,
-        block_structure,
-    ):
+        iteration: int,
+        omega: np.ndarray,
+        mu_omega: np.ndarray,
+        D_omega: np.ndarray,
+        P: control.StateSpace,
+        K: control.StateSpace,
+        block_structure: List[
+            Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
+        ],
+    ) -> Optional[Union[int, np.ndarray]]:
         d_info = []
         for fit_order in range(self.max_fit_order + 1):
             D_fit, D_fit_inv = self.d_scale_fit.fit(
diff --git a/src/dkpy/structured_singular_value.py b/src/dkpy/structured_singular_value.py
@@ -7,7 +7,7 @@
 
 import abc
 import warnings
-from typing import Any, Dict, Optional, Tuple, Union, List
+from typing import Any, Dict, Optional, Tuple, Union, List, Sequence
 
 import cvxpy
 import joblib
@@ -29,7 +29,7 @@ class StructuredSingularValue(metaclass=abc.ABCMeta):
     def compute_ssv(
         self,
         N_omega: np.ndarray,
-        block_structure: List[
+        block_structure: Sequence[
             Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
         ],
     ) -> Tuple[np.ndarray, np.ndarray, Dict[str, Any]]:
@@ -39,22 +39,15 @@ def compute_ssv(
         ----------
         N_omega : np.ndarray
             Closed-loop transfer function evaluated at each frequency.
-        block_structure : np.ndarray
-            2D array with 2 columns and as many rows as uncertainty blocks
-            in Delta. The columns represent the number of rows and columns in
-            each uncertainty block. See [#mussv]_.
-
+        block_structure : Sequence[RealDiagonalBlock | ComplexDiagonalBlock | ComplexFullBlock]
+            Sequence of uncertainty block objects.
         Returns
         -------
         Tuple[np.ndarray, np.ndarray, Dict[str, Any]]
             Structured singular value at each frequency, D-scales at each
             frequency, and solution information. If the structured singular
             value cannot be computed, the first two elements of the tuple are
             ``None``, but solution information is still returned.
-
-        References
-        ----------
-        .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html
         """
         raise NotImplementedError()
 
@@ -146,7 +139,7 @@ def __init__(
     def compute_ssv(
         self,
         N_omega: np.ndarray,
-        block_structure: List[
+        block_structure: Sequence[
             Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
         ],
     ) -> Tuple[np.ndarray, np.ndarray, Dict[str, Any]]:
@@ -329,35 +322,32 @@ def _ssv_at_omega(
 
 
 def _variable_from_block_structure(
-    block_structure: List[
+    block_structure: Sequence[
         Union[RealDiagonalBlock, ComplexDiagonalBlock, ComplexFullBlock]
     ],
 ) -> cvxpy.Variable:
     """Get optimization variable with specified block structure.
 
     Parameters
     ----------
-    block_structure : np.ndarray
-        2D array with 2 columns and as many rows as uncertainty blocks
-        in Delta. The columns represent the number of rows and columns in
-        each uncertainty block. See [#mussv]_.
+    block_structure : Sequence[RealDiagonalBlock | ComplexDiagonalBlock | ComplexFullBlock]
+        Sequence of uncertainty block objects.
 
     Returns
     -------
     cvxpy.Variable
         CVXPY variable with specified block structure.
-
-    References
-    ----------
-    .. [#mussv] https://www.mathworks.com/help/robust/ref/mussv.html
     """
     num_blocks = len(block_structure)
     X_lst = []
     for i in range(num_blocks):
         row = []
         for j in range(num_blocks):
+            # Uncertainty blocks
             block_i = block_structure[i]
             block_j = block_structure[j]
+            # Square uncertainty block condition
+            is_block_square = block_i.num_inputs == block_i.num_outputs
             if i == j:
                 # If on the block diagonal, insert variable
                 if isinstance(block_i, RealDiagonalBlock):
@@ -368,7 +358,7 @@ def _variable_from_block_structure(
                     raise NotImplementedError(
                         "Complex diagonal perturbations are not yet supported."
                     )
-                if isinstance(block_i, ComplexFullBlock) and (not block_i.is_square):
+                if isinstance(block_i, ComplexFullBlock) and (not is_block_square):
                     raise NotImplementedError(
                         "Nonsquare perturbations are not yet supported."
                     )
diff --git a/src/dkpy/uncertainty_structure.py b/src/dkpy/uncertainty_structure.py
@@ -19,7 +19,6 @@ def __init__(self, num_channels: int):
         # Uncertainty block parameters
         self.num_inputs = num_channels
         self.num_outputs = num_channels
-        self.is_square = True
 
 
 class ComplexDiagonalBlock:
@@ -40,7 +39,6 @@ def __init__(self, num_channels: int):
         # Uncertainty block parameters
         self.num_inputs = num_channels
         self.num_outputs = num_channels
-        self.is_square = True
 
 
 class ComplexFullBlock:
@@ -63,4 +61,3 @@ def __init__(self, num_inputs: int, num_outputs: int):
         # Uncertainty block parameters
         self.num_inputs = num_inputs
         self.num_outputs = num_outputs
-        self.is_square = self.num_inputs == self.num_outputs
diff --git a/tests/test_d_scale_fit.py b/tests/test_d_scale_fit.py
