Add bounded and inequality operand definition (#53)

drpaprika · web-flow · commit df4ea02bea6a · 2025-01-31T17:00:15.000+01:00
* Add bounded and inequality operand definition

* fix: tests

* update notebook 7d

* update notebook 7d

* simplify Operand : remove bounds, rename min_val, max_val

* fix: tests

* fix: convert one-element array to float properly

* add tests for intercept_lcs operands

* simplify delta_ineq() logic

* Add comments to _fun()

* fix: check that decenter is not a single element array
diff --git a/docs/examples/Tutorial_7d_Three_Mirror_Anastigmat.ipynb b/docs/examples/Tutorial_7d_Three_Mirror_Anastigmat.ipynb
diff --git a/optiland/optimization/operand/operand.py b/optiland/optimization/operand/operand.py
@@ -15,6 +15,7 @@
 
 Kramer Harrison, 2024
 """
+from dataclasses import dataclass
 from optiland.optimization.operand.paraxial import ParaxialOperand
 from optiland.optimization.operand.aberration import AberrationOperand
 from optiland.optimization.operand.ray import RayOperand
@@ -137,47 +138,79 @@ def __repr__(self):
 for name, func in METRIC_DICT.items():
     operand_registry.register(name, func)
 
-
-class Operand(object):
+@dataclass
+class Operand:
     """
     Represents an operand used in optimization calculations.
+    If no target is specified, a default is created at the current value.
 
     Attributes:
-        type (str): The type of the operand.
-        target (float): The target value for the operand.
+        operand_type (str): The type of the operand.
+        target (float): The target value of the operand (equality operand).
+        min_val (float): The operand should stay above this value (inequality operand).
+        max_val (float): The operand should stay below this value (inequality operand).
         weight (float): The weight of the operand.
-        input_data (dict): Additional input data for the operand's metric
-            function.
+        input_data (dict): Additional input data for the operand.
 
     Methods:
         value(): Get the current value of the operand.
-        delta(): Calculate the difference between the target and current value.
+        delta_target(): Calculate the difference between the value and target.
+        delta_ineq(): Calculate the difference between the value and targets.
         fun(): Calculate the objective function value.
     """
 
-    def __init__(self, operand_type, target, weight, input_data={}):
-        self.type = operand_type
-        self.target = target
-        self.weight = weight
-        self.input_data = input_data
+    operand_type: str = None
+    target: float = None
+    min_val: float = None
+    max_val: float = None
+    weight: float = None
+    input_data: dict = None
+
+    def __post_init__(self):
+        if self.min_val is not None and self.max_val is not None and self.min_val > self.max_val:
+            raise ValueError(f"{self.operand_type} operand: min_val is higher than max_val")
+        if self.target is not None and (self.min_val is not None or self.max_val is not None):
+            raise ValueError(f"{self.operand_type} operand cannot accept both equality and inequality targets")
+        if all(x is None for x in [self.target, self.min_val, self.max_val]):
+            self.target = self.value # No target has been defined, default it to the current value
 
     @property
     def value(self):
         """Get current value of the operand"""
-        metric_function = operand_registry.get(self.type)
+        metric_function = operand_registry.get(self.operand_type)
         if metric_function:
             return metric_function(**self.input_data)
         else:
-            raise ValueError(f'Unknown operand type: {self.type}')
+            raise ValueError(f'Unknown operand type: {self.operand_type}')
 
-    def delta(self):
-        """Calculate the difference between the target and current value"""
-        return (self.value - self.target)
+    def delta_target(self):
+        """Calculate the difference between the value and target"""
+        return self.value - self.target
+
+    def delta_ineq(self):
+        """Calculate the difference between the value and targets.
 
+        If the value is on the right side of the bound(s), 
+        then this operand simply is zero. 
+        Otherwise, it is the distance to the closest bound.
+        """
+        lower_penalty = max(0, self.min_val - self.value) if self.min_val is not None else 0
+        upper_penalty = max(0, self.value - self.max_val) if self.max_val is not None else 0
+        return lower_penalty + upper_penalty
+    
+    def delta(self):
+        """Calculate the difference to target"""
+        if self.target is not None:
+            return self.delta_target()
+        elif self.min_val is not None or self.max_val is not None:
+            return self.delta_ineq()
+        else:
+            raise ValueError(f"{self.operand_type} operand cannot compute delta")
+        
     def fun(self):
         """Calculate the objective function value"""
         return self.weight * self.delta()
 
     def __str__(self):
         """Return a string representation of the operand"""
-        return self.type.replace('_', ' ')
+        return self.operand_type.replace('_', ' ')
diff --git a/optiland/optimization/operand/operand_manager.py b/optiland/optimization/operand/operand_manager.py
@@ -17,16 +17,19 @@ class OperandManager:
     def __init__(self):
         self.operands = []
 
-    def add(self, operand_type, target, weight=1, input_data={}):
+    def add(self, operand_type=None, target=None, min_val=None, max_val=None, weight=1, input_data={}):
         """Add an operand to the merit function
 
         Args:
             operand_type (str): The type of the operand.
             target (float): The target value of the operand.
+            min_val (float): The operand should stay above this value (inequality operand).
+            max_val (float): The operand should stay below this value (inequality operand).
             weight (float): The weight of the operand.
             input_data (dict): Additional input data for the operand.
+        
         """
-        self.operands.append(Operand(operand_type, target, weight, input_data))
+        self.operands.append(Operand(operand_type, target, min_val, max_val, weight, input_data))
 
     def clear(self):
         """Clear all operands from the merit function"""
diff --git a/optiland/optimization/operand/ray.py b/optiland/optimization/operand/ray.py
@@ -163,6 +163,13 @@ def z_intercept_lcs(optic, surface_number, Hx, Hy, Px, Py, wavelength):
         optic.trace_generic(Hx, Hy, Px, Py, wavelength)
         intercept = optic.surface_group.z[surface_number, 0]
         decenter = optic.surface_group.surfaces[surface_number].geometry.cs.z
+        
+        # For some reason decenter can sometimes be a single-element array.
+        # In that case, retreive the float inside.
+        # This is a workaround until a solution is found.
+        if type(decenter) == np.ndarray:
+            decenter = decenter.item()
+
         return intercept - decenter
     
     @staticmethod
diff --git a/optiland/optimization/optimization.py b/optiland/optimization/optimization.py
@@ -44,9 +44,9 @@ def __init__(self):
         self.variables = VariableManager()
         self.initial_value = 0.0
 
-    def add_operand(self, operand_type, target, weight=1, input_data={}):
+    def add_operand(self, operand_type=None, target=None, min_val=None, max_val=None, weight=1, input_data={}):
         """Add an operand to the merit function"""
-        self.operands.add(operand_type, target, weight, input_data)
+        self.operands.add(operand_type, target, min_val, max_val, weight, input_data)
 
     def add_variable(self, optic, variable_type, **kwargs):
         """Add a variable to the merit function"""
@@ -84,16 +84,18 @@ def update_optics(self):
 
     def operand_info(self):
         """Print information about the operands in the merit function"""
-        data = {'Operand Type': [op.type.replace('_', ' ')
+        data = {'Operand Type': [op.operand_type.replace('_', ' ')
                                  for op in self.operands],
-                'Target': [op.target for op in self.operands],
+                'Target': [f'{op.target:+.3f}' if op.target is not None else ''  for op in self.operands],
+                'Min Bound': [op.min_val if op.min_val else '' for op in self.operands],
+                'Max Bound': [op.max_val if op.max_val else '' for op in self.operands],
                 'Weight': [op.weight for op in self.operands],
-                'Value': [op.value for op in self.operands],
-                'Delta': [op.delta() for op in self.operands]}
+                'Value': [f'{op.value:+.3f}' for op in self.operands],
+                'Delta': [f'{op.delta():+.3f}' for op in self.operands]}
 
         df = pd.DataFrame(data)
         funs = self.fun_array()
-        df['Contribution (%)'] = funs / np.sum(funs) * 100
+        df['Contrib. [%]'] = np.round(funs / np.sum(funs) * 100, 2)
 
         print(df.to_markdown(headers='keys', tablefmt='fancy_outline'))
 
@@ -119,7 +121,7 @@ def merit_info(self):
             improve_percent = ((self.initial_value - current_value) /
                                self.initial_value * 100)
 
-        data = {'Merit Function Value': [self.sum_squared()],
+        data = {'Merit Function Value': [current_value],
                 'Improvement (%)': improve_percent}
         df = pd.DataFrame(data)
         print(df.to_markdown(headers='keys', tablefmt='fancy_outline'))
@@ -211,11 +213,16 @@ def _fun(self, x):
         Returns:
             rss (float): The residual sum of squares.
         """
+
+        # Update all variables to their new values
         for idvar, var in enumerate(self.problem.variables):
             var.update(x[idvar])
-        self.problem.update_optics()  # update all optics (e.g., pickups)
-        funs = np.array([op.fun() for op in self.problem.operands])
-        rss = np.sum(funs**2)
+        
+        # Update optics (e.g., pickups and solves)
+        self.problem.update_optics()
+
+        # Compute merit function value
+        rss = self.problem.sum_squared()
         if np.isnan(rss):
             return 1e10
         else:
diff --git a/optiland/surfaces/surface_factory.py b/optiland/surfaces/surface_factory.py
@@ -155,7 +155,7 @@ def _configure_cs(self, index, **kwargs):
             elif index == 1:
                 z = 0  # first surface, always at zero
             else:
-                z = float(self._surface_group.positions[index-1]) + \
+                z = float(self._surface_group.positions[index-1].item()) + \
                     self.last_thickness
 
                 # self.last_thickness = thickness
diff --git a/optiland/tolerancing/core.py b/optiland/tolerancing/core.py
@@ -58,21 +58,20 @@ def __init__(self, optic, method='generic', tol=1e-5):
         self.compensator = CompensatorOptimizer(method=method, tol=tol)
 
     def add_operand(self, operand_type: str, input_data: dict = {},
-                    target: float = None, weight: float = 1.0):
+                    target: float = None, weight: float = 1.0,
+                    min_val: float = None,  max_val: float = None):
         """
         Add an operand to the tolerancing problem.
 
         Args:
-            operand_type: The type of the operand.
-            input_data: A dictionary of input data for the operand. Defaults to
-                an empty dictionary.
-            target: The target value for the operand. Defaults to None, in
-                which case the target is set based on the current value of
-                the operand.
-            weight: The weight of the operand for optimization during
-                compensation. Defaults to 1.0.
+            operand_type (str): The type of the operand.
+            target (float): The target value of the operand (equality operand).
+            min_val (float): The operand should stay above this value (inequality operand).
+            max_val (float): The operand should stay below this value (inequality operand).
+            weight (float): The weight of the operand.
+            input_data (dict): Additional input data for the operand.
         """
