From e2516a1b67eef881eee6a5481e6935b0b55ce6a4 Mon Sep 17 00:00:00 2001 From: petermattia Date: Tue, 9 Mar 2021 10:39:19 -0800 Subject: [PATCH 1/7] Add global optimization function - Add scipy.optimize.basinhopping global optimization algorithm - Added test_circuit_fit tests - Added docstrings - Renamed some single-letter variables --- impedance/models/circuits/fitting.py | 120 +++++++++++++++++++++++---- impedance/tests/test_fitting.py | 37 ++++++++- 2 files changed, 137 insertions(+), 20 deletions(-) diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index 96c261ac..265069fb 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -1,6 +1,7 @@ from .elements import circuit_elements, get_element_from_name import numpy as np -from scipy.optimize import curve_fit +from scipy.linalg import inv +from scipy.optimize import curve_fit, basinhopping ints = '0123456789' @@ -22,9 +23,18 @@ def rmse(a, b): def circuit_fit(frequencies, impedances, circuit, initial_guess, constants, - bounds=None, bootstrap=False, **kwargs): + bounds=None, bootstrap=False, global_opt=False, **kwargs): - """ Main function for fitting an equivalent circuit to data + """ Main function for fitting an equivalent circuit to data. + + By default, this function uses `scipy.optimize.curve_fit + ` + to fit the equivalent circuit. This function generally works well for + simple circuits. However, the final results may be sensitive to + the initial conditions for more complex circuits. In these cases, + the `scipy.optimize.basinhopping + `) + global optimization algorithm can be used to attempt a better fit. Parameters ----------------- @@ -49,8 +59,13 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants, parameters of 0 and np.inf, except the CPE alpha which has an upper bound of 1 + global_opt : bool, optional + If global optimization should be used (uses the shgo algorithm). + Defaults to False + kwargs : - Keyword arguments passed to scipy.optimize.curve_fit + Keyword arguments passed to scipy.optimize.curve_fit or + scipy.optimize.basinhopping Returns ------------ @@ -66,7 +81,6 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants, Currently, an error of -1 is returned. """ - f = frequencies Z = impedances # extract the elements from the circuit @@ -74,29 +88,56 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants, # set upper and lower bounds on a per-element basis if bounds is None: - lb, ub = [], [] + 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 + '_{}'.format(i) in constants: continue if raw_element in ['CPE', 'La'] and i == 1: - ub.append(1) + upper_bounds.append(1) else: - ub.append(np.inf) - lb.append(0) - bounds = ((lb), (ub)) + upper_bounds.append(np.inf) + lower_bounds.append(0) + bounds = ((lower_bounds), (upper_bounds)) - if 'maxfev' not in kwargs: - kwargs['maxfev'] = 100000 - if 'ftol' not in kwargs: - kwargs['ftol'] = 1e-13 + if not global_opt: + if 'maxfev' not in kwargs: + kwargs['maxfev'] = 100000 + if 'ftol' not in kwargs: + kwargs['ftol'] = 1e-13 - popt, pcov = curve_fit(wrapCircuit(circuit, constants), f, - np.hstack([Z.real, Z.imag]), p0=initial_guess, - bounds=bounds, **kwargs) + popt, pcov = curve_fit(wrapCircuit(circuit, constants), frequencies, + np.hstack([Z.real, Z.imag]), + p0=initial_guess, bounds=bounds, **kwargs) - perror = np.sqrt(np.diag(pcov)) + perror = np.sqrt(np.diag(pcov)) + + else: + 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)(frequencies, *x), + np.hstack([Z.real, Z.imag])) + + results = basinhopping(opt_function, x0=initial_guess) + popt = results.x + + # jacobian -> covariance + # https://stats.stackexchange.com/q/231868 + jac = results.lowest_optimization_result["jac"][np.newaxis] + perror = inv(np.dot(jac.T, jac)) * opt_function(popt) ** 2 return popt, perror @@ -247,6 +288,19 @@ def count_parens(string): def extract_circuit_elements(circuit): + """ Extracts circuit elements from a circuit string. + + Parameters + ---------- + circuit : str + Circuit string. + + Returns + ------- + extracted_elements : list + list of extracted elements. + + """ p_string = [x for x in circuit if x not in 'p(),-'] extracted_elements = [] current_element = [] @@ -267,6 +321,19 @@ def extract_circuit_elements(circuit): def calculateCircuitLength(circuit): + """ Calculates the number of elements in the circuit. + + Parameters + ---------- + circuit : str + Circuit string. + + Returns + ------- + length : int + Length of circuit. + + """ length = 0 if circuit: extracted_elements = extract_circuit_elements(circuit) @@ -278,6 +345,23 @@ def calculateCircuitLength(circuit): def check_and_eval(element): + """ Checks if an element is valid, then evaluates it. + + Parameters + ---------- + element : str + Circuit element. + + Raises + ------ + ValueError + Raised if an element is not in the list of allowed elements. + + Returns + ------- + Evaluated element. + + """ allowed_elements = circuit_elements.keys() if element not in allowed_elements: raise ValueError(f'{element} not in ' + diff --git a/impedance/tests/test_fitting.py b/impedance/tests/test_fitting.py index f7dc7464..60a13a80 100644 --- a/impedance/tests/test_fitting.py +++ b/impedance/tests/test_fitting.py @@ -1,5 +1,10 @@ -from impedance.models.circuits.fitting import buildCircuit, rmse, \ - extract_circuit_elements +from impedance.preprocessing import ignoreBelowX +from impedance.models.circuits.fitting import buildCircuit, \ + circuit_fit, rmse, extract_circuit_elements +from impedance.tests.test_preprocessing import frequencies \ + as example_frequencies +from impedance.tests.test_preprocessing import Z_correct + import numpy as np # def test_residuals(): @@ -12,6 +17,34 @@ # +def test_circuit_fit(): + + # 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] + results_local = np.array([1.65e-2, 8.68e-3, 3.32e0, 5.39e-3, + 6.31e-2, 2.33e2, 2.20e-1]) + results_global = np.array([1.65e-2, 8.78e-3, 3.41, 5.45e-3, + 1.32e-1, 1.10e3, 2.23e-1]) + + # Filter + example_frequencies_filtered, \ + Z_correct_filtered = ignoreBelowX(example_frequencies, Z_correct) + + # Test local fitting + assert np.isclose(circuit_fit(example_frequencies_filtered, + Z_correct_filtered, circuit, + initial_guess, constants={})[0], + results_local, rtol=1e-2).all() + + # Test global fitting + assert np.isclose(circuit_fit(example_frequencies_filtered, + Z_correct_filtered, circuit, + initial_guess, constants={}, + global_opt=True)[0], + results_global, rtol=1e-1).all() + + def test_buildCircuit(): # Test simple Randles circuit with CPE From af661b562c96a00d94612ab232f5bb3a79fc6164 Mon Sep 17 00:00:00 2001 From: petermattia Date: Tue, 9 Mar 2021 11:03:01 -0800 Subject: [PATCH 2/7] Fix typo --- impedance/models/circuits/fitting.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index 265069fb..8785e781 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -60,8 +60,8 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants, which has an upper bound of 1 global_opt : bool, optional - If global optimization should be used (uses the shgo algorithm). - Defaults to False + If global optimization should be used (uses the basinhopping + algorithm). Defaults to False kwargs : Keyword arguments passed to scipy.optimize.curve_fit or From 414c33756ef62ad087b83fa7af59119d8995ca0d Mon Sep 17 00:00:00 2001 From: Peter Attia Date: Tue, 9 Mar 2021 13:47:54 -0800 Subject: [PATCH 3/7] Update fitting.py flake8 compliance --- impedance/models/circuits/fitting.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index 8785e781..7a0dee32 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -27,12 +27,12 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants, """ Main function for fitting an equivalent circuit to data. - By default, this function uses `scipy.optimize.curve_fit + By default, this function uses `scipy.optimize.curve_fit ` to fit the equivalent circuit. This function generally works well for simple circuits. However, the final results may be sensitive to the initial conditions for more complex circuits. In these cases, - the `scipy.optimize.basinhopping + the `scipy.optimize.basinhopping `) global optimization algorithm can be used to attempt a better fit. From a044e2248be2530ccd6e7c2e4cd2782c1d8afb21 Mon Sep 17 00:00:00 2001 From: petermattia Date: Tue, 9 Mar 2021 14:21:26 -0800 Subject: [PATCH 4/7] Add seed and catch LinAlgErrors in computing perror --- impedance/models/circuits/fitting.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index 7a0dee32..947a1ba8 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -23,7 +23,8 @@ def rmse(a, b): def circuit_fit(frequencies, impedances, circuit, initial_guess, constants, - bounds=None, bootstrap=False, global_opt=False, **kwargs): + bounds=None, bootstrap=False, global_opt=False, seed=0, + **kwargs): """ Main function for fitting an equivalent circuit to data. @@ -63,6 +64,10 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants, If global optimization should be used (uses the basinhopping algorithm). Defaults to False + seed : int, optional + Random seed, only used for the basinhopping algorithm. + Defaults to 0 + kwargs : Keyword arguments passed to scipy.optimize.curve_fit or scipy.optimize.basinhopping @@ -137,7 +142,11 @@ def opt_function(x): # jacobian -> covariance # https://stats.stackexchange.com/q/231868 jac = results.lowest_optimization_result["jac"][np.newaxis] - perror = inv(np.dot(jac.T, jac)) * opt_function(popt) ** 2 + try: + perror = inv(np.dot(jac.T, jac)) * opt_function(popt) ** 2 + except np.LinAlg.LinAlgError: + print('Failed to compute perror') + perror = None return popt, perror From 70654d0014a8af36e595a490a68eb499fac5d55d Mon Sep 17 00:00:00 2001 From: petermattia Date: Tue, 9 Mar 2021 14:42:41 -0800 Subject: [PATCH 5/7] Update fitting.py typo in linalgerror --- impedance/models/circuits/fitting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index 947a1ba8..f5a44ba5 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -144,7 +144,7 @@ def opt_function(x): jac = results.lowest_optimization_result["jac"][np.newaxis] try: perror = inv(np.dot(jac.T, jac)) * opt_function(popt) ** 2 - except np.LinAlg.LinAlgError: + except np.linalg.LinAlgError: print('Failed to compute perror') perror = None From e9f07a05b95e73d703ad50ed9dfb5ca173612e89 Mon Sep 17 00:00:00 2001 From: petermattia Date: Wed, 10 Mar 2021 23:21:54 -0800 Subject: [PATCH 6/7] Replace print with warnings.warn --- impedance/models/circuits/fitting.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index f5a44ba5..52574dcb 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -1,8 +1,11 @@ -from .elements import circuit_elements, get_element_from_name +import warnings + import numpy as np from scipy.linalg import inv from scipy.optimize import curve_fit, basinhopping +from .elements import circuit_elements, get_element_from_name + ints = '0123456789' @@ -145,7 +148,7 @@ def opt_function(x): try: perror = inv(np.dot(jac.T, jac)) * opt_function(popt) ** 2 except np.linalg.LinAlgError: - print('Failed to compute perror') + warnings.warn("Failed to compute perror") perror = None return popt, perror From 6f48265d611b7e1be116713505d6de5c9a97f3c1 Mon Sep 17 00:00:00 2001 From: petermattia Date: Wed, 10 Mar 2021 23:30:08 -0800 Subject: [PATCH 7/7] (Actually) add seed to basinhopping --- impedance/models/circuits/fitting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/impedance/models/circuits/fitting.py b/impedance/models/circuits/fitting.py index 52574dcb..cfaa2720 100644 --- a/impedance/models/circuits/fitting.py +++ b/impedance/models/circuits/fitting.py @@ -139,7 +139,7 @@ def opt_function(x): return rmse(wrapCircuit(circuit, constants)(frequencies, *x), np.hstack([Z.real, Z.imag])) - results = basinhopping(opt_function, x0=initial_guess) + results = basinhopping(opt_function, x0=initial_guess, seed=seed) popt = results.x # jacobian -> covariance