Add plotting and bootstrapped confidence intervals

Included matplotlib in the requirements.txt

Bootstrapping is done using the calculated standard deviations from fitting.
A plotting.py file includes a plot_nyquist function for simple, formatted matplotlib Nyquist plots

eisfit/circuits.py

@@ -1,4 +1,6 @@
 from .fitting import circuit_fit, computeCircuit, calculateCircuitLength
+from .plotting import plot_nyquist
+import matplotlib.pyplot as plt
 import numpy as np
 
 
@@ -58,10 +60,13 @@ def fit(self, frequencies, impedance):
             'mismatch in length of input frequencies and impedances'
 
         if self.initial_guess is not None:
-            self.parameters_, _ = circuit_fit(frequencies, impedance,
-                                              self.circuit, self.initial_guess,
-                                              self.algorithm,
-                                              bounds=self.bounds)
+            parameters, conf = circuit_fit(frequencies, impedance,
+                                           self.circuit, self.initial_guess,
+                                           self.algorithm,
+                                           bounds=self.bounds)
+            self.parameters_ = parameters
+            if conf is not None:
+                self.conf_ = conf
         else:
             # TODO auto calc guess
             raise ValueError('no initial guess supplied')
@@ -134,8 +139,78 @@ def __str__(self):
         to_print += '\n-------------------------------\n'
         return to_print
 
+    def plot(self, f_data=None, Z_data=None, CI=True):
+        """ a convenience method for plotting Nyquist plots
+
+
+        Parameters
+        ----------
+        f_data: np.array of type float
+            Frequencies of input data (for Bode plots)
+        Z_data: np.array of type complex
+            Impedance data to plot
+        CI: boolean
+            Include bootstrapped confidence intervals in plot
+
+        Returns
+        -------
+        ax: matplotlib.axes
+            axes of the created nyquist plot
+        """
+
+        fig, ax = plt.subplots(figsize=(5, 5))
+
+        if Z_data is not None:
+            ax = plot_nyquist(ax, f_data, Z_data)
+
+        if self._is_fit():
+
+            if f_data is not None:
+                f_pred = f_data
+            else:
+                f_pred = np.logspace(5, -3)
+
+            Z_fit = self.predict(f_pred)
+            ax = plot_nyquist(ax, f_data, Z_fit, fit=True)
+
+            if CI:
+                N = 1000
+                n = len(self.parameters_)
+                f_pred = np.logspace(np.log10(min(f_data)),
+                                     np.log10(max(f_data)),
+                                     num=100)
+
+                params = self.parameters_
+                confs = self.conf_
+
+                full_range = np.ndarray(shape=(N, len(f_pred)), dtype=complex)
+                for i in range(N):
+                    self.parameters_ = params + \
+                                        confs*np.random.uniform(-2, 2, size=n)
+
+                    full_range[i, :] = self.predict(f_pred)
+
+                self.parameters_ = params
+
+                min_Z = []
+                max_Z = []
+                for x in np.real(Z_fit):
+                    ys = []
+                    for run in full_range:
+                        ind = np.argmin(np.abs(run.real - x))
+                        ys.append(run[ind].imag)
+
+                    min_Z.append(x + 1j*min(ys))
+                    max_Z.append(x + 1j*max(ys))
+
+                ax.fill_between(np.real(min_Z), -np.imag(min_Z),
+                                -np.imag(max_Z), alpha='.2')
+
+        plt.show()
+
 
 class Randles(BaseCircuit):
+    """ A Randles circuit model class """
     def __init__(self, CPE=False, **kwargs):
         """ Constructor for the Randles' circuit class
 

eisfit/fitting.py

@@ -73,7 +73,7 @@ def circuit_fit(frequencies, impedances, circuit,
             for i, __ in enumerate(covar):
                 p_error.append(np.absolute(p_cov[i][i])**0.5)
         else:
-            p_error = len(p_values)*[-1]
+            p_error = None
     elif algorithm in ['SLSQP', 'L-BFGS-B', 'TNC']:
         if bounds is None:
             bounds = []
@@ -94,7 +94,7 @@ def circuit_fit(frequencies, impedances, circuit,
                        method=algorithm)
         p_values = res.x
         covar = None
-        p_error = len(p_values)*[-1]
+        p_error = None
 
     return p_values, p_error
 

eisfit/plotting.py

@@ -0,0 +1,52 @@
+import numpy as np
+from matplotlib.ticker import ScalarFormatter
+
+
+class FixedOrderFormatter(ScalarFormatter):
+    """Formats axis ticks using scientific notation with a constant order of
+    magnitude"""
+    def __init__(self, order_of_mag=0, useOffset=True, useMathText=True):
+        self._order_of_mag = order_of_mag
+        ScalarFormatter.__init__(self, useOffset=useOffset,
+                                 useMathText=useMathText)
+
+    def _set_orderOfMagnitude(self, range):
+        """Over-riding this to avoid having orderOfMagnitude reset elsewhere"""
+        self.orderOfMagnitude = self._order_of_mag
+
+
+def plot_nyquist(ax, freq, Z, fit=False):
+
+    if fit:
+        fmt = '.-'
+    else:
+        fmt = 'o'
+
+    ax.plot(np.real(Z), -np.imag(Z), fmt, lw=3)
+
+    # Make the axes square
+    ax.set_aspect('equal')
+
+    # Set the labels to -imaginary vs real
+    ax.set_xlabel('$Z_{1}^{\prime}(\omega)$ $[m\Omega]$', fontsize=20)
+    ax.set_ylabel('$-Z_{1}^{\prime\prime}(\omega)$ $[m\Omega]$', fontsize=20)
+
+    # Make the tick labels larger
+    ax.tick_params(axis='both', which='major', labelsize=14)
+
+    # Change the number of labels on each axis to five
+    ax.locator_params(axis='x', nbins=5, tight=True)
+    ax.locator_params(axis='y', nbins=5, tight=True)
+
+    # Add a light grid
+    ax.grid(b=True, which='major', axis='both', alpha=.5)
+
+    # Change axis units to 10^-3 and resize the offset text
+    ax.xaxis.set_major_formatter(FixedOrderFormatter(-3))
+    ax.yaxis.set_major_formatter(FixedOrderFormatter(-3))
+    y_offset = ax.yaxis.get_offset_text()
+    y_offset.set_size(18)
+    t = ax.xaxis.get_offset_text()
+    t.set_size(18)
+
+    return ax

requirements.txt

@@ -1,4 +1,5 @@
 certifi==2018.4.16
+matplotlib==2.2.2
 nose==1.3.7
 numpy==1.14.2
 numpydoc==0.8.0