add more broad documentation

user_libs/DebyeFit/Readme.rst

@@ -8,11 +8,16 @@ In the ``user_libs`` sub-package is a module called ``DebyeFit`` which can be us
    \epsilon(\omega) = \epsilon_{\infty} + \sum_{i=1}^{N}\frac{\Delta\epsilon_{i}}{1+j\omega t_{0,i}},
 
 where :math:`\epsilon(\omega)` is frequency dependent dielectric properties, :math:`\Delta\epsilon` - difference between the real permittivity at zero and infinity frequency.
-:math:`\tau_{0}` is relaxation time,  :math:`\epsilon_{\infty}` - real part of relative permittivity at infinity frequency, and :math:`N` is number of the Debye poles.
+:math:`\tau_{0}` is relaxation time (s),  :math:`\epsilon_{\infty}` - real part of relative permittivity at infinity frequency, and :math:`N` is number of the Debye poles.
 
 The user can choose between Havriliak-Negami, Jonsher, Complex Refractive Index Mixing models, and arbitrary dielectric data derived experimentally
 or calculated using some other function.
 
+License
+=======
+
+``DebyeFit`` sub-package is currently released under the `GNU General Public License v3 or higher <http://www.gnu.org/copyleft/gpl.html>`_.
+
 Code structure
 ==============
 
@@ -21,12 +26,15 @@ The ``user_libs`` sub-package contains two main scripts:
 * ```Debye_fit.py``` with definition of all Relaxation functions classes,
 * ```optimization.py``` with definition of three choosen global optimization methods.
 
+
 Class Relaxation
 ################
 
 This class is designed for modelling different relaxation functions, like Havriliak-Negami (```Class HavriliakNegami```), Jonsher (```Class Jonsher```), Complex Refractive Index Mixing (```Class CRIM```) models, and arbitrary dielectric data derived experimentally
 or calculated using some other function (```Class Rawdata```).
 
+More about ``Class Relaxation`` structure can be found in [relaxation.md](./docs/relaxation.md).
+
 Havriliak-Negami Function
 *************************
 
@@ -56,10 +64,10 @@ The ``HavriliakNegami`` class has the following structure:
 * ``beta`` is real positive float number which varies 0 < $\beta$ < 1,
 * ``e_inf`` is a real part of relative permittivity at infinity frequency,
 * ``de`` is a difference between the real permittivity at zero and infinity frequency,
-* ``tau_0`` is a relaxation time,
+* ``tau_0`` is a relaxation time (s),
 * ``sigma`` is a conductivity (Siemens/metre),
 * ``mu`` is a relative permeability,
-* ``mu_sigma`` is a magnetic loss,
+* ``mu_sigma`` is a magnetic loss (Ohms/metre),
 * ``material_name`` is definition of material name,
 * ``number_of_debye_poles`` is choosen number of Debye poles,
 * ``f_n`` is choosen number of frequences,
@@ -157,32 +165,35 @@ Class Optimizer
 This class supports global optimization algorithms (particle swarm, duall annealing, evolutionary algorithms) for finding an optimal set of relaxation times that minimise the error between the actual and the approximated electric permittivity and calculate optimised weights for the given relaxation times.
 Code written here is mainly based on external libraries, like ```scipy``` and ```pyswarm```.
 
+More about ``Class Optimizer`` structure can be found in [optimization.md](docs/optimization.md).
+
 PSO_DLS Class
 *************
 
-Create hybrid Particle Swarm-Damped Least Squares optimisation object with predefined parameters.
+Creation of hybrid Particle Swarm-Damped Least Squares optimisation object with predefined parameters.
+The current code is a modified edition of the pyswarm package which can be found at https://pythonhosted.org/pyswarm/.
 
-DA Class
+DA_DLS Class
-********
+************
 
-Create Dual Annealing-Damped Least Squares optimisation object with predefined parameters. The current class is a modified edition of the scipy.optimize package which can be found at:
+Creation of Dual Annealing-Damped Least Squares optimisation 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
+https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.dual_annealing.html#scipy.optimize.dual_annealing.
 
