Feat/optimization algorithms (#62)

* fix: improve formatting in optimization.py

* feat: add SHGO optimizer to optimization module

* feat: implement SHGO optimizer in optimization module

* feat: add BasinHopping optimizer to optimization module

* feat: add BasinHopping optimizer to optimization module

* fix: optimize update method to conditionally call update_paraxial

* docs: update tutorial to include notes on SHGO and BasinHopping optimizers

* docs: update tutorial to include SHGO and basin-hopping in global optimization examples

* docs: improve formatting of docstring in EvenAsphere class

* docs: add SHGO optimization example to gallery documentation

* docs: update optimization gallery to include basin-hopping and improve descriptions

* docs: update optimization framework documentation to include global and local optimizers

Author: HarrisonKramer

docs/developers_guide/optimization_framework.rst

@@ -30,9 +30,12 @@ Components Explained
    - A base `Optimizer` class wraps `scipy.optimize.minimize` and provides a unified interface.
    - Built-in optimizers include:
 
-     - **Dual Annealing**
-     - **Differential Evolution**
-     - **Least Squares** (and more)
+     - **Dual Annealing** (global)
+     - **Differential Evolution** (global)
+     - **Basin Hopping** (global)
+     - **SHGO** (global)
+     - **Least Squares** (local)
+     - **Nelder-Mead**, **Powell**, **BFGS**, **L-BFGS-B**, **COBYLA**, etc. (local optimization, from `scipy.optimize.minimize`)
    - Users can subclass the base optimizer for custom methods.
 
 3. **Operands and Variables**:

docs/examples/Tutorial_5b_Advanced_Optimization.ipynb

@@ -26,7 +26,7 @@
     "\n",
     "- Various operand types (paraxial, real ray-based, aberrations)\n",
     "- Various variable types (radii, thickness, conic)\n",
-    "- Global optimization (e.g. dual annealing or differential evolution)"
+    "- Global optimization (e.g. dual annealing, differential evolution, SHGO, basin-hopping)"
    ]
   },
   {
@@ -224,6 +224,9 @@
    "outputs": [],
    "source": [
     "optimizer = optimization.DifferentialEvolution(problem)\n",
+    "# optimizer = optimization.SHGO(problem)  # note SHGO requires bounds\n",
+    "# optimizer = optimization.BasinHopping(problem)  # note BasinHopping requires no bounds\n",
+    "\n",
     "res = optimizer.optimize(maxiter=256, disp=False, workers=-1)"
    ]
   },

docs/gallery/optimization.rst

@@ -13,6 +13,8 @@ These examples demonstrate the optimization capabilities of optical systems usin
     optimization/constrained
     optimization/bounded_operands
     optimization/global
+    optimization/basin_hopping
+    optimization/shgo
     optimization/asphere
     optimization/beam_expander
     optimization/freeform

docs/gallery/optimization/basin_hopping.ipynb

@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Global Optimization"
+    "# Global Optimization (Differential Evolution)"
    ]
   },
   {
@@ -14,6 +14,8 @@
     "The following global optimizers are implemented in Optiland:\n",
     "1. Differential Evolution\n",
     "2. Dual Annealing\n",
+    "3. SHGO\n",
+    "4. Basin-hopping\n",
     "\n",
     "These optimizers wrap the scipy.optimize implementations."
    ]

docs/gallery/optimization/global.ipynb

@@ -11,9 +11,9 @@
 - k is the conic constant
 - Ci are the aspheric coefficients
 
-Even-order aspheric surfaces are commonly used to correct specific aberrations 
-while maintaining rotational symmetry. These surfaces are defined by polynomial 
-terms with even exponents, ensuring symmetry about the optical axis. 
+Even-order aspheric surfaces are commonly used to correct specific aberrations
+while maintaining rotational symmetry. These surfaces are defined by polynomial
+terms with even exponents, ensuring symmetry about the optical axis.
 
 Kramer Harrison, 2024
 """
@@ -34,9 +34,10 @@ class EvenAsphere(NewtonRaphsonGeometry):
     - k is the conic constant
     - Ci are the aspheric coefficients
 
-    Even-order aspheric surfaces are commonly used to correct specific aberrations 
-    while maintaining rotational symmetry. These surfaces are defined by polynomial 
-    terms with even exponents, ensuring symmetry about the optical axis. 
+    Even-order aspheric surfaces are commonly used to correct specific
+    aberrations while maintaining rotational symmetry. These surfaces are
+    defined by polynomial terms with even exponents, ensuring symmetry about
+    the optical axis.
 
     Args:
         coordinate_system (str): The coordinate system used for the geometry.

docs/gallery/optimization/shgo.ipynb

@@ -417,13 +417,16 @@ def update_normalization(self, surface) -> None:
         if surface.surface_type == 'zernike':
             surface.geometry.norm_radius = surface.semi_aperture*1.1
 
-    def update(self)->None:
+    def update(self) -> None:
         """
         Update the surfaces based on the pickup operations.
         """
         self.pickups.apply()
         self.solves.apply()
-        self.update_paraxial()
+
+        if any(surface.surface_type in ['chebyshev', 'zernike']
+           for surface in self.surface_group.surfaces):
+            self.update_paraxial()
 
     def image_solve(self):
         """Update the image position such that the marginal ray crosses the

