add hybrid Differential Evolution alghoritm

user_libs/DebyeFit/Debye_Fit.py

@@ -52,9 +52,10 @@ class Relaxation(object):
                         (not neglected as default: True).
     :type save: bool, optional, default:True
     :param optimizer: chosen optimization method:
-                                        Particle Swarm, Genetic or Dual Annealing
-                                        (Default: Particle_swarm).
-    :type optimizer: Optimizer class, optional, default:Particle_swarm
+                                        Hybrid Particle Swarm-Damped Least-Squares,
+                                        Genetic or Dual Annealing (DA)
+                                        (Default: PSO_DLS).
+    :type optimizer: Optimizer class, optional
     :param optimizer_options: Additional keyword arguments passed to
                                      optimizer class (Default: empty dict).
     :type optimizer_options: dict, optional, default: empty dict
@@ -64,7 +65,7 @@ class Relaxation(object):
                  material_name, f_n=50,
                  number_of_debye_poles=-1,
                  plot=True, save=True,
-                 optimizer=Particle_swarm,
+                 optimizer=PSO_DLS,
                  optimizer_options={}):
         self.name = 'Relaxation function'
         self.params = {}
@@ -139,13 +140,13 @@ class Relaxation(object):
         """ Calling the main optimisation module with defined lower and upper boundaries of search.
 
         Returns:
-            xpm (ndarray): The logarithm with base 10 of relaxation times of the Debyes poles.
-            mx (ndarray): Resulting optimised weights for the given relaxation times.
+            tau (ndarray): The optimised relaxation times.
+            weights (ndarray): Resulting optimised weights for the given relaxation times.
             ee (float): Average error between the actual and the approximated real part.
-            rp (ndarray): The real part of the permittivity for the optimised relaxation
-                          times and weights for the frequnecies included in freq.
-            ip (ndarray): The imaginary part of the permittivity for the optimised
-                          relaxation times and weights for the frequnecies included in freq.
+            rl (ndarray): Real parts of chosen relaxation function
+                          for given frequency points.
+            im (ndarray): Imaginary parts of chosen relaxation function
+                          for given frequency points.
         """
         # Define the lower and upper boundaries of search
         lb = np.full(self.number_of_debye_poles,
@@ -153,14 +154,13 @@ class Relaxation(object):
         ub = np.full(self.number_of_debye_poles,
                      -np.log10(np.min(self.freq)) + 3)
         # Call optimizer to minimize the cost function
-        xmp, _ = self.optimizer.fit(func=cost_function,
-                                    lb=lb, ub=ub,
-                                    funckwargs={'rl_g': self.rl,
-                                                'im_g': self.im,
-                                                'freq_g': self.freq}
-                                    )
-        _, _, mx, ee, rp, ip = DLS(self.rl, self.im, xmp, self.freq)
-        return xmp, mx, ee, rp, ip
+        tau, weights, ee, rl, im = self.optimizer.fit(func=self.optimizer.cost_function,
+                                                      lb=lb, ub=ub,
+                                                      funckwargs={'rl': self.rl,
+                                                                  'im': self.im,
+                                                                  'freq': self.freq}
+                                                     )
+        return tau, weights, ee, rl, im
 
     def run(self):
         """ Solve the problem described by the given relaxation function
@@ -182,6 +182,10 @@ class Relaxation(object):
         xmp, mx, ee, rp, ip = self.optimize()
         # Print the results in gprMax format style
         properties = self.print_output(xmp, mx, ee)
+        err_real, err_imag = self.error(rp + ee, ip)
+        print(f'The average fractional error for:\n'
+              f'- real part: {err_real}\n'
+              f'- imaginary part: {err_imag}\n')
         if self.save:
             self.save_result(properties)
         # Plot the actual and the approximate dielectric properties
@@ -265,6 +269,25 @@ class Relaxation(object):
         ax.set_ylabel("Approximation error (%)")
         plt.show()
 
+    def error(self, rl_exp, im_exp):
+        """ Calculate the average fractional error separately for
+        relative permittivity (real part) and conductivity (imaginary part)
+
+        Args:
+            rl_exp (ndarray): Real parts of optimised Debye expansion
+                              for given frequency points (plus average error).
+            im_exp (ndarray): Imaginary parts of optimised Debye expansion
+                              for given frequency points.
+        Returns:
+            avg_err_real (float): average fractional error
+                                  for relative permittivity (real part)
+            avg_err_imag (float): average fractional error
+                                  for conductivity (imaginary part)
+        """
+        avg_err_real = np.sum(np.abs((rl_exp - self.rl)/self.rl) * 100)/len(rl_exp)
+        avg_err_imag = np.sum(np.abs((im_exp - self.im)/self.im) * 100)/len(im_exp)
+        return avg_err_real, avg_err_imag
+
     @staticmethod
     def save_result(output, fdir="../materials"):
         """ Save the resulting Debye parameters in a gprMax format.
@@ -328,7 +351,7 @@ class HavriliakNegami(Relaxation):
                  sigma, mu, mu_sigma, material_name,
                  number_of_debye_poles=-1, f_n=50,
                  plot=False, save=True,
-                 optimizer=Particle_swarm,
+                 optimizer=PSO_DLS,
                  optimizer_options={}):
         super(HavriliakNegami, self).__init__(sigma=sigma, mu=mu, mu_sigma=mu_sigma,
                                               material_name=material_name, f_n=f_n,
@@ -398,7 +421,7 @@ class Jonscher(Relaxation):
                  sigma, mu, mu_sigma,
                  material_name, number_of_debye_poles=-1,
                  f_n=50, plot=False, save=True,
-                 optimizer=Particle_swarm,
+                 optimizer=PSO_DLS,
                  optimizer_options={}):
         super(Jonscher, self).__init__(sigma=sigma, mu=mu, mu_sigma=mu_sigma,
                                        material_name=material_name, f_n=f_n,
@@ -431,7 +454,7 @@ class Jonscher(Relaxation):
         if self.n_p > 1:
             sys.exit("n_p value must range between 0-1 (0 <= n_p <= 1)")
         if self.f_min == self.f_max:
-            sys.exit("Error: Null frequency range")
+            sys.exit("Error: Null frequency range!")
 
     def calculation(self):
         """Calculates the Q function for the given parameters"""
@@ -462,7 +485,7 @@ class Crim(Relaxation):
                  materials, sigma, mu, mu_sigma, material_name, 
                  number_of_debye_poles=-1, f_n=50,
                  plot=False, save=True,
-                 optimizer=Particle_swarm,
+                 optimizer=PSO_DLS,
                  optimizer_options={}):
 
         super(Crim, self).__init__(sigma=sigma, mu=mu, mu_sigma=mu_sigma,
@@ -553,7 +576,7 @@ class Rawdata(Relaxation):
                  material_name, number_of_debye_poles=-1,
                  f_n=50, delimiter =',',
                  plot=False, save=True,
-                 optimizer=Particle_swarm,
+                 optimizer=PSO_DLS,
                  optimizer_options={}):
 
         super(Rawdata, self).__init__(sigma=sigma, mu=mu, mu_sigma=mu_sigma,
@@ -619,10 +642,20 @@ if __name__ == "__main__":
                             material_name="Kelley",
                             number_of_debye_poles=5, f_n=100,
                             plot=True, save=False,
-                            optimizer=Dual_annealing,
+                            optimizer=DA,
                             optimizer_options={'seed': 111})
     setup.run()