-DE Class
+DE_DLS Class
-********
+************
 
-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:
+Creation of 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
+https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.differential_evolution.html#scipy.optimize.differential_evolution.
 
 DLS function
 ************
 
-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.
+Finding 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.
 
 Examples
 ########
 
-In directory ```examples```, we provided jupyter notebooks, scripts and data to show how use stand alone script ```DebyeFit.py```:
+In directory [examples](./examles), we provided jupyter notebooks, scripts and data to show how use stand alone script ```DebyeFit.py```:
 
 * ```example_DebyeFitting.ipynb```: simple cases of using all available implemented relaxation functions,
 * ```example_BiologicalTissues.ipynb```: simple cases of using Cole-Cole function for biological tissues,
@@ -219,13 +230,13 @@ Allows you to model dielectric properties by fitting multi-Debye functions to Ha
 
     #havriliak_negami: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 i1 str1 [i2]
 
-* ``f1`` is the lower frequency bound (in Hz).
+* ``f1`` is the lower frequency bound (Hz).
-* ``f2`` is the upper frequency bound (in Hz).
+* ``f2`` is the upper frequency bound (Hz).
 * ``f3`` is the :math:`\alpha` parameter beetwen bonds :math:`\left(0 < \alpha < 1 \right)`.
 * ``f4`` is the :math:`\beta` parameter beetwen bonds :math:`\left(0 < \beta < 1 \right)`.
 * ``f5`` is the real relative permittivity at infinity frequency, :math:`\epsilon_{\infty}`.
 * ``f6`` is the difference between the real permittivity at zero and infinity frequency, :math:`\Delta\epsilon`.
-* ``f7`` is the relaxation time, :math:`t_{0}`.
+* ``f7`` is the relaxation time (s), :math:`t_{0}`.
 * ``f8`` is the conductivity (Siemens/metre), :math:`\sigma`
 * ``f9`` is the relative permeability, :math:`\mu_r`
 * ``f10`` is the magnetic loss (Ohms/metre), :math:`\sigma_*`

user_libs/DebyeFit/docs/optimization.md

@@ -0,0 +1,48 @@
+# Optimization methods of multi-Debye fitting
+
+``Class Optimizer`` supports global optimization algorithms (particle swarm, duall annealing, evolutionary algorithms) for finding an optimal set of relaxation times that minimise the error between the actual and the approximated electric permittivity and calculate optimised weights for the given relaxation times.
+Code written here is mainly based on external libraries, like ```scipy``` and ```pyswarm```.
+
+Supported methods:
+- [x] hybrid Particle Swarm-Damped Least Squares
+- [x] hybrid Dual Annealing-Damped Least Squares
+- [x] hybrid Differential Evolution-Damped Least Squares
+
+## Methods
+
+1. __constructor__ - it is called in all childern classes. 
+
+    It takes the following arguments:
+    - `maxiter`: maximum number of iterations for the optimizer,
+    - `seed`: Seed for RandomState.
+
+    In constructor the attributes:
+    - `maxiter`,
+    - `seed`,
+    - `calc_weights` (used to fit weight, non-linear least squares (LS) method is used as a default)
+    are set.
+
+2. __fit__ - it is inherited by all children classes. It calls 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 and calculate optimised weights for the given relaxation times.
+    It takes the following arguments:
+    - `func`: objective function to be optimized,
+    - `lb`: the lower bounds of the design variable(s),
+    - `ub`: the upper bounds of the design variable(s),
+    - `funckwargs`: optional arguments takien by objective function.
+
+3. __cost_function__ - it is inherited by all children classes. It calculate cost function as the average error between the actual and the approximated electric permittivity (sum of real and imaginary part).
+    It takes the following arguments:
+    - `x`: the logarithm with base 10 of relaxation times of the Debyes poles,
+    - `rl`: real parts of chosen relaxation function for given frequency points,
+    - `im`: imaginary parts of chosen relaxation function for given frequency points,
+    - `freq`: the frequencies vector for defined grid.
+
+4. __calc_relaxation_times__ - it find an optimal set of relaxation times that minimise an objective function using appropriate optimization procedure.
+    It takes the following arguments:
+    - `func`: objective function to be optimized,
+    - `lb`: the lower bounds of the design variable(s),
+    - `ub`: the upper bounds of the design variable(s),
+    - `funckwargs`: optional arguments takien by objective function.
+
+Each new class of optimizer should:
+- define constructor with appropriate arguments,
+- overload __calc_relaxation_times__ method (and optional define __calc_weights__ function in case of hybrid optimization procedure).