optiland/geometries/even_asphere.py

@@ -22,5 +22,7 @@
     OptimizerGeneric,
     LeastSquares,
     DualAnnealing,
-    DifferentialEvolution
+    DifferentialEvolution,
+    SHGO,
+    BasinHopping
 )

optiland/optic.py

@@ -212,7 +212,7 @@ def undo(self):
                 var.update(x0[idvar])
             self._x.pop(-1)
 
-    def _fun(self, x)->float:
+    def _fun(self, x) -> float:
         """
         Internal function to evaluate the objective function.
 
@@ -231,7 +231,7 @@ def _fun(self, x)->float:
         self.problem.update_optics()
 
         # Compute merit function value
-        try:   
+        try:
             rss = self.problem.sum_squared()
             if np.isnan(rss):
                 return 1e10
@@ -405,3 +405,108 @@ def optimize(self, maxiter=1000, disp=True, workers=-1):
                                                      updating=updating,
                                                      workers=workers)
         return result
+
+
+class SHGO(OptimizerGeneric):
+    """
+    Simplicity Homology Global Optimization (SHGO).
+
+    Args:
+        problem (OptimizationProblem): The optimization problem to be solved.
+
+    Methods:
+        optimize(workers=-1, *args, **kwargs): Runs the SHGO algorithm.
+    """
+
+    def __init__(self, problem: OptimizationProblem):
+        """
+        Initializes a new instance of the SHGO class.
+
+        Args:
+            problem (OptimizationProblem): The optimization problem to be
+                solved.
+        """
+        super().__init__(problem)
+
+    def optimize(self, workers=-1, *args, **kwargs):
+        """
+        Runs the SHGO algorithm. Note that the SHGO algorithm accepts the same
+        arguments as the scipy.optimize.shgo function.
+
+        Args:
+            workers (int): Number of parallel workers to use. Set to -1 to use
+                all available CPU processors. Default is -1.
+            *args: Variable length argument list.
+            **kwargs: Arbitrary keyword arguments.
+
+        Returns:
+            result (OptimizeResult): The optimization result.
+
+        Raises:
+            ValueError: If any variable in the problem does not have bounds.
+        """
+        x0 = [var.value for var in self.problem.variables]
+        self._x.append(x0)
+        bounds = tuple([var.bounds for var in self.problem.variables])
+        if any(None in bound for bound in bounds):
+            raise ValueError('SHGO requires all variables have bounds.')
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=RuntimeWarning)
+
+            result = optimize.shgo(self._fun, bounds=bounds,
+                                   workers=workers, *args, **kwargs)
+        return result
+
+
+class BasinHopping(OptimizerGeneric):
+    """
+    Basin-hopping optimizer for solving optimization problems.
+
+    Args:
+        problem (OptimizationProblem): The optimization problem to be solved.
+
+    Methods:
+        optimize(maxiter=1000, disp=True, workers=-1): Runs the basin-hopping
+            optimization algorithm.
+    """
+
+    def __init__(self, problem: OptimizationProblem):
+        """
+        Initializes a new instance of the BasinHopping class.
+
+        Args:
+            problem (OptimizationProblem): The optimization problem to be
+                solved.
+        """
+        super().__init__(problem)
+
+    def optimize(self, niter=100, *args, **kwargs):
+        """
+        Runs the basin-hopping algorithm. Note that the basin-hopping
+        algorithm accepts the same arguments as the
+        scipy.optimize.basinhopping function.
+
+        Args:
+            niter (int): Number of iterations to perform. Default is 100.
+            *args: Variable length argument list.
+            **kwargs: Arbitrary keyword arguments.
+
+        Returns:
+            result (OptimizeResult): The optimization result.
+
+        Raises:
+            ValueError: If any variable in the problem does not have bounds.
+        """
+        x0 = [var.value for var in self.problem.variables]
+        self._x.append(x0)
+        bounds = tuple([var.bounds for var in self.problem.variables])
+        if not all(x is None for pair in bounds for x in pair):
+            raise ValueError('Basin-hopping does not accept bounds.')
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=RuntimeWarning)
+
+            result = optimize.basinhopping(self._fun, x0=x0,
+                                           niter=niter, *args, **kwargs)
+        return result