From dc9836128d38955b586ac8cf9855465cee082568 Mon Sep 17 00:00:00 2001 From: Jake Anderson Date: Thu, 5 Sep 2024 10:36:16 -0600 Subject: [PATCH 1/2] use execution graph and add additional fitting methods --- impedance/models/circuits/circuits.py | 302 ++++++------ impedance/models/circuits/elements.py | 37 +- impedance/models/circuits/fitting.py | 563 ++++++++++++----------- impedance/tests/test_circuit_elements.py | 22 +- impedance/tests/test_circuits.py | 16 +- impedance/tests/test_fitting.py | 380 +++++++++------ requirements.txt | 2 + setup.py | 3 +- 8 files changed, 734 insertions(+), 591 deletions(-) diff --git a/impedance/models/circuits/circuits.py b/impedance/models/circuits/circuits.py index b670d19..a801378 100644 --- a/impedance/models/circuits/circuits.py +++ b/impedance/models/circuits/circuits.py @@ -1,7 +1,15 @@ -from .fitting import circuit_fit, buildCircuit -from .fitting import calculateCircuitLength, check_and_eval +from .fitting import ( + circuit_fit, + CircuitGraph, + calculateCircuitLength, + extract_circuit_elements, +) from impedance.visualization import plot_altair, plot_bode, plot_nyquist -from .elements import circuit_elements, get_element_from_name +from .elements import ( + circuit_elements, + get_element_from_name, + get_element_parameter_names_and_units_from_name, +) import json import matplotlib.pyplot as plt @@ -10,13 +18,17 @@ class BaseCircuit: - """ Base class for equivalent circuit models """ - def __init__(self, initial_guess=[], constants=None, name=None): - """ Base constructor for any equivalent circuit model + """Base class for equivalent circuit models""" + + def __init__(self, circuit="", initial_guess=None, constants=None, name=None): + """Base constructor for any equivalent circuit model Parameters ---------- - initial_guess: numpy array + circuit: str, optional + String describing this circuit + + initial_guess: numpy array, optional Initial guess of the circuit values constants : dict, optional @@ -26,12 +38,13 @@ def __init__(self, initial_guess=[], constants=None, name=None): name : str, optional Name for the circuit """ + if initial_guess is None: + initial_guess = [] # if supplied, check that initial_guess is valid and store initial_guess = [x for x in initial_guess if x is not None] - for i in initial_guess: - if not isinstance(i, (float, int, np.int32, np.float64)): - raise TypeError(f'value {i} in initial_guess is not a number') + if not np.all(np.isfinite(initial_guess)): + raise TypeError(f"initial_guess must be numeric or None ({initial_guess})") # initalize class attributes self.initial_guess = initial_guess @@ -44,6 +57,21 @@ def __init__(self, initial_guess=[], constants=None, name=None): # initialize fit parameters and confidence intervals self.parameters_ = None self.conf_ = None + self.cov_ = None + + self.circuit = CircuitGraph._whitespce.sub("", circuit) + + circuit_len = calculateCircuitLength(self.circuit) + + if len(self.initial_guess) + len(self.constants) != circuit_len: + raise ValueError( + "The number of initial guesses " + + f"({len(self.initial_guess)}) + " + + "the number of constants " + + f"({len(self.constants)})" + + " must be equal to " + + f"the circuit length ({circuit_len})" + ) def __eq__(self, other): if self.__class__ == other.__class__: @@ -55,11 +83,12 @@ def __eq__(self, other): matches.append(value == other.__dict__[key]) return np.array(matches).all() else: - raise TypeError('Comparing object is not of the same type.') + raise TypeError("Comparing object is not of the same type.") - def fit(self, frequencies, impedance, bounds=None, - weight_by_modulus=False, **kwargs): - """ Fit the circuit model + def fit( + self, frequencies, impedance, bounds=None, weight_by_modulus=False, **kwargs + ): + """Fit the circuit model Parameters ---------- @@ -90,36 +119,42 @@ def fit(self, frequencies, impedance, bounds=None, self: returns an instance of self """ - frequencies = np.array(frequencies, dtype=float) - impedance = np.array(impedance, dtype=complex) + frequencies = np.asarray(frequencies, dtype=float) + impedance = np.asarray(impedance, dtype=complex) if len(frequencies) != len(impedance): - raise TypeError('length of frequencies and impedance do not match') + raise TypeError("length of frequencies and impedance do not match") if self.initial_guess != []: - parameters, conf = circuit_fit(frequencies, impedance, - self.circuit, self.initial_guess, - constants=self.constants, - bounds=bounds, - weight_by_modulus=weight_by_modulus, - **kwargs) + parameters, conf, cov = circuit_fit( + frequencies, + impedance, + self.circuit, + self.initial_guess, + constants=self.constants, + bounds=bounds, + weight_by_modulus=weight_by_modulus, + return_covariance=True, + **kwargs, + ) self.parameters_ = parameters if conf is not None: self.conf_ = conf + self.cov_ = cov else: - raise ValueError('No initial guess supplied') + raise ValueError("No initial guess supplied") return self def _is_fit(self): - """ check if model has been fit (parameters_ is not None) """ + """check if model has been fit (parameters_ is not None)""" if self.parameters_ is not None: return True else: return False def predict(self, frequencies, use_initial=False): - """ Predict impedance using an equivalent circuit model + """Predict impedance using an equivalent circuit model Parameters ---------- @@ -134,81 +169,65 @@ def predict(self, frequencies, use_initial=False): Predicted impedance at each frequency """ frequencies = np.array(frequencies, dtype=float) - + cg = CircuitGraph(self.circuit, self.constants) if self._is_fit() and not use_initial: - return eval(buildCircuit(self.circuit, frequencies, - *self.parameters_, - constants=self.constants, eval_string='', - index=0)[0], - circuit_elements) + return cg.compute(frequencies, *self.parameters_) else: warnings.warn("Simulating circuit based on initial parameters") - return eval(buildCircuit(self.circuit, frequencies, - *self.initial_guess, - constants=self.constants, eval_string='', - index=0)[0], - circuit_elements) + return cg.compute(frequencies, *self.initial_guess) def get_param_names(self): - """ Converts circuit string to names and units """ + """Converts circuit string to names and units""" # parse the element names from the circuit string - names = self.circuit.replace('p', '').replace('(', '').replace(')', '') - names = names.replace(',', '-').replace(' ', '').split('-') + names = extract_circuit_elements(self.circuit) full_names, all_units = [], [] for name in names: - elem = get_element_from_name(name) - num_params = check_and_eval(elem).num_params - units = check_and_eval(elem).units - if num_params > 1: - for j in range(num_params): - full_name = '{}_{}'.format(name, j) - if full_name not in self.constants.keys(): - full_names.append(full_name) - all_units.append(units[j]) - else: - if name not in self.constants.keys(): - full_names.append(name) - all_units.append(units[0]) + p_names, p_units = get_element_parameter_names_and_units_from_name(name) + full_names.extend([n for n in p_names if n not in self.constants]) + all_units.extend( + [u for u, n in zip(p_units, p_names) if n not in self.constants] + ) return full_names, all_units def __str__(self): - """ Defines the pretty printing of the circuit""" + """Defines the pretty printing of the circuit""" - to_print = '\n' + to_print = "\n" if self.name is not None: - to_print += 'Name: {}\n'.format(self.name) - to_print += 'Circuit string: {}\n'.format(self.circuit) + to_print += "Name: {}\n".format(self.name) + to_print += "Circuit string: {}\n".format(self.circuit) to_print += "Fit: {}\n".format(self._is_fit()) if len(self.constants) > 0: - to_print += '\nConstants:\n' + to_print += "\nConstants:\n" for name, value in self.constants.items(): elem = get_element_from_name(name) - units = check_and_eval(elem).units - if '_' in name and len(units) > 1: - unit = units[int(name.split('_')[-1])] + units = circuit_elements[elem].units + name_parts = name.split("_") + if len(name_parts) > 1: + unit = units[int(name_parts[-1])] else: unit = units[0] - to_print += ' {:>5} = {:.2e} [{}]\n'.format(name, value, unit) + to_print += " {:>5} = {:.2e} [{}]\n".format(name, value, unit) names, units = self.get_param_names() - to_print += '\nInitial guesses:\n' + to_print += "\nInitial guesses:\n" for name, unit, param in zip(names, units, self.initial_guess): - to_print += ' {:>5} = {:.2e} [{}]\n'.format(name, param, unit) + to_print += " {:>5} = {:.2e} [{}]\n".format(name, param, unit) if self._is_fit(): params, confs = self.parameters_, self.conf_ - to_print += '\nFit parameters:\n' + to_print += "\nFit parameters:\n" for name, unit, param, conf in zip(names, units, params, confs): - to_print += ' {:>5} = {:.2e}'.format(name, param) - to_print += ' (+/- {:.2e}) [{}]\n'.format(conf, unit) + to_print += " {:>5} = {:.2e}".format(name, param) + to_print += " (+/- {:.2e}) [{}]\n".format(conf, unit) return to_print - def plot(self, ax=None, f_data=None, Z_data=None, kind='altair', **kwargs): - """ visualizes the model and optional data as a nyquist, + def plot(self, ax=None, f_data=None, Z_data=None, kind="altair", **kwargs): + """visualizes the model and optional data as a nyquist, bode, or altair (interactive) plots Parameters @@ -236,12 +255,12 @@ def plot(self, ax=None, f_data=None, Z_data=None, kind='altair', **kwargs): axes of the created nyquist plot """ - if kind == 'nyquist': + if kind == "nyquist": if ax is None: _, ax = plt.subplots(figsize=(5, 5)) if Z_data is not None: - ax = plot_nyquist(Z_data, ls='', marker='s', ax=ax, **kwargs) + ax = plot_nyquist(Z_data, ls="", marker="s", ax=ax, **kwargs) if self._is_fit(): if f_data is not None: @@ -250,9 +269,9 @@ def plot(self, ax=None, f_data=None, Z_data=None, kind='altair', **kwargs): f_pred = np.logspace(5, -3) Z_fit = self.predict(f_pred) - ax = plot_nyquist(Z_fit, ls='-', marker='', ax=ax, **kwargs) + ax = plot_nyquist(Z_fit, ls="-", marker="", ax=ax, **kwargs) return ax - elif kind == 'bode': + elif kind == "bode": if ax is None: _, ax = plt.subplots(nrows=2, figsize=(5, 5)) @@ -263,21 +282,20 @@ def plot(self, ax=None, f_data=None, Z_data=None, kind='altair', **kwargs): if Z_data is not None: if f_data is None: - raise ValueError('f_data must be specified if' + - ' Z_data for a Bode plot') - ax = plot_bode(f_data, Z_data, ls='', marker='s', - axes=ax, **kwargs) + raise ValueError( + "f_data must be specified if" + " Z_data for a Bode plot" + ) + ax = plot_bode(f_data, Z_data, ls="", marker="s", axes=ax, **kwargs) if self._is_fit(): Z_fit = self.predict(f_pred) - ax = plot_bode(f_pred, Z_fit, ls='-', marker='', - axes=ax, **kwargs) + ax = plot_bode(f_pred, Z_fit, ls="-", marker="", axes=ax, **kwargs) return ax - elif kind == 'altair': + elif kind == "altair": plot_dict = {} if Z_data is not None and f_data is not None: - plot_dict['data'] = {'f': f_data, 'Z': Z_data} + plot_dict["data"] = {"f": f_data, "Z": Z_data} if self._is_fit(): if f_data is not None: @@ -289,17 +307,19 @@ def plot(self, ax=None, f_data=None, Z_data=None, kind='altair', **kwargs): if self.name is not None: name = self.name else: - name = 'fit' - plot_dict[name] = {'f': f_pred, 'Z': Z_fit, 'fmt': '-'} + name = "fit" + plot_dict[name] = {"f": f_pred, "Z": Z_fit, "fmt": "-"} chart = plot_altair(plot_dict, **kwargs) return chart else: - raise ValueError("Kind must be one of 'altair'," + - f"'nyquist', or 'bode' (received {kind})") + raise ValueError( + "Kind must be one of 'altair'," + + f"'nyquist', or 'bode' (received {kind})" + ) def save(self, filepath): - """ Exports a model to JSON + """Exports a model to JSON Parameters ---------- @@ -315,27 +335,34 @@ def save(self, filepath): if self._is_fit(): parameters_ = list(self.parameters_) model_conf_ = list(self.conf_) - - data_dict = {"Name": model_name, - "Circuit String": model_string, - "Initial Guess": initial_guess, - "Constants": self.constants, - "Fit": True, - "Parameters": parameters_, - "Confidence": model_conf_, - } + model_cov_ = None + if self.cov_ is not None: + model_cov_ = self.cov_.tolist() + + data_dict = { + "Name": model_name, + "Circuit String": model_string, + "Initial Guess": initial_guess, + "Constants": self.constants, + "Fit": True, + "Parameters": parameters_, + "Confidence": model_conf_, + "Covariance": model_cov_, + } else: - data_dict = {"Name": model_name, - "Circuit String": model_string, - "Initial Guess": initial_guess, - "Constants": self.constants, - "Fit": False} - - with open(filepath, 'w') as f: + data_dict = { + "Name": model_name, + "Circuit String": model_string, + "Initial Guess": initial_guess, + "Constants": self.constants, + "Fit": False, + } + + with open(filepath, "w") as f: json.dump(data_dict, f) def load(self, filepath, fitted_as_initial=False): - """ Imports a model from JSON + """Imports a model from JSON Parameters ---------- @@ -350,7 +377,7 @@ def load(self, filepath, fitted_as_initial=False): fitted parameters as a completed model """ - json_data_file = open(filepath, 'r') + json_data_file = open(filepath, "r") json_data = json.load(json_data_file) model_name = json_data["Name"] @@ -366,16 +393,19 @@ def load(self, filepath, fitted_as_initial=False): if json_data["Fit"]: if fitted_as_initial: - self.initial_guess = np.array(json_data['Parameters']) + self.initial_guess = np.array(json_data["Parameters"]) else: self.parameters_ = np.array(json_data["Parameters"]) self.conf_ = np.array(json_data["Confidence"]) + if json_data.get("Covariance", None) is not None: + self.cov_ = np.array(json_data["Covariance"]) class Randles(BaseCircuit): - """ A Randles circuit model class """ + """A Randles circuit model class""" + def __init__(self, CPE=False, **kwargs): - """ Constructor for the Randles' circuit class + """Constructor for the Randles' circuit class Parameters ---------- @@ -385,60 +415,18 @@ def __init__(self, CPE=False, **kwargs): CPE: boolean Use a constant phase element instead of a capacitor """ - super().__init__(**kwargs) - if CPE: - self.name = 'Randles w/ CPE' - self.circuit = 'R0-p(R1-Wo1,CPE1)' + name = "Randles w/ CPE" + circuit = "R0-p(R1-Wo1,CPE1)" else: - self.name = 'Randles' - self.circuit = 'R0-p(R1-Wo1,C1)' + name = "Randles" + circuit = "R0-p(R1-Wo1,C1)" - circuit_len = calculateCircuitLength(self.circuit) - - if len(self.initial_guess) + len(self.constants) != circuit_len: - raise ValueError('The number of initial guesses ' + - f'({len(self.initial_guess)}) + ' + - 'the number of constants ' + - f'({len(self.constants)})' + - ' must be equal to ' + - f'the circuit length ({circuit_len})') - - -class CustomCircuit(BaseCircuit): - def __init__(self, circuit='', **kwargs): - """ Constructor for a customizable equivalent circuit model - - Parameters - ---------- - initial_guess: numpy array - Initial guess of the circuit values - - circuit: string - A string that should be interpreted as an equivalent circuit - - Notes - ----- - A custom circuit is defined as a string comprised of elements in series - (separated by a `-`) and elements in parallel (grouped as (x,y)). - Each element can be appended with an integer (e.g. R0) or an underscore - and an integer (e.g. CPE_1) to make keeping track of multiple elements - of the same type easier. - - Example: - Randles circuit is given by 'R0-p(R1-Wo1,C1)' - - """ + kwargs["circuit"] = circuit + if "name" not in kwargs: + kwargs["name"] = name super().__init__(**kwargs) - self.circuit = circuit.replace(" ", "") - circuit_len = calculateCircuitLength(self.circuit) - if len(self.initial_guess) + len(self.constants) != circuit_len: - raise ValueError('The number of initial guesses ' + - f'({len(self.initial_guess)}) + ' + - 'the number of constants ' + - f'({len(self.constants)})' + - ' must be equal to ' + - f'the circuit length ({circuit_len})') +CustomCircuit = BaseCircuit diff --git a/impedance/models/circuits/elements.py b/impedance/models/circuits/elements.py index 643ff69..0a95b03 100644 --- a/impedance/models/circuits/elements.py +++ b/impedance/models/circuits/elements.py @@ -45,12 +45,6 @@ def wrapper(p, f): ) else: circuit_elements[func.__name__] = wrapper - # Adding numpy to circuit_elements for proper evaluation with - # numpy>=2.0.0 because the scalar representation was changed. - # "Scalars are now printed as np.float64(3.0) rather than just 3.0." - # https://numpy.org/doc/2.0/release/2.0.0-notes.html - # #representation-of-numpy-scalars-changed - circuit_elements["np"] = np return wrapper @@ -412,17 +406,20 @@ def get_element_from_name(name): def typeChecker(p, f, name, length): - assert isinstance(p, list), \ - "in {}, input must be of type list".format(name) - for i in p: - assert isinstance( - i, (float, int, np.int32, np.float64) - ), "in {}, value {} in {} is not a number".format(name, i, p) - for i in f: - assert isinstance( - i, (float, int, np.int32, np.float64) - ), "in {}, value {} in {} is not a number".format(name, i, f) - assert len(p) == length, "in {}, input list must be length {}".format( - name, length - ) - return + if not np.all(np.isreal(p)): + raise TypeError(f"In {name} all of the parameters are not real ({p})") + if not np.all(np.isreal(f)): + raise TypeError(f"In {name} all of the frequencies are not ({f})") + if len(p) != length: + raise TypeError(f"In {name} length of parameters is not {length} ") + +def get_element_parameter_names_and_units_from_name(name): + elem = get_element_from_name(name) + n_params = circuit_elements[elem].num_params + return [ + format_parameter_name(name, j, n_params) + for j in range(n_params) + ], circuit_elements[elem].units + +def format_parameter_name(name, j, n_params): + return f"{name}_{j}" if n_params > 1 else f"{name}" \ No newline at end of file diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index 5fe1c30..7bc1e25 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -1,12 +1,40 @@ -import warnings +import networkx as nx +import re import numpy as np -from scipy.linalg import inv -from scipy.optimize import curve_fit, basinhopping +import numdifftools as nd +from scipy.linalg import inv, norm +from scipy.optimize import curve_fit, dual_annealing, minimize, basinhopping, Bounds -from .elements import circuit_elements, get_element_from_name +from .elements import circuit_elements, get_element_from_name, format_parameter_name -ints = '0123456789' +ints = "0123456789" +_LBFGSB_OPTIONS = { + "ftol": 1e-14, + "gtol": 1e-13, + "maxfun": 10_000, +} + +_KWARGS = { + "method": "L-BFGS-B", + "options": _LBFGSB_OPTIONS, +} + + +class BoundsCheck(object): + def __init__(self, lb, ub): + self.lb = lb + self.ub = ub + + if np.any(self.lb > self.ub): + raise ValueError( + "the lower bound must be less than or equal to the upper bound." + ) + + def __call__(self, f_new, x_new, f_old, x_old): + return np.all( + np.logical_and(x_new >= self.lb, x_new <= self.ub) + ) and np.isfinite(f_new) def rmse(a, b): @@ -22,11 +50,41 @@ def rmse(a, b): """ n = len(a) - return np.linalg.norm(a - b) / np.sqrt(n) + return np.real(np.sqrt(sum_square_error(a, b) / n)) + + +def sum_square_error(a, b, w=1): + amb = a - b + return np.sum(amb * np.conjugate(amb) * w) + + +def objective_function(graph, f, Z, weight_by_modulus=False): + if weight_by_modulus: + weights = 1 / norm(Z) + else: + weights = 1.0 + + def cost_function(x): + """Short function to optimize over. + We want to minimize the square error between the model and the data. + + Parameters + ---------- + x : args + Parameters for optimization. + + Returns + ------- + function + Returns a function + """ + return sum_square_error(graph(f, *x), np.stack([Z.real, Z.imag]), weights) + + return cost_function def set_default_bounds(circuit, constants={}): - """ This function sets default bounds for optimization. + """This function sets default bounds for optimization. set_default_bounds sets bounds of 0 and np.inf for all parameters, except the CPE and La alphas which have an upper bound of 1. @@ -53,10 +111,10 @@ def set_default_bounds(circuit, constants={}): lower_bounds, upper_bounds = [], [] for elem in extracted_elements: raw_element = get_element_from_name(elem) - for i in range(check_and_eval(raw_element).num_params): - if elem in constants or elem + f'_{i}' in constants: + for i in range(circuit_elements[raw_element].num_params): + if elem in constants or elem + f"_{i}" in constants: continue - if raw_element in ['CPE', 'La'] and i == 1: + if raw_element in ["CPE", "La"] and i == 1: upper_bounds.append(1) else: upper_bounds.append(np.inf) @@ -66,11 +124,20 @@ def set_default_bounds(circuit, constants={}): return bounds -def circuit_fit(frequencies, impedances, circuit, initial_guess, constants={}, - bounds=None, weight_by_modulus=False, global_opt=False, - **kwargs): - - """ Main function for fitting an equivalent circuit to data. +def circuit_fit( + frequencies, + impedances, + circuit, + initial_guess, + constants={}, + bounds=None, + weight_by_modulus=False, + global_opt=False, + opt_method="curve_fit", + return_covariance=False, + **kwargs, +): + """Main function for fitting an equivalent circuit to data. By default, this function uses `scipy.optimize.curve_fit `_ @@ -113,8 +180,18 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants={}, If global optimization should be used (uses the basinhopping algorithm). Defaults to False + opt_method : str, optional + What optimization method to use defaults to "curve_fit" when global_opt + is False, but it could also be "minimize". + When global_opt is True the default is "basinhopping", but + "dual_annealing" is also a valid choice. + + return_covariance : bool, option + Return the covariance matrix as well as the errors. Defaults to False. + kwargs : - Keyword arguments passed to scipy.optimize.curve_fit or + Keyword arguments passed to scipy.optimize.curve_fit, + scipy.optimize.minimize, scipy.optimize.dual_annealing or scipy.optimize.basinhopping Returns @@ -123,247 +200,123 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants={}, best fit parameters for specified equivalent circuit p_errors : list of floats - one standard deviation error estimates for fit parameters + one standard deviation error estimates for fit parameters derived from + an estimate of the covariance matrix - Notes - --------- - Need to do a better job of handling errors in fitting. - Currently, an error of -1 is returned. + p_cov : (when return_covariance is True) numpy array of floats + estimated covariance matrix among the fit parameters """ - f = np.array(frequencies, dtype=float) - Z = np.array(impedances, dtype=complex) + f = np.asarray(frequencies, dtype=float) + Z = np.asarray(impedances, dtype=complex) + popt = None + perror = None + pcov = None # set upper and lower bounds on a per-element basis if bounds is None: bounds = set_default_bounds(circuit, constants=constants) + cg = CircuitGraph(circuit, constants) + sumsq_errors = objective_function(cg, f, Z, weight_by_modulus=weight_by_modulus) if not global_opt: - if 'maxfev' not in kwargs: - kwargs['maxfev'] = 1e5 - if 'ftol' not in kwargs: - kwargs['ftol'] = 1e-13 - - # weighting scheme for fitting - if weight_by_modulus: - abs_Z = np.abs(Z) - kwargs['sigma'] = np.hstack([abs_Z, abs_Z]) - - popt, pcov = curve_fit(wrapCircuit(circuit, constants), f, - np.hstack([Z.real, Z.imag]), - p0=initial_guess, bounds=bounds, **kwargs) - - # Calculate one standard deviation error estimates for fit parameters, - # defined as the square root of the diagonal of the covariance matrix. - # https://stackoverflow.com/a/52275674/5144795 - perror = np.sqrt(np.diag(pcov)) + if opt_method == "minimize": + if "options" not in kwargs: + kwargs["options"] = _LBFGSB_OPTIONS + bounds = [(lb, ub) for lb, ub in zip(*bounds)] + result = minimize( + sumsq_errors, + initial_guess, + bounds=bounds, + **kwargs, + ) + popt = result.x + pcov = result.get("hess_inv", None) + if (pcov is None) and ("hess" in result): + pcov = inv(result.hess) + try: + pcov = pcov.todense() + except AttributeError: + pass + else: + if "maxfev" not in kwargs: + kwargs["maxfev"] = 1e5 + if "ftol" not in kwargs: + kwargs["ftol"] = 1e-13 + + # weighting scheme for fitting + if weight_by_modulus: + abs_Z = np.abs(Z) + kwargs["sigma"] = np.concatenate([abs_Z, abs_Z]) + + popt, pcov = curve_fit( + cg.compute_long, + f, + np.concatenate([Z.real, Z.imag]), + p0=initial_guess, + bounds=bounds, + **kwargs, + ) else: - if 'seed' not in kwargs: - kwargs['seed'] = 0 - - def opt_function(x): - """ Short function for basinhopping to optimize over. - We want to minimize the RMSE between the fit and the data. - - Parameters - ---------- - x : args - Parameters for optimization. - - Returns - ------- - function - Returns a function (RMSE as a function of parameters). - """ - return rmse(wrapCircuit(circuit, constants)(f, *x), - np.hstack([Z.real, Z.imag])) - - class BasinhoppingBounds(object): - """ Adapted from the basinhopping documetation - https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.basinhopping.html - """ - - def __init__(self, xmin, xmax): - self.xmin = np.array(xmin) - self.xmax = np.array(xmax) - - def __call__(self, **kwargs): - x = kwargs['x_new'] - tmax = bool(np.all(x <= self.xmax)) - tmin = bool(np.all(x >= self.xmin)) - return tmax and tmin - - basinhopping_bounds = BasinhoppingBounds(xmin=bounds[0], - xmax=bounds[1]) - results = basinhopping(opt_function, x0=initial_guess, - accept_test=basinhopping_bounds, **kwargs) - popt = results.x - - # Calculate perror - jac = results.lowest_optimization_result['jac'][np.newaxis] - try: - # jacobian -> covariance - # https://stats.stackexchange.com/q/231868 - pcov = inv(np.dot(jac.T, jac)) * opt_function(popt) ** 2 - # covariance -> perror (one standard deviation - # error estimates for fit parameters) - perror = np.sqrt(np.diag(pcov)) - except (ValueError, np.linalg.LinAlgError): - warnings.warn('Failed to compute perror') - perror = None - - return popt, perror - - -def wrapCircuit(circuit, constants): - """ wraps function so we can pass the circuit string """ - def wrappedCircuit(frequencies, *parameters): - """ returns a stacked array of real and imaginary impedance - components - - Parameters - ---------- - circuit : string - constants : dict - parameters : list of floats - frequencies : list of floats - - Returns - ------- - array of floats - - """ - - x = eval(buildCircuit(circuit, frequencies, *parameters, - constants=constants, eval_string='', - index=0)[0], - circuit_elements) - y_real = np.real(x) - y_imag = np.imag(x) - - return np.hstack([y_real, y_imag]) - return wrappedCircuit - - -def buildCircuit(circuit, frequencies, *parameters, - constants=None, eval_string='', index=0): - """ recursive function that transforms a circuit, parameters, and - frequencies into a string that can be evaluated - - Parameters - ---------- - circuit: str - frequencies: list/tuple/array of floats - parameters: list/tuple/array of floats - constants: dict - - Returns - ------- - eval_string: str - Python expression for calculating the resulting fit - index: int - Tracks parameter index through recursive calling of the function - """ - - parameters = np.array(parameters).tolist() - frequencies = np.array(frequencies).tolist() - circuit = circuit.replace(' ', '') - - def parse_circuit(circuit, parallel=False, series=False): - """ Splits a circuit string by either dashes (series) or commas - (parallel) outside of any paranthesis. Removes any leading 'p(' - or trailing ')' when in parallel mode """ - - assert parallel != series, \ - 'Exactly one of parallel or series must be True' - - def count_parens(string): - return string.count('('), string.count(')') - - if parallel: - special = ',' - if circuit.endswith(')') and circuit.startswith('p('): - circuit = circuit[2:-1] - if series: - special = '-' - - split = circuit.split(special) - - result = [] - skipped = [] - for i, sub_str in enumerate(split): - if i not in skipped: - if '(' not in sub_str and ')' not in sub_str: - result.append(sub_str) - else: - open_parens, closed_parens = count_parens(sub_str) - if open_parens == closed_parens: - result.append(sub_str) - else: - uneven = True - while i < len(split) - 1 and uneven: - sub_str += special + split[i+1] - - open_parens, closed_parens = count_parens(sub_str) - uneven = open_parens != closed_parens - - i += 1 - skipped.append(i) - result.append(sub_str) - return result - - parallel = parse_circuit(circuit, parallel=True) - series = parse_circuit(circuit, series=True) - - if series is not None and len(series) > 1: - eval_string += "s([" - split = series - elif parallel is not None and len(parallel) > 1: - eval_string += "p([" - split = parallel - elif series == parallel: # only single element - split = series - - for i, elem in enumerate(split): - if ',' in elem or '-' in elem: - eval_string, index = buildCircuit(elem, frequencies, - *parameters, - constants=constants, - eval_string=eval_string, - index=index) + if "seed" not in kwargs: + kwargs["seed"] = 0 + + if opt_method == "dual_annealing": + if "local_search_options" not in kwargs: + kwargs["local_search_options"] = _KWARGS + kwargs["local_search_options"]["bounds"] = Bounds( + lb=bounds[0], ub=bounds[1] + ) + + bounds = (np.maximum(-1e9, bounds[0]), np.minimum(1e9, bounds[1])) + + results = dual_annealing( + sumsq_errors, + x0=initial_guess, + bounds=np.asarray(bounds).T, + **kwargs, + ) + popt = results.x + pcov = None else: - param_string = "" - raw_elem = get_element_from_name(elem) - elem_number = check_and_eval(raw_elem).num_params - param_list = [] - for j in range(elem_number): - if elem_number > 1: - current_elem = elem + '_{}'.format(j) - else: - current_elem = elem - - if current_elem in constants.keys(): - param_list.append(constants[current_elem]) - else: - param_list.append(parameters[index]) - index += 1 - - param_string += str(param_list) - new = raw_elem + '(' + param_string + ',' + str(frequencies) + ')' - eval_string += new - - if i == len(split) - 1: - if len(split) > 1: # do not add closing brackets if single element - eval_string += '])' - else: - eval_string += ',' + if "minimizer_kwargs" not in kwargs: + kwargs["minimizer_kwargs"] = _KWARGS + + results = basinhopping( + sumsq_errors, + x0=initial_guess, + accept_test=BoundsCheck(*bounds), + **kwargs, + ) + + popt = results.x + pcov = results.get("lowest_optimization_result", {}).get("hess_inv", None) + try: + pcov = pcov.todense() + except AttributeError: + pass + + if (pcov is None) or np.any(np.diag(pcov) < 0): + for m in ["central", "backward", "forward"]: + hess = nd.Hessian(sumsq_errors, method=m) + hmat = hess(popt) + pcov = inv(hmat) + if np.all(np.diag(pcov) >= 0): + break + pcov = None + + # Calculate one standard deviation error estimates for fit parameters, + # defined as the square root of the diagonal of the covariance matrix. + # https://stackoverflow.com/a/52275674/5144795 + if pcov is not None: + perror = np.sqrt(np.diag(pcov)) - return eval_string, index + return (popt, perror) if not return_covariance else (popt, perror, pcov) def extract_circuit_elements(circuit): - """ Extracts circuit elements from a circuit string. + """Extracts circuit elements from a circuit string. Parameters ---------- @@ -376,7 +329,7 @@ def extract_circuit_elements(circuit): list of extracted elements. """ - p_string = [x for x in circuit if x not in 'p(),-'] + p_string = [x for x in circuit if x not in "p(),-"] extracted_elements = [] current_element = [] length = len(p_string) @@ -385,18 +338,18 @@ def extract_circuit_elements(circuit): current_element.append(char) else: # min to prevent looking ahead past end of list - if p_string[min(i+1, length-1)] not in ints: + if p_string[min(i + 1, length - 1)] not in ints: current_element.append(char) - extracted_elements.append(''.join(current_element)) + extracted_elements.append("".join(current_element)) current_element = [] else: current_element.append(char) - extracted_elements.append(''.join(current_element)) + extracted_elements.append("".join(current_element)) return extracted_elements def calculateCircuitLength(circuit): - """ Calculates the number of elements in the circuit. + """Calculates the number of elements in the circuit. Parameters ---------- @@ -414,32 +367,108 @@ def calculateCircuitLength(circuit): extracted_elements = extract_circuit_elements(circuit) for elem in extracted_elements: raw_element = get_element_from_name(elem) - num_params = check_and_eval(raw_element).num_params + num_params = circuit_elements[raw_element].num_params length += num_params return length -def check_and_eval(element): - """ Checks if an element is valid, then evaluates it. +class CircuitGraph: + _parallel_block_expression = re.compile(r"p\([^\(\)]*\)") + _whitespce = re.compile(r"\s+") + + def __init__(self, circuit, constants=None): + self.circuit = self._whitespce.sub("", circuit) + + self.parse_circuit() + self.execution_order = list(nx.topological_sort(self.graph)) + self.constants = constants if constants is not None else dict() + + def parse_circuit(self): + self.snum = 1 + self.pnum = 1 + + parsing_circuit = self.circuit + + # determine all of the base elements, their functions and add them to the graph + element_name = extract_circuit_elements(parsing_circuit) + element_func = [ + circuit_elements[get_element_from_name(e)] for e in element_name + ] + self.graph = nx.DiGraph() + for e, f in zip(element_name, element_func): + self.graph.add_node(e, Z=f) + + # find unnested parallel blocks + pblocks = self._parallel_block_expression.findall(parsing_circuit) + while len(pblocks) > 0: + # add parallel blocks to the graph unnesting each time around the loop + for p in pblocks: + pelem = p[2:-1].split(",") + pnode = f"p{self.pnum}" + self.pnum += 1 + self.graph.add_node(pnode, Z=circuit_elements["p"]) + for elem in pelem: + elem = self.add_series_elements(elem) + self.graph.add_edge(elem, pnode) + parsing_circuit = parsing_circuit.replace(p, pnode) + pblocks = self._parallel_block_expression.findall(parsing_circuit) + + # pick up any top line series connections + self.add_series_elements(parsing_circuit) + + for layer, nodes in enumerate(nx.topological_generations(self.graph)): + for n in nodes: + self.graph.nodes[n]["layer"] = layer + + def add_series_elements(self, elem): + selem = elem.split("-") + if len(selem) > 1: + node = f"s{self.snum}" + self.snum += 1 + self.graph.add_node(node, Z=circuit_elements["s"]) + for n in selem: + self.graph.add_edge(n, node) + return node + + # if there isn't a series connection in elem just return it unchanged + return selem[0] + + def visualize_graph(self, **kwargs): + pos = nx.multipartite_layout(self.graph, subset_key="layer") + nx.draw_networkx(self.graph, pos=pos, **kwargs) + + def compute(self, f, *parameters): + node_results = {} + pindex = 0 + for node in self.execution_order: + Zfunc = self.graph.nodes[node]["Z"] + plist = [node_results[pred] for pred in self.graph.predecessors(node)] + if len(plist) < 1: + n_params = Zfunc.num_params + for j in range(n_params): + p_name = format_parameter_name(node, j, n_params) + if p_name in self.constants: + plist.append(self.constants[p_name]) + else: + plist.append(parameters[pindex]) + pindex += 1 + node_results[node] = Zfunc(plist, f) + else: + node_results[node] = Zfunc(plist) - Parameters - ---------- - element : str - Circuit element. + return np.squeeze(node_results[node]) - Raises - ------ - ValueError - Raised if an element is not in the list of allowed elements. + def __call__(self, f, *parameters): + Z = self.compute(f, *parameters) + return np.stack([Z.real, Z.imag]) - Returns - ------- - Evaluated element. + def compute_long(self, f, *parameters): + Z = self.compute(f, *parameters) + return np.concatenate([Z.real, Z.imag]) - """ - allowed_elements = circuit_elements.keys() - if element not in allowed_elements: - raise ValueError(f'{element} not in ' + - f'allowed elements ({allowed_elements})') - else: - return eval(element, circuit_elements) + def calculate_circuit_length(self): + n_params = [ + getattr(Zfunc, "num_params", 0) + for node, Zfunc in self.graph.nodes(data="Z") + ] + return np.sum(n_params) diff --git a/impedance/tests/test_circuit_elements.py b/impedance/tests/test_circuit_elements.py index cd1e3f2..7057c99 100644 --- a/impedance/tests/test_circuit_elements.py +++ b/impedance/tests/test_circuit_elements.py @@ -83,14 +83,20 @@ def test_each_element(): val = f(input_vals[:num_inputs], freqs) assert np.isclose(val, correct_vals[key]).all() - # check for typing: - with pytest.raises(AssertionError): - f = circuit_elements["R"] - f(1, 2) - - # test for handling more wrong inputs - with pytest.raises(AssertionError): - f(["hi"], ["yes", "hello"]) + # check for typing: + f = circuit_elements["R"] + with pytest.raises(TypeError): + f(1, 2) + + # test for handling more wrong inputs + with pytest.raises(TypeError): + f(["hi"], ["yes", "hello"]) + with pytest.raises(TypeError): + f([0.1, 0.2], [1]) + with pytest.raises(TypeError): + f([1 + 1j], [1]) + with pytest.raises(TypeError): + f([0.1], [1 + 1j]) # Test no overflow in T at high frequencies with warnings.catch_warnings(): diff --git a/impedance/tests/test_circuits.py b/impedance/tests/test_circuits.py index ec0e5c3..68251b7 100644 --- a/impedance/tests/test_circuits.py +++ b/impedance/tests/test_circuits.py @@ -16,30 +16,32 @@ def test_BaseCircuit(): - initial_guess = [0.01, 0.02, 50] + initial_guess = [] # __init__() # check initial_guess is loaded in correctly - base_circuit = BaseCircuit(initial_guess) + base_circuit = BaseCircuit(initial_guess=initial_guess) assert base_circuit.initial_guess == initial_guess # improper input_guess types raise an TypeError with pytest.raises(TypeError): r = BaseCircuit(initial_guess=['hi', 0.1]) + with pytest.raises(TypeError): + r = BaseCircuit(initial_guess=[0.1, np.inf]) # __eq__() # incorrect comparisons raise a TypeError with pytest.raises(TypeError): - r = BaseCircuit(initial_guess=[.01, .005, .1, .0001, 200]) + r = BaseCircuit(initial_guess=[]) r == 8 # fit() with pytest.raises(TypeError): - r = BaseCircuit(initial_guess=[.01, .005, .1, .0001, 200]) + r = BaseCircuit(initial_guess=[]) r.fit(np.array([42 + 42j]), []) # frequencies are complex with pytest.raises(TypeError): - r = BaseCircuit(initial_guess=[.01, .005, .1, .0001, 200]) + r = BaseCircuit(initial_guess=[]) r.fit(np.array([42, 4.2]), np.array([42 + 42j])) # mismatched lengths # plot() @@ -125,7 +127,7 @@ def test_CustomCircuit(): # check predictions from initial_guesses high_f = np.array([1e9]) - assert np.isclose(np.real(custom_circuit.predict(high_f)[0]), + assert np.isclose(np.real(custom_circuit.predict(high_f)), initial_guess[0]) # check complex frequencies raise TypeError @@ -183,7 +185,7 @@ def test_CustomCircuit(): custom_circuit.fit(f, Z) # incorrect circuit element in circuit - with pytest.raises(ValueError): + with pytest.raises(KeyError): custom_circuit = CustomCircuit('R0-NotAnElement', initial_guess=[1, 2]) # test single element circuit diff --git a/impedance/tests/test_fitting.py b/impedance/tests/test_fitting.py index 4bf9ca9..1152d13 100644 --- a/impedance/tests/test_fitting.py +++ b/impedance/tests/test_fitting.py @@ -1,34 +1,40 @@ from impedance.preprocessing import ignoreBelowX -from impedance.models.circuits.fitting import buildCircuit, \ - circuit_fit, rmse, extract_circuit_elements, \ - set_default_bounds -from impedance.tests.test_preprocessing import frequencies \ - as example_frequencies +from impedance.models.circuits.fitting import ( + circuit_fit, + rmse, + extract_circuit_elements, + set_default_bounds, + objective_function, + CircuitGraph, + BoundsCheck, +) +from impedance.tests.test_preprocessing import frequencies as example_frequencies from impedance.tests.test_preprocessing import Z_correct import numpy as np +import pytest def test_set_default_bounds(): # Test example circuit from "Getting Started" page - circuit = 'R0-p(R1,C1)-p(R2-Wo1,C2)' + circuit = "R0-p(R1,C1)-p(R2-Wo1,C2)" # Test with no constants - default_bounds = (np.zeros(7), np.inf*np.ones(7)) + default_bounds = (np.zeros(7), np.inf * np.ones(7)) bounds_from_func = set_default_bounds(circuit) assert np.allclose(default_bounds, bounds_from_func) # Test with constants - constants = {'R0': 1} - default_bounds = (np.zeros(6), np.inf*np.ones(6)) + constants = {"R0": 1} + default_bounds = (np.zeros(6), np.inf * np.ones(6)) bounds_from_func = set_default_bounds(circuit, constants=constants) assert np.allclose(default_bounds, bounds_from_func) # Test with CPEs - circuit = 'R0-p(R1,CPE1)-p(R2-CPE2)' - default_bounds = (np.zeros(7), np.inf*np.ones(7)) + circuit = "R0-p(R1,CPE1)-p(R2-CPE2)" + default_bounds = (np.zeros(7), np.inf * np.ones(7)) default_bounds[1][3] = 1 default_bounds[1][6] = 1 bounds_from_func = set_default_bounds(circuit) @@ -37,9 +43,8 @@ def test_set_default_bounds(): def test_circuit_fit(): - # Test trivial model (10 Ohm resistor) - circuit = 'R0' + circuit = "R0" initial_guess = [10] results_simple = [10] @@ -47,176 +52,289 @@ def test_circuit_fit(): frequencies = np.array([10, 100, 1000]) Z_data = np.array([10, 10, 10]) # impedance is real - assert np.allclose(circuit_fit(frequencies, Z_data, circuit, - initial_guess, constants={}, - global_opt=True)[0], - results_simple, rtol=1e-1) + assert np.allclose( + circuit_fit( + frequencies, Z_data, circuit, initial_guess, constants={}, global_opt=True + )[0], + results_simple, + rtol=1e-1, + ) # check that list inputs work frequency_list = [10, 100, 1000] Z_data_list = [10, 10, 10] - assert np.allclose(circuit_fit(frequency_list, Z_data_list, circuit, - initial_guess, constants={}, - global_opt=True)[0], - results_simple, rtol=1e-1) + assert np.allclose( + circuit_fit( + frequency_list, + Z_data_list, + circuit, + initial_guess, + constants={}, + global_opt=True, + )[0], + results_simple, + rtol=1e-1, + ) # Test example circuit from "Getting Started" page - circuit = 'R0-p(R1,C1)-p(R2-Wo1,C2)' - initial_guess = [.01, .01, 100, .01, .05, 100, 1] - bounds = [(0, 0, 0, 0, 0, 0, 0), - (10, 1, 1e3, 1, 1, 1e4, 100)] + circuit = "R0-p(R1,C1)-p(R2-Wo1,C2)" + initial_guess = [0.01, 0.01, 1, 0.001, 0.1, 100, 0.1] + bounds = [(0, 0, 0, 0, 0, 0, 0), (10, 1, 1e3, 1, 1, 1e4, 100)] # results change slightly using predefined bounds - results_local = np.array([1.65e-2, 8.68e-3, 3.32, 5.39e-3, - 6.31e-2, 2.33e2, 2.20e-1]) + results_local = np.array( + [1.65e-2, 8.78e-3, 3.41, 5.45e-3, 1.43e-1, 1.30e3, 2.23e-1] + ) results_local_bounds = results_local.copy() - results_local_bounds[5] = 2.38e2 - results_local_weighted = np.array([1.64e-2, 9.06e-3, 3.06, - 5.29e-3, 1.45e-1, 1.32e3, 2.02e-1]) - - results_global = np.array([1.65e-2, 5.34e-3, 0.22, 9.15e-3, - 1.31e-1, 1.10e3, 2.78]) + results_local_weighted = np.array( + [1.64e-2, 9.06e-3, 3.06, 5.29e-3, 1.45e-1, 1.32e3, 2.02e-1] + ) + + results_minimize = np.array( + [1.65e-2, 8.67e-3, 3.32, 5.39e-3, 6.37e-2, 2.38e2, 2.20e-1] + ) + results_powell = np.array( + [1.65e-2, 8.69e-3, 3.31, 5.39e-3, 6.45e-2, 2.45e2, 2.19e-1] + ) + results_global = np.array( + [1.65e-2, 8.68e-03, 3.32, 5.39e-3, 6.37e-2, 2.38e2, 2.20e-1] + ) # Filter - example_frequencies_filtered, \ - Z_correct_filtered = ignoreBelowX(example_frequencies, Z_correct) - + example_frequencies_filtered, Z_correct_filtered = ignoreBelowX( + example_frequencies, Z_correct + ) + + graph = CircuitGraph(circuit) + sum_sq_error = objective_function( + graph, example_frequencies_filtered, Z_correct_filtered + ) # Test local fitting - assert np.allclose(circuit_fit(example_frequencies_filtered, - Z_correct_filtered, circuit, - initial_guess, constants={})[0], - results_local, rtol=1e-2) + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + constants={}, + )[0] + assert np.allclose(popt, results_local, rtol=1e-2) + assert sum_sq_error(popt) < 2e-5 # Test local fitting with predefined bounds - assert np.allclose(circuit_fit(example_frequencies_filtered, - Z_correct_filtered, circuit, - initial_guess, bounds=bounds, - constants={})[0], - results_local_bounds, rtol=1e-2) + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + bounds=bounds, + constants={}, + )[0] + assert np.allclose(popt, results_local_bounds, rtol=1e-2) + assert sum_sq_error(popt) < 2e-5 # Test local fitting with predefined weights # Use abs(Z), stacked in order of (Re, Im) components - sigma = np.hstack((np.abs(Z_correct_filtered), - np.abs(Z_correct_filtered))) - assert np.allclose(circuit_fit(example_frequencies_filtered, - Z_correct_filtered, circuit, - initial_guess, sigma=sigma, - constants={})[0], - results_local_weighted, rtol=1e-2) + sigma = np.hstack((np.abs(Z_correct_filtered), np.abs(Z_correct_filtered))) + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + sigma=sigma, + constants={}, + )[0] + assert np.allclose(popt, results_local_weighted, rtol=1e-2) + assert sum_sq_error(popt) < 2e-5 # Test if using weight_by_modulus=True produces the same results - assert np.allclose(circuit_fit(example_frequencies_filtered, - Z_correct_filtered, circuit, - initial_guess, weight_by_modulus=True, - constants={})[0], - results_local_weighted, rtol=1e-2) + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + weight_by_modulus=True, + constants={}, + )[0] + assert np.allclose(popt, results_local_weighted, rtol=1e-2) + assert sum_sq_error(popt) < 2e-5 + + # Test if using method "minimize" produces the same results + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + opt_method="minimize", + constants={}, + )[0] + assert np.allclose(popt, results_minimize, rtol=1e-2) + assert sum_sq_error(popt) < 2e-5 + + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + opt_method="minimize", + constants={}, + method="Powell", + )[0] + assert np.allclose(popt, results_powell, rtol=1e-2) + assert sum_sq_error(popt) < 2e-5 + + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + opt_method="minimize", + constants={}, + method="TNC", + options={"scale": initial_guess, "maxfun": 10_000}, + )[0] + assert sum_sq_error(popt) < 1e-4 # Test global fitting on multiple seeds # All seeds should converge to the same parameter values # seed = 0 (default) - assert np.allclose(circuit_fit(example_frequencies_filtered, - Z_correct_filtered, circuit, - initial_guess, constants={}, - global_opt=True)[0], - results_global, rtol=1e-1) + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + constants={}, + global_opt=True, + )[0] + assert np.allclose(popt, results_global, rtol=1e-1) + assert sum_sq_error(popt) < 2e-5 + + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + constants={}, + global_opt=True, + minimizer_kwargs={ + "method": "TNC", + "options": { + "scale": initial_guess, + "maxfun": 10_000, + } + }, + )[0] + assert sum_sq_error(popt) < 5e-5 # seed = 0, with predefined bounds - assert np.allclose(circuit_fit(example_frequencies_filtered, - Z_correct_filtered, circuit, - initial_guess, constants={}, - global_opt=True, bounds=bounds, - seed=0)[0], - results_global, rtol=1e-1) + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + constants={}, + global_opt=True, + bounds=bounds, + seed=0, + )[0] + assert np.allclose(popt, results_global, rtol=1e-1) + assert sum_sq_error(popt) < 2e-5 # seed = 1 - assert np.allclose(circuit_fit(example_frequencies_filtered, - Z_correct_filtered, circuit, - initial_guess, constants={}, - global_opt=True, seed=1)[0], - results_global, rtol=1e-1) + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + constants={}, + global_opt=True, + seed=1, + )[0] + assert np.allclose(popt, results_global, rtol=1e-1) + assert sum_sq_error(popt) < 2e-5 # seed = 42 - assert np.allclose(circuit_fit(example_frequencies_filtered, - Z_correct_filtered, circuit, - initial_guess, constants={}, - global_opt=True, seed=42)[0], - results_global, rtol=1e-1) - - -def test_buildCircuit(): - + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + constants={}, + global_opt=True, + seed=42, + )[0] + assert np.allclose(popt, results_global, rtol=1e-1) + assert sum_sq_error(popt) < 2e-5 + + # Test dual annealing + popt = circuit_fit( + example_frequencies_filtered, + Z_correct_filtered, + circuit, + initial_guess, + constants={}, + global_opt=True, + opt_method="dual_annealing", + )[0] + assert np.allclose(popt, results_minimize, rtol=1e-1) + assert sum_sq_error(popt) < 5e-5 + + +def test_CircuitGraph(): # Test simple Randles circuit with CPE - circuit = 'R0-p(R1-Wo1,CPE1)' - params = [.1, .01, 1, 1000, 15, .9] + circuit = "R0-p(R1-Wo1,CPE1)" + params = [0.1, 0.01, 1, 1000, 15, 0.9] frequencies = [1000.0, 5.0, 0.01] - assert buildCircuit(circuit, frequencies, *params, - constants={})[0].replace(' ', '') == \ - 's([R([0.1],[1000.0,5.0,0.01]),' + \ - 'p([s([R([0.01],[1000.0,5.0,0.01]),' + \ - 'Wo([1.0,1000.0],[1000.0,5.0,0.01])]),' + \ - 'CPE([15.0,0.9],[1000.0,5.0,0.01])])])' + cg = CircuitGraph(circuit) + assert len(cg.execution_order) == 7 + assert len(cg.compute(frequencies, *params)) == len(frequencies) + assert cg.calculate_circuit_length() == 6 + + import matplotlib.pyplot as plt + + f, ax = plt.subplots() + cg.visualize_graph(ax=ax) + plt.close(f) # Test multiple parallel elements - circuit = 'R0-p(C1,R1,R2)' - params = [.1, .01, .2, .3] - frequencies = [1000.0, 5.0, 0.01] + circuit = "R0-p(C1,R1,R2)" - assert buildCircuit(circuit, frequencies, *params, - constants={})[0].replace(' ', '') == \ - 's([R([0.1],[1000.0,5.0,0.01]),' + \ - 'p([C([0.01],[1000.0,5.0,0.01]),' + \ - 'R([0.2],[1000.0,5.0,0.01]),' + \ - 'R([0.3],[1000.0,5.0,0.01])])])' + assert len(CircuitGraph(circuit).execution_order) == 6 # Test nested parallel groups - circuit = 'R0-p(p(R1, C1)-R2, C2)' - params = [1, 2, 3, 4, 5] - frequencies = [1000.0, 5.0, 0.01] + circuit = "R0-p(p(R1, C1)-R2, C2)" - assert buildCircuit(circuit, frequencies, *params, - constants={})[0].replace(' ', '') == \ - 's([R([1],[1000.0,5.0,0.01]),' + \ - 'p([s([p([R([2],[1000.0,5.0,0.01]),' + \ - 'C([3],[1000.0,5.0,0.01])]),' + \ - 'R([4],[1000.0,5.0,0.01])]),' + \ - 'C([5],[1000.0,5.0,0.01])])])' + assert len(CircuitGraph(circuit).execution_order) == 9 # Test parallel elements at beginning and end - circuit = 'p(C1,R1)-p(C2,R2)' - params = [.1, .01, .2, .3] - frequencies = [1000.0, 5.0, 0.01] + circuit = "p(C1,R1)-p(C2,R2)" - assert buildCircuit(circuit, frequencies, *params, - constants={})[0].replace(' ', '') == \ - 's([p([C([0.1],[1000.0,5.0,0.01]),' + \ - 'R([0.01],[1000.0,5.0,0.01])]),' + \ - 'p([C([0.2],[1000.0,5.0,0.01]),' + \ - 'R([0.3],[1000.0,5.0,0.01])])])' + assert len(CircuitGraph(circuit).execution_order) == 7 # Test single element circuit - circuit = 'R1' - params = [100] - frequencies = [1000.0, 5.0, 0.01] + circuit = "R1" - assert buildCircuit(circuit, frequencies, *params, - constants={})[0].replace(' ', '') == \ - 'R([100],[1000.0,5.0,0.01])' + assert len(CircuitGraph(circuit).execution_order) == 1 def test_RMSE(): - a = np.array([2 + 4*1j, 3 + 2*1j]) - b = np.array([2 + 4*1j, 3 + 2*1j]) + a = np.array([2 + 4 * 1j, 3 + 2 * 1j]) + b = np.array([2 + 4 * 1j, 3 + 2 * 1j]) assert rmse(a, b) == 0.0 - c = np.array([2 + 4*1j, 1 + 4*1j]) - d = np.array([4 + 2*1j, 3 + 2*1j]) - assert np.isclose(rmse(c, d), 2*np.sqrt(2)) + c = np.array([2 + 4 * 1j, 1 + 4 * 1j]) + d = np.array([4 + 2 * 1j, 3 + 2 * 1j]) + assert np.isclose(rmse(c, d), 2 * np.sqrt(2)) def test_element_extraction(): - circuit = 'R0-p(RR0,C1)-p(R1,C2032478)-W1' + circuit = "R0-p(RR0,C1)-p(R1,C2032478)-W1" extracted_elements = extract_circuit_elements(circuit) - assert extracted_elements == ['R0', 'RR0', 'C1', 'R1', 'C2032478', 'W1'] + assert extracted_elements == ["R0", "RR0", "C1", "R1", "C2032478", "W1"] + + +def test_bounds(): + with pytest.raises(ValueError): + BoundsCheck(np.array([1, 2, 8]), np.array([2, 3, 4])) diff --git a/requirements.txt b/requirements.txt index ea8a8d3..a83f7fc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,8 @@ altair>=3.0 coveralls==3.2.0 matplotlib>=3.5 +networkx>=2.6.3 +numdifftools>=0.9.41 numpy>=1.14 pytest>=4.6 pytest-cov diff --git a/setup.py b/setup.py index 610ccfb..68cc777 100644 --- a/setup.py +++ b/setup.py @@ -16,7 +16,8 @@ packages=setuptools.find_packages(), python_requires="~=3.7", install_requires=['altair>=3.0', 'matplotlib>=3.5', - 'numpy>=1.14', 'scipy>=1.0'], + 'numdifftools>=0.9.41', 'networkx>=2.6.3', + 'numpy>=1.14', 'scipy>=1.0',], classifiers=( "Programming Language :: Python :: 3", "License :: OSI Approved :: MIT License", From a323a8a4b7cd788486b1383bb4b3796783cee780 Mon Sep 17 00:00:00 2001 From: Jake Anderson Date: Thu, 12 Sep 2024 09:55:19 -0600 Subject: [PATCH 2/2] update formatting to pass linting --- impedance/models/circuits/circuits.py | 36 +++++++++++++++++++++------ impedance/models/circuits/elements.py | 32 ++++++++++++------------ impedance/models/circuits/fitting.py | 32 ++++++++++++++++++------ impedance/tests/test_fitting.py | 13 +++++++--- setup.py | 13 +++++++--- 5 files changed, 88 insertions(+), 38 deletions(-) diff --git a/impedance/models/circuits/circuits.py b/impedance/models/circuits/circuits.py index a801378..534415f 100644 --- a/impedance/models/circuits/circuits.py +++ b/impedance/models/circuits/circuits.py @@ -20,7 +20,9 @@ class BaseCircuit: """Base class for equivalent circuit models""" - def __init__(self, circuit="", initial_guess=None, constants=None, name=None): + def __init__( + self, circuit="", initial_guess=None, constants=None, name=None + ): """Base constructor for any equivalent circuit model Parameters @@ -44,7 +46,9 @@ def __init__(self, circuit="", initial_guess=None, constants=None, name=None): # if supplied, check that initial_guess is valid and store initial_guess = [x for x in initial_guess if x is not None] if not np.all(np.isfinite(initial_guess)): - raise TypeError(f"initial_guess must be numeric or None ({initial_guess})") + raise TypeError( + f"initial_guess must be numeric or None ({initial_guess})" + ) # initalize class attributes self.initial_guess = initial_guess @@ -86,7 +90,12 @@ def __eq__(self, other): raise TypeError("Comparing object is not of the same type.") def fit( - self, frequencies, impedance, bounds=None, weight_by_modulus=False, **kwargs + self, + frequencies, + impedance, + bounds=None, + weight_by_modulus=False, + **kwargs, ): """Fit the circuit model @@ -184,10 +193,16 @@ def get_param_names(self): full_names, all_units = [], [] for name in names: - p_names, p_units = get_element_parameter_names_and_units_from_name(name) + p_names, p_units = get_element_parameter_names_and_units_from_name( + name + ) full_names.extend([n for n in p_names if n not in self.constants]) all_units.extend( - [u for u, n in zip(p_units, p_names) if n not in self.constants] + [ + u + for u, n in zip(p_units, p_names) + if n not in self.constants + ] ) return full_names, all_units @@ -283,13 +298,18 @@ def plot(self, ax=None, f_data=None, Z_data=None, kind="altair", **kwargs): if Z_data is not None: if f_data is None: raise ValueError( - "f_data must be specified if" + " Z_data for a Bode plot" + "f_data must be specified if" + + " Z_data for a Bode plot" ) - ax = plot_bode(f_data, Z_data, ls="", marker="s", axes=ax, **kwargs) + ax = plot_bode( + f_data, Z_data, ls="", marker="s", axes=ax, **kwargs + ) if self._is_fit(): Z_fit = self.predict(f_pred) - ax = plot_bode(f_pred, Z_fit, ls="-", marker="", axes=ax, **kwargs) + ax = plot_bode( + f_pred, Z_fit, ls="-", marker="", axes=ax, **kwargs + ) return ax elif kind == "altair": plot_dict = {} diff --git a/impedance/models/circuits/elements.py b/impedance/models/circuits/elements.py index 0a95b03..ae98351 100644 --- a/impedance/models/circuits/elements.py +++ b/impedance/models/circuits/elements.py @@ -1,12 +1,10 @@ import numpy as np -class ElementError(Exception): - ... +class ElementError(Exception): ... -class OverwriteError(ElementError): - ... +class OverwriteError(ElementError): ... def element(num_params, units, overwrite=False): @@ -35,13 +33,14 @@ def wrapper(p, f): global circuit_elements if func.__name__ in ["s", "p"]: - raise ElementError("cannot redefine elements 's' (series)" + - "or 'p' (parallel)") + raise ElementError( + "cannot redefine elements 's' (series)" + "or 'p' (parallel)" + ) elif func.__name__ in circuit_elements and not overwrite: raise OverwriteError( - f"element {func.__name__} already exists. " + - "If you want to overwrite the existing element," + - "use `overwrite=True`." + f"element {func.__name__} already exists. " + + "If you want to overwrite the existing element," + + "use `overwrite=True`." ) else: circuit_elements[func.__name__] = wrapper @@ -313,9 +312,9 @@ def K(p, f): return Z -@element(num_params=3, units=['Ohm', 'sec', '']) +@element(num_params=3, units=["Ohm", "sec", ""]) def Zarc(p, f): - """ An RQ element rewritten with resistance and + """An RQ element rewritten with resistance and and time constant as paramenters. Equivalent to a Cole-Cole relaxation in dielectrics. @@ -326,9 +325,9 @@ def Zarc(p, f): Z = \\frac{R}{1 + (j \\omega \\tau_k)^\\gamma } """ - omega = 2*np.pi*np.array(f) + omega = 2 * np.pi * np.array(f) R, tau_k, gamma = p[0], p[1], p[2] - Z = R/(1 + ((1j*omega*tau_k)**gamma)) + Z = R / (1 + ((1j * omega * tau_k) ** gamma)) return Z @@ -413,13 +412,14 @@ def typeChecker(p, f, name, length): if len(p) != length: raise TypeError(f"In {name} length of parameters is not {length} ") + def get_element_parameter_names_and_units_from_name(name): elem = get_element_from_name(name) n_params = circuit_elements[elem].num_params return [ - format_parameter_name(name, j, n_params) - for j in range(n_params) + format_parameter_name(name, j, n_params) for j in range(n_params) ], circuit_elements[elem].units + def format_parameter_name(name, j, n_params): - return f"{name}_{j}" if n_params > 1 else f"{name}" \ No newline at end of file + return f"{name}_{j}" if n_params > 1 else f"{name}" diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index 7bc1e25..473f359 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -4,9 +4,19 @@ import numpy as np import numdifftools as nd from scipy.linalg import inv, norm -from scipy.optimize import curve_fit, dual_annealing, minimize, basinhopping, Bounds - -from .elements import circuit_elements, get_element_from_name, format_parameter_name +from scipy.optimize import ( + curve_fit, + dual_annealing, + minimize, + basinhopping, + Bounds, +) + +from .elements import ( + circuit_elements, + get_element_from_name, + format_parameter_name, +) ints = "0123456789" _LBFGSB_OPTIONS = { @@ -78,7 +88,9 @@ def cost_function(x): function Returns a function """ - return sum_square_error(graph(f, *x), np.stack([Z.real, Z.imag]), weights) + return sum_square_error( + graph(f, *x), np.stack([Z.real, Z.imag]), weights + ) return cost_function @@ -218,7 +230,9 @@ def circuit_fit( bounds = set_default_bounds(circuit, constants=constants) cg = CircuitGraph(circuit, constants) - sumsq_errors = objective_function(cg, f, Z, weight_by_modulus=weight_by_modulus) + sumsq_errors = objective_function( + cg, f, Z, weight_by_modulus=weight_by_modulus + ) if not global_opt: if opt_method == "minimize": if "options" not in kwargs: @@ -291,7 +305,9 @@ def circuit_fit( ) popt = results.x - pcov = results.get("lowest_optimization_result", {}).get("hess_inv", None) + pcov = results.get("lowest_optimization_result", {}).get( + "hess_inv", None + ) try: pcov = pcov.todense() except AttributeError: @@ -442,7 +458,9 @@ def compute(self, f, *parameters): pindex = 0 for node in self.execution_order: Zfunc = self.graph.nodes[node]["Z"] - plist = [node_results[pred] for pred in self.graph.predecessors(node)] + plist = [ + node_results[pred] for pred in self.graph.predecessors(node) + ] if len(plist) < 1: n_params = Zfunc.num_params for j in range(n_params): diff --git a/impedance/tests/test_fitting.py b/impedance/tests/test_fitting.py index 1152d13..a47d4a2 100644 --- a/impedance/tests/test_fitting.py +++ b/impedance/tests/test_fitting.py @@ -8,7 +8,9 @@ CircuitGraph, BoundsCheck, ) -from impedance.tests.test_preprocessing import frequencies as example_frequencies +from impedance.tests.test_preprocessing import ( + frequencies as example_frequencies, +) from impedance.tests.test_preprocessing import Z_correct import numpy as np @@ -54,7 +56,12 @@ def test_circuit_fit(): assert np.allclose( circuit_fit( - frequencies, Z_data, circuit, initial_guess, constants={}, global_opt=True + frequencies, + Z_data, + circuit, + initial_guess, + constants={}, + global_opt=True, )[0], results_simple, rtol=1e-1, @@ -221,7 +228,7 @@ def test_circuit_fit(): "options": { "scale": initial_guess, "maxfun": 10_000, - } + }, }, )[0] assert sum_sq_error(popt) < 5e-5 diff --git a/setup.py b/setup.py index 68cc777..ba4ca8c 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ from impedance import __version__ import setuptools -with open("README.md", "r", encoding='utf8') as fh: +with open("README.md", "r", encoding="utf8") as fh: long_description = fh.read() setuptools.setup( @@ -15,9 +15,14 @@ url="https://impedancepy.readthedocs.io/en/latest/", packages=setuptools.find_packages(), python_requires="~=3.7", - install_requires=['altair>=3.0', 'matplotlib>=3.5', - 'numdifftools>=0.9.41', 'networkx>=2.6.3', - 'numpy>=1.14', 'scipy>=1.0',], + install_requires=[ + "altair>=3.0", + "matplotlib>=3.5", + "numdifftools>=0.9.41", + "networkx>=2.6.3", + "numpy>=1.14", + "scipy>=1.0", + ], classifiers=( "Programming Language :: Python :: 3", "License :: OSI Approved :: MIT License",