-    ### Testing setup
+    setup = HavriliakNegami(f_min=1e7, f_max=1e11,
+                            alpha=1-0.09, beta=0.45,
+                            e_inf=2.7, de=8.6-2.7, tau_0=9.4e-10,
+                            sigma=0, mu=0, mu_sigma=0,
+                            material_name="Kelley",
+                            number_of_debye_poles=5, f_n=100,
+                            plot=True, save=False,
+                            optimizer=DE,
+                            optimizer_options={'seed': 111})
+    setup.run()
+    '''### Testing setup
     setup = Rawdata("Test.txt", 0.1, 1, 0.1, "M1", 3, plot=True,
                     optimizer_options={'seed': 111,
                                        'pflag': True})
@@ -640,4 +673,4 @@ if __name__ == "__main__":
     materials = np.array([material1, material2])
     setup = Crim(1*1e-1, 1e-9, 0.5, f, materials, 0.1,
                  1, 0, "M4", 2, plot=True)
-    setup.run()
+    setup.run()'''

user_libs/DebyeFit/optimization.py

@@ -24,7 +24,7 @@ from tqdm import tqdm
 
 class Optimizer(object):
     """
-    Create particle swarm optimisation object.
+    Create choosen optimisation object.
 
     :param maxiter: The maximum number of iterations for the
                     optimizer (Default: 1000).
@@ -36,14 +36,57 @@ class Optimizer(object):
     def __init__(self, maxiter=1000, seed=None):
         self.maxiter = maxiter
         self.seed = seed
+        self.calc_weights = DLS
 
-    def fit(self):
-        """
-        Call the optimization function that tries to find an optimal set
+    def fit(self, func, lb, ub, funckwargs={}):
+        """Call the optimization function that tries to find an optimal set
         of relaxation times that minimise the error
-        between the actual and the approximated electric permittivity.
+        between the actual and the approximated electric permittivity
+        and calculate optimised weights for the given relaxation times.
+
+        Args:
+            func (function): The function to be minimized.
+            lb (ndarray): The lower bounds of the design variable(s).
+            ub (ndarray): The upper bounds of the design variable(s).
+            funckwargs (dict): Additional keyword arguments passed to
+                               objective and constraint function
+                               (Default: empty dict).
+
+        Returns:
+            tau (ndarray): The the best relaxation times.
+            weights (ndarray): Resulting optimised weights for the given relaxation times.
+            ee (float): Average error between the actual and the approximated real part.
+            rl (ndarray): Real parts of chosen relaxation function
+                          for given frequency points.
+            im (ndarray): Imaginary parts of chosen relaxation function
+                          for given frequency points.
         """