user_libs/DebyeFit/docs/relaxation.md

@@ -0,0 +1,68 @@
+# Relaxation classes for multi-Debye fitting
+
+This class is designed for modelling different relaxation functions, like Havriliak-Negami (```Class HavriliakNegami```), Jonsher (```Class Jonsher```), Complex Refractive Index Mixing (```Class CRIM```) models, and arbitrary dielectric data derived experimentally or calculated using some other function (```Class Rawdata```).
+
+Supported relaxation classes:
+- [x] Havriliak-Negami,
+- [x] Jonsher,
+- [x] Complex Refractive Index Mixing,
+- [x] Experimental data,
+
+## Methods
+
+1. __constructor__ - it is called in all childern classes, creates Relaxation function object for complex material.
+
+    It takes the following arguments:
+    - ``sigma`` is a conductivity (Siemens/metre),
+    - ``mu`` is a relative permeability,
+    - ``mu_sigma`` is a magnetic loss (Ohms/metre),
+    - ``material_name`` is definition of material name,
+    - ``number_of_debye_poles`` is choosen number of Debye poles,
+    - ``f_n`` is choosen number of frequences,
+    - ``plot`` is a switch to turn on the plotting,
+    - ``save`` is a switch to turn on the saving final material properties,
+    - ``optimizer`` is a choosen optimizer to fit model to dielectric data,
+    - ``optimizer_options`` is a dict for options of choosen optimizer.
+
+    Additional parameters:
+    - ``rl`` calculated real part of chosen relaxation function for given frequency points,
+    - ``im`` calculated imaginary part of chosen relaxation function for given frequency points.
+
+2. __set_freq__ - it is inherited by all childern classes, interpolates frequency vector using n equally logarithmicaly spaced frequencies.
+    It takes the following arguments:
+    - `f_min_`: first bound of the frequency range used to approximate the given function (Hz),
+    - `f_max`: second bound of the frequency range used to approximate the given function (Hz),
+    - `f_n`: the number of frequency points in frequency grid (Default: 50).
+
+3. __run__ - it is inherited by all childern classes, solves the problem described by the given relaxation function (main operational method).
+    It consists of following steps:
+    1) Check the validity of the inputs using ```check_inputs``` method.
+    2) Print information about chosen approximation settings using ```print_info``` method.
+    3) Calculate both real and imaginary parts using ```calculation``` method, and then set ```self.rl``` and ```self.im``` properties.
+    4) Calling the main optimisation module using ```optimize``` method and calculate error based on ```error``` method.
+       a) [OPTIONAL] If number of debye poles is set to -1 optimization procedure is repeted until percentage error goes smaller than 5% or 20 Debye poles is reached.
+    5) Print the results in gprMax format style using ```print_output``` method.
+    6) [OPTIONAL] Save results in gprMax style using ```save_result``` method.
+    7) [OPTIONAL] Plot the actual and the approximate dielectric properties using ```plot_result``` method.
+
+4. __check_inputs__ - it is called in all childern classes, finds an optimal set of relaxation times that minimise an objective function using appropriate optimization procedure.
+
+5. __calculation__ - it is inherited by all childern classes, should be definied in all new children classes, approximates the given relaxation function.
+
+6. __print_info__ - it is inherited by all childern classes, prints readable string of parameters for given approximation settings.
+
+7. __optimize__ - it is inherited by all childern classes, calls the main optimisation module with defined lower and upper boundaries of search.
+
+8. __print_output__ - it is inherited by all childern classes, prints out the resulting Debye parameters in a gprMax format.
+
+9. __plot_result__ - it is inherited by all childern classes, plots the actual and the approximated electric permittivity, along with relative error for real and imaginary parts using a semilogarithm X axes.
+
+10. __save_result__ - it is inherited by all childern classes, saves the resulting Debye parameters in a gprMax format.
+
+11. __error__ - it is inherited by all childern classes, calculates the average fractional error separately for relative permittivity (real part) and conductivity (imaginary part).
+
+
+Each new class of relaxation object should:
+- define constructor with appropriate arguments,
+- define __check_inputs__ method to check relaxation class specific parameters,
+- overload __calculation__ method.