-        new_operand = Operand(operand_type, target, weight, input_data)
+        new_operand = Operand(operand_type, target, min_val, max_val, weight, input_data)
         if target is None:
             new_operand.target = new_operand.value
         self.operands.append(new_operand)
diff --git a/tests/test_compensator.py b/tests/test_compensator.py
@@ -32,7 +32,7 @@ def test_run_optimizer_generic():
     optic = Edmund_49_847()
     optimizer = CompensatorOptimizer(method='generic')
     optimizer.add_variable(optic, 'radius', surface_number=1)
-    optimizer.add_operand('f2', 25, 1, {'optic': optic})
+    optimizer.add_operand(operand_type='f2', target=25, weight=1, input_data={'optic': optic})
     result = optimizer.run()
     assert result is not None
 
@@ -41,7 +41,7 @@ def test_run_optimizer_least_squares():
     optic = Edmund_49_847()
     optimizer = CompensatorOptimizer(method='least_squares')
     optimizer.add_variable(optic, 'radius', surface_number=1)
-    optimizer.add_operand('f2', 25, 1, {'optic': optic})
+    optimizer.add_operand(operand_type='f2', target=2525, weight=1, input_data={'optic': optic})
     result = optimizer.run()
     assert result is not None
 
diff --git a/tests/test_operand.py b/tests/test_operand.py
@@ -158,6 +158,24 @@ def test_z_intercept(self, hubble):
         assert np.isclose(operand.RayOperand.z_intercept(**data),
                           6347.146837237045)
 