-        raise NotImplementedError()
+        np.random.seed(self.seed)
+        # find the relaxation frequencies using Differential Evolution
+        tau, _ = self.calc_relaxation_times(func, lb, ub, funckwargs)
+        # find the weights using a damped least squares
+        _, _, weights, ee, rl_exp, im_exp = self.calc_weights(tau, **funckwargs)
+        return tau, weights, ee, rl_exp, im_exp
+
+    @staticmethod
+    def cost_function(x, rl, im, freq):
+        """
+        The cost function is the average error between
+        the actual and the approximated electric permittivity.
+
+        Args:
+            x (ndarray): The logarithm with base 10 of relaxation times of the Debyes poles.
+            rl (ndarray): Real parts of chosen relaxation function
+                          for given frequency points.
+            im (ndarray): Imaginary parts of chosen relaxation function
+                          for given frequency points.
+            freq (ndarray): The frequencies vector for defined grid.
+
+        Returns:
+            cost (float): Sum of mean absolute errors for real and imaginary parts.
+        """
+        cost_i, cost_r, _, _, _, _ = DLS(x, rl, im, freq)
+        return cost_i + cost_r
 
     @staticmethod
     def plot(x, y):
@@ -66,9 +109,9 @@ class Optimizer(object):
         plt.pause(0.0001)
 
 
-class Particle_swarm(Optimizer):
-    """
-    Create particle swarm optimisation object with predefined parameters.
+class PSO_DLS(Optimizer):
+    """ Create hybrid Particle Swarm-Damped Least Squares optimisation
+    object with predefined parameters.
 
     :param swarmsize: The number of particles in the swarm (Default: 40).
     :type swarmsize: int, optional