user_libs/DebyeFit/examples/Readme.md

@@ -0,0 +1,33 @@
+# Getting started
+
+Please see provided examples for the basic usage of ``DebyeFit`` module. We provide jupyter tutorials, and full guidance for quick run with existing relaxation functions and optimizers:
+
+* ```example_DebyeFitting.ipynb```: simple cases of using all available implemented relaxation functions,
+* ```example_BiologicalTissues.ipynb```: simple cases of using Cole-Cole function for biological tissues,
+* ```example_ColeCole.py```: simple cases of using Cole-Cole function in case of 3, 5 and automatically chosen number of Debye poles.
+
+Main usage of the specific relaxation fucntion based on creation of choosen relaxation model and then calling run method.
+
+```python
+
+    # set Havrilak-Negami function with initial parameters
+    setup = HavriliakNegami(f_min=1e4, f_max=1e11,
+                            alpha=0.3, beta=1,
+                            e_inf=3.4, de=2.7, tau_0=.8e-10,
+                            sigma=0.45e-3, mu=1, mu_sigma=0,
+                            material_name="dry_sand", f_n=100,
+                            plot=True, save=False,
+                            number_of_debye_poles=3,
+                            optimizer_options={'swarmsize':30,
+                                               'maxiter':100,
+                                               'omega':0.5,
+                                               'phip':1.4,
+                                               'phig':1.4,
+                                               'minstep':1e-8,
+                                               'minfun':1e-8,
+                                               'seed': 111,
+                                               'pflag': True})
+    # run optimization
+    setup.run()
+
+```

user_libs/DebyeFit/optimization.py

@@ -36,7 +36,7 @@ class Optimizer(object):
     def __init__(self, maxiter=1000, seed=None):
         self.maxiter = maxiter
         self.seed = seed
-        self.calc_weights = DLS
+        self.calc_weights = None
 
     def fit(self, func, lb, ub, funckwargs={}):
         """Call the optimization function that tries to find an optimal set
@@ -64,10 +64,19 @@ class Optimizer(object):
         np.random.seed(self.seed)
         # find the relaxation frequencies using choosen optimization alghoritm
         tau, _ = self.calc_relaxation_times(func, lb, ub, funckwargs)
-        # find the weights using a damped least squares
+        # find the weights using a calc_weights method
+        if self.calc_weights is None:
+            raise NotImplementedError()
         _, _, weights, ee, rl_exp, im_exp = self.calc_weights(tau, **funckwargs)
         return tau, weights, ee, rl_exp, im_exp
 
+    def calc_relaxation_times(self):
+        """ Optimisation method that tries to find an optimal set
+        of relaxation times that minimise the error
+        between the actual and the approximated electric permittivity.
+        """
+        raise NotImplementedError()
+
     @staticmethod
     def cost_function(x, rl, im, freq):
         """