+    def test_x_intercept_lcs(self, hubble):
+        data = {'optic': hubble, 'surface_number': -1, 'Hx': 1.0, 'Hy': 0.0,
+                'Px': 0.0, 'Py': 0.0, 'wavelength': 0.55}
+        assert np.isclose(operand.RayOperand.x_intercept_lcs(**data),
+                          -150.42338010762842)
+
+    def test_y_intercept_lcs(self, hubble):
+        data = {'optic': hubble, 'surface_number': -1, 'Hx': 0.0, 'Hy': 1.0,
+                'Px': 0.0, 'Py': 0.0, 'wavelength': 0.55}
+        assert np.isclose(operand.RayOperand.y_intercept_lcs(**data),
+                          150.42338010762842)
+
+    def test_z_intercept_lcs(self, hubble):
+        data = {'optic': hubble, 'surface_number': -1, 'Hx': 0.0, 'Hy': 1.0,
+                'Px': 0.0, 'Py': 0.0, 'wavelength': 0.55}
+        assert np.isclose(operand.RayOperand.z_intercept_lcs(**data), 
+                          -18.062712762936826)  # Because Hubble's image is curved, otherwise it would be 0
+        
     def test_L(self, hubble):
         data = {'optic': hubble, 'surface_number': -1, 'Hx': 0.0, 'Hy': 1.0,
                 'Px': 0.0, 'Py': 0.0, 'wavelength': 0.55}