@@ -98,7 +141,7 @@ class Particle_swarm(Optimizer):
                  minstep=1e-8, minfun=1e-8,
                  pflag=False, seed=None):
 
-        super(Particle_swarm, self).__init__(maxiter, seed)
+        super(PSO_DLS, self).__init__(maxiter, seed)
         self.swarmsize = swarmsize
         self.omega = omega
         self.phip = phip
@@ -107,7 +150,7 @@ class Particle_swarm(Optimizer):
         self.minfun = minfun
         self.pflag = pflag
 
-    def fit(self, func, lb, ub, funckwargs={}):
+    def calc_relaxation_times(self, func, lb, ub, funckwargs={}):
         """
         A particle swarm optimisation that tries to find an optimal set
         of relaxation times that minimise the error
@@ -124,7 +167,7 @@ class Particle_swarm(Optimizer):
                                (Default: empty dict).
 
         Returns:
-            g (ndarray): The swarm's best known position (optimal design).
+            g (ndarray): The swarm's best known position (relaxation times).
             fg (float): The objective value at ``g``.
         """
         np.random.seed(self.seed)
@@ -217,75 +260,23 @@ class Particle_swarm(Optimizer):
                 else:
                     xpp.append(it)
                     ypp.append(fg)
-                Particle_swarm.plot(xpp, ypp)
+                PSO_DLS.plot(xpp, ypp)
         return g, fg
 
 
-class Dual_annealing(Optimizer):
-    """
-    Create dual annealing object with predefined parameters.
-
-    :param maxiter: The maximum number of iterations for the swarm
-                    to search (Default: 1000).
-    :type maxiter: int, optional
-    :param local_search_options: Extra keyword arguments to be passed
-                                 to the local minimizer, reffer to
-                                 scipy.optimize.minimize() function
-                                 (Default: empty dict).
-    :type local_search_options: dict, optional
-    :param initial_temp (float): The initial temperature, use higher values to
-                                facilitates a wider search of the energy
-                                landscape, allowing dual_annealing to escape
-                                local minima that it is trapped in.
-                                Range is (0.01, 5.e4] (Default: 5230).
-    :type initial_temp: float, optional
-    :param restart_temp_ratio: During the annealing process,
-                               temperature is decreasing, when it
-                               reaches initial_temp * restart_temp_ratio,
-                               the reannealing process is triggered.
-                               Range is (0, 1) (Default: 2e-5).
-    :type restart_temp_ratio: float, optional
-    :param visit: Parameter for visiting distribution. The value range is (1, 3]
-                  (Default: 2.62).
-    :type visit: float, optional
-    :param accept: Parameter for acceptance distribution. It is used to control
-                   the probability of acceptance. The lower the acceptance parameter,
-                   the smaller the probability of acceptance. The value range (-1e4, -5]
-                   (Default: -5.0).
-    :type accept: float, optional
-    :param no_local_search (bool): If no_local_search is set to True, a traditional
-                                   Generalized Simulated Annealing will be performed
-                                   with no local search strategy applied (Default: False).
-    :type no_local_search: bool, optional
-    :param maxfun: Soft limit for the number of objective function calls.
-                        (Default: 1e7).
-    :type maxfun: int, optional
-    :param callback: A callback function with signature callback(x, f, context),
-                     which will be called for all minima found.
-                     x and f are the coordinates and function value of
-                     the latest minimum found, and context has value in [0, 1, 2],
-                     with the following meaning:
-                     0: minimum detected in the annealing process.
-                     1: detection occurred in the local search process.
-                     2: detection done in the dual annealing process.
-                     If the callback implementation returns True,
-                     the algorithm will stop.
-    :type callback: None, callable, optional
-    :param x0: Coordinates of a single N-D starting point, shape(n,).
-                        (Default: None).
-    :type x0: None, ndarray, optional
-    :param seed: Specify seed for repeatable minimizations.
-                 The random numbers generated with this seed only
-                 affect the visiting distribution function and
-                 new coordinates generation (Default: None).
-    :type seed: None, int, optional
+class DA(Optimizer):
+    """ Create Dual Annealing object with predefined parameters.
+    The current class is a modified edition of the scipy.optimize
+    package which can be found at:
+    https://docs.scipy.org/doc/scipy/reference/generated/
+    scipy.optimize.dual_annealing.html#scipy.optimize.dual_annealing
     """
     def __init__(self, maxiter=1000,
                  local_search_options={}, initial_temp=5230.0,
                  restart_temp_ratio=2e-05, visit=2.62, accept=-5.0,
                  maxfun=1e7, no_local_search=False,
                  callback=None, x0=None, seed=None):
-        super(Dual_annealing, self).__init__(maxiter, seed)
+        super(DA, self).__init__(maxiter, seed)
         self.local_search_options = local_search_options
         self.initial_temp = initial_temp
         self.restart_temp_ratio = restart_temp_ratio
@@ -296,7 +287,7 @@ class Dual_annealing(Optimizer):
         self.callback = callback
         self.x0 = x0
 
-    def fit(self, func, lb, ub, funckwargs={}):
+    def calc_relaxation_times(self, func, lb, ub, funckwargs={}):
         """
         Find the global minimum of a function using Dual Annealing.
         The current class is a modified edition of the scipy.optimize