@@ -88,28 +97,6 @@ class Optimizer(object):
         cost_i, cost_r, _, _, _, _ = DLS(x, rl, im, freq)
         return cost_i + cost_r
 
-    @staticmethod
-    def plot(x, y):
-        """
-        Dynamically plot the error as the optimisation takes place.
-
-        Args:
-            x (array): The number of current iterations.
-            y (array): The objective value at for all x points.
-        """
-        # it clears an axes
-        plt.cla()
-        plt.plot(x, y, "b-", linewidth=1.0)
-        plt.ylim(min(y) - 0.1 * min(y),
-                 max(y) + 0.1 * max(y))
-        plt.xlim(min(x) - 0.1, max(x) + 0.1)
-        plt.grid(b=True, which="major", color="k",
-                 linewidth=0.2, linestyle="--")
-        plt.suptitle("Debye fitting process")
-        plt.xlabel("Iteration")
-        plt.ylabel("Average Error")
-        plt.pause(0.0001)
-
 
 class PSO_DLS(Optimizer):
     """ Create hybrid Particle Swarm-Damped Least Squares optimisation
@@ -151,6 +138,7 @@ class PSO_DLS(Optimizer):
         self.minstep = minstep
         self.minfun = minfun
         self.pflag = pflag
+        self.calc_weights = DLS
 
     def calc_relaxation_times(self, func, lb, ub, funckwargs={}):
         """
@@ -265,8 +253,30 @@ class PSO_DLS(Optimizer):
                 PSO_DLS.plot(xpp, ypp)
         return g, fg
 
+    @staticmethod
+    def plot(x, y):
+        """
+        Dynamically plot the error as the optimisation takes place.
 
-class DA(Optimizer):
+        Args:
+            x (array): The number of current iterations.
+            y (array): The objective value at for all x points.
+        """
+        # it clears an axes
+        plt.cla()
+        plt.plot(x, y, "b-", linewidth=1.0)
+        plt.ylim(min(y) - 0.1 * min(y),
+                 max(y) + 0.1 * max(y))
+        plt.xlim(min(x) - 0.1, max(x) + 0.1)
+        plt.grid(b=True, which="major", color="k",
+                 linewidth=0.2, linestyle="--")
+        plt.suptitle("Debye fitting process")
+        plt.xlabel("Iteration")
+        plt.ylabel("Average Error")
+        plt.pause(0.0001)
+
+
+class DA_DLS(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:
@@ -278,7 +288,7 @@ class DA(Optimizer):
                  restart_temp_ratio=2e-05, visit=2.62, accept=-5.0,
                  maxfun=1e7, no_local_search=False,
                  callback=None, x0=None, seed=None):
-        super(DA, self).__init__(maxiter, seed)
+        super(DA_DLS, self).__init__(maxiter, seed)
         self.local_search_options = local_search_options
         self.initial_temp = initial_temp
         self.restart_temp_ratio = restart_temp_ratio
@@ -288,6 +298,7 @@ class DA(Optimizer):
         self.no_local_search = no_local_search
         self.callback = callback
         self.x0 = x0
+        self.calc_weights = DLS
 
     def calc_relaxation_times(self, func, lb, ub, funckwargs={}):
         """
@@ -327,7 +338,7 @@ class DA(Optimizer):
         return result.x, result.fun
 
 
-class DE(Optimizer):
+class DE_DLS(Optimizer):
     """
     Create Differential Evolution-Damped Least Squares object with predefined parameters.
     The current class is a modified edition of the scipy.optimize
@@ -340,7 +351,7 @@ class DE(Optimizer):
                  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)
+        super(DE_DLS, self).__init__(maxiter, seed)
         self.strategy = strategy
         self.popsize = popsize
         self.tol = tol
@@ -354,6 +365,7 @@ class DE(Optimizer):
         self.updating = updating
         self.workers = workers
         self.constraints = constraints
+        self.calc_weights = DLS
 
     def calc_relaxation_times(self, func, lb, ub, funckwargs={}):
         """