feature: Improved Progress Bar to iteration based + showing convergence as % (#10)

feature: Improved Progress Bar: iteration based + showing convergence as %

Author: elenafuengar

.github/workflows/nightly_tests.yml

@@ -28,7 +28,7 @@ jobs:
     - name: pip install
       run:  |
         pip install wakis['notebook']
-        pip install iddefix
+        pip install iddefix pymoo
 
     - name: Print versions
       run: conda list

examples/005_cmaes_examples.ipynb

@@ -15,6 +15,9 @@ class SmartBoundDetermination:
 
     def __init__(self, frequency_data, impedance_data,
                  minimum_peak_height=1.0,
+                 threshold=None,
+                 distance=None,
+                 prominence=None,
                  Rs_bounds=[0.8, 10],
                  Q_bounds=[0.5, 5],
                  fres_bounds=[-0.01e9, +0.01e9]):
@@ -34,6 +37,16 @@ def __init__(self, frequency_data, impedance_data,
             Impedance magnitude data in Ohms.
         minimum_peak_height : float, optional
             Minimum height for a peak to be considered a resonance. Default is 1.0.
+        threshold : float, optional
+            Required vertical distance between a peak and its neighboring values
+            to be considered a peak. Passed to `scipy.signal.find_peaks`. Default is None.
+        distance : float, optional
+            Required minimum horizontal distance (in indices) between peaks.
+            Passed to `scipy.signal.find_peaks`. Default is None.
+        prominence : float, optional
+            Required prominence of peaks. The prominence measures how much a peak
+            stands out compared to its surrounding values. Passed to `scipy.signal.find_peaks`.
+            Default is None.
         Rs_bounds : list, optional
             Scaling factors [min, max] for Rs bounds. Default is [0.8, 10].
         Q_bounds : list, optional
@@ -78,11 +91,34 @@ def __init__(self, frequency_data, impedance_data,
         - Computed parameter bounds are stored in `self.parameterBounds`.
         - The `inspect()` method visualizes peak detection results.
         - The `to_table()` method prints a structured table of parameter ranges.
+
+        Returns
+        -------
+        parameterBounds : list of tuples
+            A list of parameter bounds for fitting. Each resonance contributes
+            three sets of bounds:
+            - `(Rs_min, Rs_max)`: Bounds for resistance Rs.
+            - `(Q_min, Q_max)`: Bounds for quality factor Q.
+            - `(freq_min, freq_max)`: Bounds for the resonant frequency.
+
+        Notes
+        -----
+        - The peak-finding algorithm is implemented using `scipy.signal.find_peaks`.
+        See the official documentation for more details:
+        https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html
+        - The 3dB bandwidth method is used to estimate initial Q factors and
+        define frequency bounds.
+        - The detected peaks and their heights are stored in instance attributes
+        `self.peaks` and `self.peaks_height`, respectively.
+        - The number of detected resonances is stored in `self.Nres`.
         """
 
         self.frequency_data = frequency_data
         self.impedance_data = impedance_data
         self.minimum_peak_height = minimum_peak_height
+        self.threshold = threshold
+        self.distance = distance
+        self.prominence = prominence
         self.Rs_bounds = Rs_bounds
         self.Q_bounds = Q_bounds
         self.fres_bounds = fres_bounds
@@ -159,7 +195,12 @@ def find(self, frequency_data=None, impedance_data=None,
             impedance_data = self.impedance_data
         if minimum_peak_height is None:
             minimum_peak_height = self.minimum_peak_height
-
+        if threshold is None:
+            threshold = self.threshold
+        if distance is None:
+            distance = self.distance
+        if prominence is None:
+            prominence = self.prominence
 
         # Find the peaks of the impedance data
         peaks, peaks_height = find_peaks(impedance_data,
@@ -197,6 +238,8 @@ def find(self, frequency_data=None, impedance_data=None,
             Q_bounds = (initial_Qs[i]*self.Q_bounds[0] , initial_Qs[i]*self.Q_bounds[1])
             freq_bounds = (frequency_data[peaks[i]]+self.fres_bounds[0], frequency_data[peaks[i]]+self.fres_bounds[1])
 
+            if peaks_height['peak_heights'][i] < 0:
+                Rs_bounds = (Rs_bounds[1], Rs_bounds[0])  # Swap for negative peaks
             parameterBounds.extend([Rs_bounds, Q_bounds, freq_bounds])
 
         # Store peaks and peaks_height as instance attributes
@@ -223,7 +266,7 @@ def inspect(self):
 
         return None
 
-    def to_table(self, to_markdown=False):
+    def to_table(self, parameterBounds=None, to_markdown=False):
         """
         Displays resonance parameters in a formatted ASCII table.
 
@@ -240,7 +283,7 @@ def to_table(self, to_markdown=False):
         2      |  85.61 to 864.12       |  120.55 to 200.23|  5.30e+08 to 7.23e+08
         ------------------------------------------------------------
         """
-        params = self.parameterBounds
+        params = self.parameterBounds if parameterBounds is None else parameterBounds
         N_resonators = len(params) // 3  # Compute number of resonators
 
         # Define formatting

iddefix/smartBoundDetermination.py

@@ -12,19 +12,45 @@
 from scipy.optimize import differential_evolution
 
 class ProgressBarCallback:
-    def __init__(self, max_generations):
+    def __init__(self, max_generations, desc="Optimization"):
         self.max_generations = max_generations
         self.current_generation = 0
-        self.previous_convergence = 0
-        self.pbar = tqdm(total=100, desc="Optimization Progress %")
-
-    def __call__(self, xk, convergence):
-        self.current_generation += 1
-        self.pbar.update((convergence-self.previous_convergence)*100)  # Update the progress bar
-        self.previous_convergence = convergence
-        if convergence > 1.:  # Convergence threshold
-            self.pbar.close()
-            return True  # Stop optimization early
+        self.pbar = tqdm(total=max_generations, desc=desc, unit="gen")
+
+    def __call__(self, *args):
+
+        # --- scipy differential_evolution: (xk, convergence) ---
+        if len(args) == 2 and not hasattr(args[0], "evaluator"):
+            xk, convergence = args
+            self.current_generation += 1
+            self.pbar.update(1)
+            self.pbar.set_postfix({
+                "conv": f"{(100*(convergence)):6.1f} %"
+            })
+            # optional early stopping
+            if convergence >= 1.0:
+                self.pbar.close()
+                return True
+        
+        # --- pymoo minimize: (algorithm) ---
+        elif len(args) == 1:
+            algorithm = args[0]
+            self.current_generation += 1
+            self.pbar.update(1)
+
+            # compute convergence metric
+            F = np.atleast_1d(algorithm.pop.get("F"))
+            gap = float(np.mean(F) - np.min(F))
+            self._gap0 = getattr(self, "_gap0", gap if gap > 0 else 1.0)
+            convergence = float(np.clip(1.0 - gap / (self._gap0 + 1e-12), 0.0, 1.0))
+
+            self.pbar.set_postfix({
+                "conv": f"{(100*(convergence)):6.1f} %"
+            })
+            # optional early stopping
+            if convergence >= 1.0 and self.current_generation > 1:
+                self.pbar.close()
+                return True
 
     def close(self):
         self.pbar.close()
@@ -73,7 +99,7 @@ def run_scipy_solver(parameterBounds,
                 - The solution found by the solver.
                 - A message indicating the solver's status.
         """
-        pbar = ProgressBarCallback(maxiter)
+        pbar = ProgressBarCallback(maxiter, desc='Differential Evolution')
         result = differential_evolution(minimization_function,
                                         parameterBounds,
                                         popsize=popsize,
@@ -91,22 +117,6 @@ def run_scipy_solver(parameterBounds,
                                     )
         pbar.close()
 
-        # Need to be reworked to use the last population as the new initial population to speed up convergence
-        """while ((result.message == 'Maximum number of iterations has been exceeded.') and (iteration_convergence)):
-            warning = 'Increased number of iterations by 10% to reach convergence. \n'
-            maxiter = int(1.1*maxiter)
-            result = differential_evolution(minimization_function,parameterBounds,
-                                            popsize=popsize, tol=tol, maxiter=maxiter,
-                                            mutation=mutation, recombination=crossover_rate, polish=False,
-                                            init='latinhypercube',
-                                            callback=show_progress_bar,
-                                            updating='deferred', workers=-1, #vectorized=vectorized
-                                        )
-
-        else:
-            warning = ''
-        """
-
         solution, message = result.x, result.message
 
         return solution, message
@@ -206,8 +216,9 @@ def run_pyfde_jade_solver(parameterBounds,
     def run_pymoo_cmaes_solver(parameterBounds,
                             minimization_function,
                             sigma=0.1,
-                            maxiter=1000,
-                            popsize=50,
+                            maxiter=None, # default: 100 + 150 * (N+3)**2 // popsize**0.5
+                            popsize=None, # defaul: 4 + int(3 * np.log(len(parameterBounds)))
+                            verbose=False,
                             **kwargs):
         """
         Runs the pymoo CMAES solver to minimize a given function.
@@ -229,10 +240,10 @@ def run_pymoo_cmaes_solver(parameterBounds,
             from pymoo.core.problem import Problem
             from pymoo.optimize import minimize
             from pymoo.termination import get_termination
-        except:
-            ImportError('''Please install the pymoo package to use the CMA-ES solver:
-                           >>> pip install pymoo
-                        ''')
+        except ImportError:
+            raise ImportError('''Please install the pymoo package to use the CMA-ES solver:
+                            >>> pip install pymoo
+                            ''')
 
         class OptimizationProblem(Problem):
             def __init__(self, objective_function, n_var, n_obj, xl, xu):
@@ -262,14 +273,25 @@ def _evaluate(self, x, out):
             sigma=sigma,
             popsize=popsize,
             seed=42,
+            restarts=3,
+            restart_from_best=True,
             **kwargs,
         )
         # use ftol and n_gen as stopping criteria
         termination_criteria = get_termination("n_gen", maxiter)
 
-        res = minimize(problem, solver, termination_criteria, seed=42, verbose=True)
+        if not verbose:
+            cb = ProgressBarCallback(maxiter, desc="CMA-ES evolution")  
+        else: cb = None
+
+        res = minimize(problem, solver, termination_criteria, 
+                       seed=42, 
+                       callback=cb,
+                       verbose=verbose)
+
+        if not verbose: cb.close()
 
         solution = res.X
-        message = "Convergence achieved" #if res. < maxiter else "Maximum iterations reached"
+        message = "Convergence achieved" if res.algorithm.n_gen < maxiter else "Maximum iterations reached"
 
         return solution, message, res