@@ -306,15 +297,15 @@ class Dual_annealing(Optimizer):
 
         Args:
             func (function): The function to be minimized
-            lb (array): The lower bounds of the design variable(s)
-            ub (array): The upper bounds of the design variable(s)
+            lb (ndarray): The lower bounds of the design variable(s)
+            ub (ndarray): The upper bounds of the design variable(s)
             funckwargs (dict): Additional keyword arguments passed to
                                objective and constraint function
                                (Default: empty dict)
 
         Returns:
-            g (array): The solution array (optimal design).
-            fg (float): The objective value at the solution.
+            x (ndarray): The solution array (relaxation times).
+            fun (float): The objective value at the best solution.
         """
         np.random.seed(self.seed)
         result = scipy.optimize.dual_annealing(func,
@@ -330,22 +321,92 @@ class Dual_annealing(Optimizer):
                                                no_local_search=self.no_local_search,
                                                callback=self.callback,
                                                x0=self.x0)
+        print(result.message)
         return result.x, result.fun
 
 
-def DLS(rl, im, logt, freq):
+class DE(Optimizer):
+    """
+    Create Differential Evolution-Damped Least Squares object with predefined parameters.
+    The current class is a modified edition of the scipy.optimize
+    package which can be found at:
+    https://docs.scipy.org/doc/scipy/reference/generated/
+    scipy.optimize.differential_evolution.html#scipy.optimize.differential_evolution
+    """
+    def __init__(self, maxiter=1000,
+                 strategy='best1bin', popsize=15, tol=0.01, mutation=(0.5, 1),
+                 recombination=0.7, callback=None, disp=False, polish=True,
+                 init='latinhypercube', atol=0, updating='immediate', workers=1,
+                 constraints=(), seed=None):
+        super(DE, self).__init__(maxiter, seed)
+        self.strategy = strategy
+        self.popsize = popsize
+        self.tol = tol
+        self.mutation = mutation
+        self.recombination = recombination
+        self.callback = callback
+        self.disp = disp
+        self.polish = polish
+        self.init= init
+        self.atol = atol
+        self.updating = updating
+        self.workers = workers
+        self.constraints = constraints
+
+    def calc_relaxation_times(self, func, lb, ub, funckwargs={}):
+        """
+        Find the global minimum of a function using Differential Evolution.
+        The current class is a modified edition of the scipy.optimize
+        package which can be found at:
+        https://docs.scipy.org/doc/scipy/reference/generated/
+        scipy.optimize.differential_evolution.html#scipy.optimize.differential_evolution
+
+        Args:
+            func (function): The function to be minimized
+            lb (ndarray): The lower bounds of the design variable(s)
+            ub (ndarray): The upper bounds of the design variable(s)
+            funckwargs (dict): Additional keyword arguments passed to
+                               objective and constraint function
+                               (Default: empty dict)
+
+        Returns:
+            x (ndarray): The solution array (relaxation times).
+            fun (float): The objective value at the best solution.
+        """
+        np.random.seed(self.seed)
+        result = scipy.optimize.differential_evolution(func,
+                                                       bounds=list(zip(lb, ub)),
+                                                       args=funckwargs.values(),
+                                                       strategy=self.strategy,
+                                                       popsize=self.popsize,
+                                                       tol=self.tol,
+                                                       mutation=self.mutation,
+                                                       recombination=self.recombination,
+                                                       callback=self.callback,
+                                                       disp=self.disp,
+                                                       polish=self.polish,
+                                                       init=self.init,
+                                                       atol=self.atol,
+                                                       updating=self.updating,
+                                                       workers=self.workers,
+                                                       constraints=self.constraints)
+        print(result.message)
+        return result.x, result.fun
+
+
+def DLS(logt, rl, im, freq):
     """
     Find the weights using a non-linear least squares (LS) method,
     the Levenberg–Marquardt algorithm (LMA or just LM),
     also known as the damped least-squares (DLS) method.
 
     Args:
+        logt (ndarray): The best known position form optimization module (optimal design),
+                        the logarithm with base 10 of relaxation times of the Debyes poles.
         rl (ndarray): Real parts of chosen relaxation function
                       for given frequency points.
         im (ndarray): Imaginary parts of chosen relaxation function
                       for given frequency points.
-        logt (ndarray): The best known position form optimization module (optimal design),
-                        the logarithm with base 10 of relaxation times of the Debyes poles.
         freq (ndarray): The frequencies vector for defined grid.
 
     Returns:
@@ -364,36 +425,16 @@ def DLS(rl, im, logt, freq):
     # logt=log10(t0) for efficiency during the optimisation
     # Here they are transformed back t0=10**logt
     tt = 10**logt
-    # y = Ax , here the A matrix for the real and the imaginary part is builded
+    # y = Ax, here the A matrix for the real and the imaginary part is builded
     d = 1 / (1 + 1j * 2 * np.pi * np.repeat(
              freq, len(tt)).reshape((-1, len(tt))) * tt)
     # Adding dumping (Levenberg–Marquardt algorithm)
     # Solving the overdetermined system y=Ax
-    x = np.abs(np.linalg.lstsq(d.imag, im, rcond=-1)[0])  # absolute damped least-squares solution
+    x = np.abs(np.linalg.lstsq(d.imag, im, rcond=None)[0])  # absolute damped least-squares solution
     rp, ip = np.matmul(d.real, x[np.newaxis].T).T[0], np.matmul(d.imag, x[np.newaxis].T).T[0]
     cost_i = np.sum(np.abs(ip-im))/len(im)
     ee = np.mean(rl - rp)
     if ee < 1:
         ee = 1
-    cost_r = np.sum(np.abs(rp - rl + ee))/len(im)
+    cost_r = np.sum(np.abs(rp + ee - rl))/len(im)
     return cost_i, cost_r, x, ee, rp, ip
-
-
-def cost_function(x, rl_g, im_g, freq_g):
-    """
-    The cost function is the average error between
-    the actual and the approximated electric permittivity.
-
-    Args:
-        x (ndarray): The logarithm with base 10 of relaxation times of the Debyes poles.
-        rl_g (ndarray): Real parts of chosen relaxation function
-                        for given frequency points.
-        im_g (ndarray): Imaginary parts of chosen relaxation function
-                        for given frequency points.
-        freq (ndarray): The frequencies vector for defined grid.
-
-    Returns:
-        cost (float): Sum of mean absolute errors for real and imaginary part.
-    """
-    cost_i, cost_r, _, _, _, _ = DLS(rl_g, im_g, x, freq_g)
-    return cost_i + cost_r