@@ -208,23 +226,23 @@ def test_opd_diff_new_dist(self, hubble):
 class TestOperand:
     def test_get_value(self, hubble):
         input_data = {'optic': hubble}
-        op = operand.Operand('f2', 1, 1, input_data)
+        op = operand.Operand(operand_type='f2', target=1, weight=1, input_data=input_data)
         assert np.isclose(op.value, 57600.080998403595)
 
     def test_invalid_operand(self, hubble):
         input_data = {'optic': hubble}
-        op = operand.Operand('f3', 1, 1, input_data)
+        op = operand.Operand(operand_type='f3', target=1, weight=1, input_data=input_data)
         with pytest.raises(ValueError):
             op.value
 
     def test_delta(self, hubble):
         input_data = {'optic': hubble}
-        op = operand.Operand('f2', 5000, 1, input_data)
+        op = operand.Operand(operand_type='f2', target=5000, weight=1, input_data=input_data)
         assert np.isclose(op.delta(), 52600.080998403595)
 
     def test_fun(self, hubble):
         input_data = {'optic': hubble}
-        op = operand.Operand('f2', 1e5, 1.5, input_data)
+        op = operand.Operand(operand_type='f2', target=1e5, weight=1.5, input_data=input_data)
         assert np.isclose(op.fun(), -63599.87850239461)
 
     def test_reregister_operand(self):
@@ -268,7 +286,7 @@ def test_getitem(self):
     def test_setitem(self):
         manager = operand.OperandManager()
         manager.add('f1', 1)
-        manager[0] = operand.Operand('f2', 1, 1, {})
+        manager[0] = operand.Operand(operand_type='f2', target=1, weight=1, input_data={})
         assert len(manager) == 1
 
     def test_len(self):
diff --git a/tests/test_optimization.py b/tests/test_optimization.py
diff --git a/tests/test_tolerancing.py b/tests/test_tolerancing.py
