From 9de9069b8d7d33af6b5fd81b0cc34db3935009e1 Mon Sep 17 00:00:00 2001 From: Yuefan Ji Date: Sun, 20 Oct 2024 13:49:55 -0700 Subject: [PATCH 1/6] added supports for nonlinear RC with constant phase element --- nleis/fitting.py | 11 +- nleis/nleis_elements_pair.py | 179 +++++++++++++++--- nleis/nleis_tests/test_nleis_elements_pair.py | 38 ++++ 3 files changed, 203 insertions(+), 25 deletions(-) diff --git a/nleis/fitting.py b/nleis/fitting.py index 7393ab1..4a12c2d 100644 --- a/nleis/fitting.py +++ b/nleis/fitting.py @@ -110,7 +110,6 @@ def seq_fit_param(input_dic, target_arr, output_arr): def set_default_bounds(circuit, constants={}): - """ Set default bounds for optimization. @@ -170,10 +169,12 @@ def set_default_bounds(circuit, constants={}): elif raw_element in ['RCn'] and i == 2: upper_bounds.append(0.5) lower_bounds.append(-0.5) - elif raw_element in ['TDSn', 'TDPn', 'TDCn'] and (i == 5): + elif raw_element in ['TDSn', 'TDPn', 'TDCn', 'RCSQn', + 'RCDQn'] and (i == 5): upper_bounds.append(np.inf) lower_bounds.append(-np.inf) - elif raw_element in ['TDSn', 'TDPn', 'TDCn'] and i == 6: + elif raw_element in ['TDSn', 'TDPn', 'TDCn', 'RCSQn', + 'RCDQn'] and i == 6: upper_bounds.append(0.5) lower_bounds.append(-0.5) elif raw_element in ['RCDn', 'RCSn'] and (i == 4): @@ -191,6 +192,10 @@ def set_default_bounds(circuit, constants={}): elif raw_element in ['TLMSn', 'TLMDn'] and (i == 8): upper_bounds.append(np.inf) lower_bounds.append(-np.inf) + elif raw_element in ['RCSQ', 'RCSQn', + 'RCDQ', 'RCDQn'] and (i == 2): + upper_bounds.append(1) + lower_bounds.append(0) else: upper_bounds.append(np.inf) lower_bounds.append(0) diff --git a/nleis/nleis_elements_pair.py b/nleis/nleis_elements_pair.py index 66b98b8..864c4c5 100644 --- a/nleis/nleis_elements_pair.py +++ b/nleis/nleis_elements_pair.py @@ -117,7 +117,7 @@ def RC(p, f): return Rct / (1 + ω_star*1j) -@element(num_params=3, units=['Ohm', 'F', '']) +@element(num_params=3, units=['Ohm', 'F', '-']) def RCn(p, f): ''' @@ -215,7 +215,7 @@ def RCD(p, f): return (Z) -@element(num_params=6, units=['Ohms', 'F', 'Ohms', 's', '1/V', '']) +@element(num_params=6, units=['Ohms', 'F', 'Ohms', 's', '1/V', '-']) def RCDn(p, f): ''' @@ -352,7 +352,7 @@ def RCS(p, f): return (Z) -@element(num_params=6, units=['Ohms', 'F', 'Ohms', 's', '1/V', '']) +@element(num_params=6, units=['Ohms', 'F', 'Ohms', 's', '1/V', '-']) def RCSn(p, f): ''' @@ -485,7 +485,7 @@ def TP(p, f): return Z -@element(num_params=4, units=['Ohms', 'Ohms', 'F', '']) +@element(num_params=4, units=['Ohms', 'Ohms', 'F', '-']) def TPn(p, f): """ @@ -634,7 +634,7 @@ def TDP(p, f): return Z -@element(num_params=7, units=['Ohms', 'Ohms', 'F', 'Ohms', 's', '1/V', '']) +@element(num_params=7, units=['Ohms', 'Ohms', 'F', 'Ohms', 's', '1/V', '-']) def TDPn(p, f): """ @@ -803,7 +803,7 @@ def TDS(p, f): return Z -@element(num_params=7, units=['Ohms', 'Ohms', 'F', 'Ohms', 's', '1/V', '']) +@element(num_params=7, units=['Ohms', 'Ohms', 'F', 'Ohms', 's', '1/V', '-']) def TDSn(p, f): """ @@ -979,7 +979,7 @@ def TDC(p, f): return Z -@element(num_params=7, units=['Ohms', 'Ohms', 'F', 'Ohms', 's', '1/V', '']) +@element(num_params=7, units=['Ohms', 'Ohms', 'F', 'Ohms', 's', '1/V', '-']) def TDCn(p, f): """ @@ -1107,7 +1107,7 @@ def TDCn(p, f): # TLM Model # -@element(num_params=6, units=['Ohm', 'Ohm', 'F', 'Ohm', 'F', '']) +@element(num_params=6, units=['Ohm', 'Ohm', 'F', 'Ohm', 'F', '-']) def TLM(p, f): """ @@ -1157,7 +1157,7 @@ def TLM(p, f): ### -@element(num_params=8, units=['Ohm', 'Ohm', 'F', '', 'Ohm', 'F', '', '']) +@element(num_params=8, units=['Ohm', 'Ohm', 'F', '-', 'Ohm', 'F', '-', '-']) def TLMn(p, f): """ @@ -1262,7 +1262,7 @@ def TLMn(p, f): return (Z) -@element(num_params=6, units=['Ohm', 'Ohm', 'F', 'Ohm', 'F', '']) +@element(num_params=6, units=['Ohm', 'Ohm', 'F', 'Ohm', 'F', '-']) def mTi(p, f): """ @@ -1343,7 +1343,7 @@ def mTi(p, f): return (I1) -@element(num_params=8, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '']) +@element(num_params=8, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '-']) def TLMS(p, f): """ @@ -1407,8 +1407,8 @@ def TLMS(p, f): return (Req) -@element(num_params=11, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '', - '1/V', '', '']) +@element(num_params=11, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '-', + '1/V', '-', '-']) def TLMSn(p, f): """ @@ -1520,7 +1520,7 @@ def TLMSn(p, f): return (Z) -@element(num_params=8, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '']) +@element(num_params=8, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '-']) def mTiS(p, f): """ @@ -1606,8 +1606,8 @@ def mTiS(p, f): return (I1) -@element(num_params=11, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '', - '1/V', '', '']) +@element(num_params=11, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '-', + '1/V', '-', '-']) def mTiSn(p, f): """ @@ -1721,7 +1721,7 @@ def mTiSn(p, f): return (I2) -@element(num_params=8, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '']) +@element(num_params=8, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '-']) def TLMD(p, f): """ @@ -1783,8 +1783,8 @@ def TLMD(p, f): return (Req) -@element(num_params=11, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '', - '1/V', '', '']) +@element(num_params=11, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '-', + '1/V', '-', '-']) def TLMDn(p, f): """ @@ -1894,7 +1894,7 @@ def TLMDn(p, f): return (Z) -@element(num_params=8, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '']) +@element(num_params=8, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '-']) def mTiD(p, f): """ @@ -1980,8 +1980,8 @@ def mTiD(p, f): return (I1) -@element(num_params=11, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '', - '1/V', '', '']) +@element(num_params=11, units=['Ohm', 'Ohm', 'F', 'Ohm', 's', 'Ohm', 'F', '-', + '1/V', '-', '-']) def mTiDn(p, f): """ @@ -2094,6 +2094,141 @@ def mTiDn(p, f): return (I2) +@element(num_params=5, units=['Ohms', 'F', '-', 'Ohms', 's']) +def RCSQ(p, f): + ''' + Beta element with CPE implementation + EIS: Randles circuit (CPE element) with spherical diffusion + + Notes + ----- + + .. math:: + + p[0] = Rct + p[1] = Qdl + p[2] = α + p[3] = Aw + p[4] = τd + + ''' + ω = np.array(f)*2*np.pi + Rct, Qdl, alpha, Aw, τd = p[0], p[1], p[2], p[3], p[4] + + Zd = Aw*np.tanh(np.sqrt(1j*ω*τd)) / \ + (np.sqrt(1j*ω*τd)-np.tanh(np.sqrt(1j*ω*τd))) + + tau = Rct*Qdl + Z = Rct/(Rct/(Rct+Zd) + tau*(1j*ω)**alpha) + return (Z) + + +@element(num_params=7, units=['Ohms', 'F', '-', 'Ohms', 's', '1/V', '-']) +def RCSQn(p, f): + ''' + Beta element with CPE implementation + 2nd-NLEIS: Randles circuit (CPE element) with spherical diffusion + + Notes + ----- + + .. math:: + + p[0] = Rct + p[1] = Qdl + p[2] = alpha + p[3] = Aw + p[4] = τd + p[5] = κ + p[6] = ε + ''' + + ω = np.array(f)*2*np.pi + Rct, Qdl, alpha, Aw, τd, κ, e = p[0], p[1], p[2], p[3], p[4], p[5], p[6] + + Zd1 = Aw*np.tanh(np.sqrt(1j*ω*τd)) / \ + (np.sqrt(1j*ω*τd)-np.tanh(np.sqrt(1j*ω*τd))) + Zd2 = Aw*np.tanh(np.sqrt(1j*2*ω*τd)) / \ + (np.sqrt(1j*2*ω*τd)-np.tanh(np.sqrt(1j*2*ω*τd))) + + tau = Rct*Qdl + y1 = Rct/(Zd1+Rct) + y2 = (Zd1/(Zd1+Rct)) + + Z1 = Rct/(y1+tau*(1j*ω)**alpha) + const = ((Rct*κ*y2**2)-Rct*e*F/(R*T)*y1**2)/(Zd2+Rct) + + Z2 = (const*Z1**2)/(tau*(1j*2*ω)**alpha+Rct/(Zd2+Rct)) + + return (Z2) + + +@element(num_params=5, units=['Ohms', 'F', '-', 'Ohms', 's']) +def RCDQ(p, f): + ''' + Beta element with CPE implementation + EIS: Randles circuit (CPE element) with spherical diffusion + + Notes + ----- + + .. math:: + + p[0] = Rct + p[1] = Qdl + p[2] = alpha + p[3] = Aw + p[4] = τd + + ''' + ω = np.array(f)*2*np.pi + Rct, Qdl, alpha, Aw, τd = p[0], p[1], p[2], p[3], p[4] + + Zd = Aw / (np.sqrt(1j*ω*τd) * np.tanh(np.sqrt(1j*ω*τd))) + + tau = Rct*Qdl + Z = Rct/(Rct/(Rct+Zd) + tau*(1j*ω)**alpha) + return (Z) + + +@element(num_params=7, units=['Ohms', 'F', '-', 'Ohms', 's', '1/V', '-']) +def RCDQn(p, f): + ''' + Beta element with CPE implementation + 2nd-NLEIS: Randles circuit (CPE element) with spherical diffusion + + Notes + ----- + + .. math:: + + p[0] = Rct + p[1] = Qdl + p[2] = alpha + p[3] = Aw + p[4] = τd + p[5] = κ + p[6] = ε + ''' + + ω = np.array(f)*2*np.pi + Rct, Qdl, alpha, Aw, τd, κ, e = p[0], p[1], p[2], p[3], p[4], p[5], p[6] + + Zd1 = Aw / (np.sqrt(1j*ω*τd) * np.tanh(np.sqrt(1j*ω*τd))) + Zd2 = Aw / (np.sqrt(1j*2*ω*τd) * np.tanh(np.sqrt(1j*2*ω*τd))) + + tau = Rct*Qdl + y1 = Rct/(Zd1+Rct) + y2 = (Zd1/(Zd1+Rct)) + + Z1 = Rct/(y1+tau*(1j*ω)**alpha) + const = ((Rct*κ*y2**2)-Rct*e*F/(R*T)*y1**2)/(Zd2+Rct) + + Z2 = (const*Z1**2)/(tau*(1j*2*ω)**alpha+Rct/(Zd2+Rct)) + + return (Z2) + + def get_element_from_name(name): excluded_chars = '0123456789_' return ''.join(char for char in name if char not in excluded_chars) diff --git a/nleis/nleis_tests/test_nleis_elements_pair.py b/nleis/nleis_tests/test_nleis_elements_pair.py index 447921a..5c69afc 100644 --- a/nleis/nleis_tests/test_nleis_elements_pair.py +++ b/nleis/nleis_tests/test_nleis_elements_pair.py @@ -114,6 +114,44 @@ def test_RC(): assert np.allclose(Z1, Z2) + # Test the convergence between + # Randles and Randles with CPE element (spherical diffusion) + + # EIS + freqs = [0.001, 1.0, 1000] + circuit_1 = CustomCircuit('RCS', initial_guess=[1, 1, 1, 1]) + circuit_2 = CustomCircuit('RCSQ', initial_guess=[1, 1, 1, 1, 1]) + Z1 = circuit_1.predict(freqs) + Z2 = circuit_2.predict(freqs) + assert np.allclose(Z1, Z2, atol=1e-3) + + # 2nd-NLEIS + circuit_1 = NLEISCustomCircuit('RCSn', initial_guess=[1, 1, 1, 1, 1, 1]) + circuit_2 = NLEISCustomCircuit('RCSQn', + initial_guess=[1, 1, 1, 1, 1, 1, 1]) + Z1 = circuit_1.predict(freqs) + Z2 = circuit_2.predict(freqs) + assert np.allclose(Z1, Z2, atol=1e-3) + + # Test the convergence between + # Randles and Randles with CPE element (Planar diffusion) + + # EIS + freqs = [0.001, 1.0, 1000] + circuit_1 = CustomCircuit('RCD', initial_guess=[1, 1, 1, 1]) + circuit_2 = CustomCircuit('RCDQ', initial_guess=[1, 1, 1, 1, 1]) + Z1 = circuit_1.predict(freqs) + Z2 = circuit_2.predict(freqs) + assert np.allclose(Z1, Z2, atol=1e-3) + + # 2nd-NLEIS + circuit_1 = NLEISCustomCircuit('RCDn', initial_guess=[1, 1, 1, 1, 1, 1]) + circuit_2 = NLEISCustomCircuit('RCDQn', + initial_guess=[1, 1, 1, 1, 1, 1, 1]) + Z1 = circuit_1.predict(freqs) + Z2 = circuit_2.predict(freqs) + assert np.allclose(Z1, Z2, atol=1e-3) + def test_d(): a = np.array([5 + 6 * 1j, 2 + 3 * 1j]) From 5f674758faca92b757979570f52f9a5e78f79cb3 Mon Sep 17 00:00:00 2001 From: Yuefan Ji Date: Thu, 9 Jan 2025 19:14:59 -0800 Subject: [PATCH 2/6] Code Improvement: Vectorized code for fast computation Test Fixes: update expected output in NLEIS tests to reflect correct units --- nleis/nleis_elements_pair.py | 901 +++++++++++++++----------------- nleis/nleis_tests/test_nleis.py | 8 +- 2 files changed, 414 insertions(+), 495 deletions(-) diff --git a/nleis/nleis_elements_pair.py b/nleis/nleis_elements_pair.py index 864c4c5..d2ff679 100644 --- a/nleis/nleis_elements_pair.py +++ b/nleis/nleis_elements_pair.py @@ -109,7 +109,7 @@ def RC(p, f): p[1] = C_{dl}; \\; """ - ω = np.array(f)*2*np.pi + ω = 2*np.pi * np.array(f) Rct, Cdl = p[0], p[1] ω_star = ω*Rct*Cdl @@ -154,7 +154,7 @@ def RCn(p, f): `_. ''' - ω = np.array(f)*2*np.pi + ω = 2*np.pi * np.array(f) Rct, Cdl, ε = p[0], p[1], p[2] ω_star = ω*Rct*Cdl @@ -206,12 +206,15 @@ def RCD(p, f): `_. ''' - ω = np.array(f)*2*np.pi + ω = 2*np.pi * np.array(f) Rct, Cdl, Aw, τd = p[0], p[1], p[2], p[3] - Zd = Aw / (np.sqrt(1j*ω*τd) * np.tanh(np.sqrt(1j*ω*τd))) - tau = ω*Rct*Cdl - Z = Rct / (Rct / (Rct + Zd) + 1j*tau) + sqrt_1j_ω_τd = np.sqrt(1j*ω*τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + Zd = Aw / (sqrt_1j_ω_τd * tanh_1j_ω_τd) + + ω_star = ω*Rct*Cdl + Z = Rct / (Rct / (Rct + Zd) + 1j*ω_star) return (Z) @@ -277,22 +280,27 @@ def RCDn(p, f): ''' - ω = np.array(f)*2*np.pi + ω = 2 * np.pi * np.array(f) Rct, Cdl, Aw, τd, κ, ε = p[0], p[1], p[2], p[3], p[4], p[5] - Zd1 = Aw / (np.sqrt(1j*ω*τd) * np.tanh(np.sqrt(1j*ω*τd))) - Zd2 = Aw / (np.sqrt(1j*2*ω*τd) * np.tanh(np.sqrt(1j*2*ω*τd))) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + Zd1 = Aw / (sqrt_1j_ω_τd * tanh_1j_ω_τd) - ω_star = ω*Rct*Cdl + sqrt_1j_2ω_τd = np.sqrt(1j * 2 * ω * τd) + tanh_1j_2ω_τd = np.tanh(sqrt_1j_2ω_τd) + Zd2 = Aw / (sqrt_1j_2ω_τd * tanh_1j_2ω_τd) + + ω_star = ω * Rct * Cdl y1 = Rct / (Zd1 + Rct) y2 = Zd1 / (Zd1 + Rct) - Z1 = Rct / (y1 + 1j*ω_star) - const = ((Rct*κ*y2**2) - Rct*ε*F/(R*T)*y1**2) / (Zd2 + Rct) + Z1 = Rct / (y1 + 1j * ω_star) + const = ((Rct * κ * y2**2) - Rct * ε * F / (R * T) * y1**2) / (Zd2 + Rct) - Z2 = (const*Z1**2) / (2*ω_star*1j + Rct/(Zd2+Rct)) + Z2 = (const * Z1**2) / (2 * ω_star * 1j + Rct / (Zd2 + Rct)) - return (Z2) + return Z2 @element(num_params=4, units=['Ohms', 'F', 'Ohms', 's']) @@ -341,15 +349,16 @@ def RCS(p, f): `_. ''' - ω = np.array(f)*2*np.pi + ω = 2 * np.pi * np.array(f) Rct, Cdl, Aw, τd = p[0], p[1], p[2], p[3] - Zd = Aw*np.tanh(np.sqrt(1j*ω*τd)) / \ - (np.sqrt(1j*ω*τd)-np.tanh(np.sqrt(1j*ω*τd))) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + Zd = Aw*tanh_1j_ω_τd / (sqrt_1j_ω_τd - tanh_1j_ω_τd) - ω_star = ω*Rct*Cdl - Z = Rct/(Rct/(Rct+Zd)+1j*ω_star) - return (Z) + ω_star = ω * Rct * Cdl + Z = Rct / (Rct / (Rct + Zd) + 1j * ω_star) + return Z @element(num_params=6, units=['Ohms', 'F', 'Ohms', 's', '1/V', '-']) @@ -417,24 +426,27 @@ def RCSn(p, f): ''' - ω = np.array(f)*2*np.pi + ω = 2 * np.pi * np.array(f) Rct, Cdl, Aw, τd, κ, ε = p[0], p[1], p[2], p[3], p[4], p[5] - Zd1 = Aw*np.tanh(np.sqrt(1j*ω*τd)) / \ - (np.sqrt(1j*ω*τd)-np.tanh(np.sqrt(1j*ω*τd))) - Zd2 = Aw*np.tanh(np.sqrt(1j*2*ω*τd)) / \ - (np.sqrt(1j*2*ω*τd)-np.tanh(np.sqrt(1j*2*ω*τd))) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + Zd1 = Aw * tanh_1j_ω_τd / (sqrt_1j_ω_τd - tanh_1j_ω_τd) - ω_star = ω*Rct*Cdl - y1 = Rct/(Zd1+Rct) - y2 = (Zd1/(Zd1+Rct)) + sqrt_1j_2ω_τd = np.sqrt(1j * 2 * ω * τd) + tanh_1j_2ω_τd = np.tanh(sqrt_1j_2ω_τd) + Zd2 = Aw * tanh_1j_2ω_τd / (sqrt_1j_2ω_τd - tanh_1j_2ω_τd) + + ω_star = ω * Rct * Cdl + y1 = Rct / (Zd1 + Rct) + y2 = Zd1 / (Zd1 + Rct) - Z1 = Rct/(y1+1j*ω_star) - const = ((Rct*κ*y2**2) - Rct*ε*F/(R*T)*y1**2)/(Zd2 + Rct) + Z1 = Rct / (y1 + 1j * ω_star) + const = ((Rct * κ * y2**2) - Rct * ε * F / (R * T) * y1**2) / (Zd2 + Rct) - Z2 = (const*Z1**2)/(2*ω_star*1j+Rct/(Zd2+Rct)) + Z2 = (const * Z1**2) / (2 * ω_star * 1j + Rct / (Zd2 + Rct)) - return (Z2) + return Z2 @element(num_params=3, units=['Ohms', 'Ohms', 'F']) @@ -475,13 +487,12 @@ def TP(p, f): ''' - ω = 2*np.pi*np.array(f) - + ω = 2 * np.pi * np.array(f) Rpore, Rct, Cdl = p[0], p[1], p[2] - beta = (1j*ω*Rpore*Cdl+Rpore/Rct)**(1/2) + beta = np.sqrt(1j * ω * Rpore * Cdl + Rpore / Rct) - Z = Rpore/(beta*np.tanh(beta)) + Z = Rpore / (beta * np.tanh(beta)) return Z @@ -536,40 +547,19 @@ def TPn(p, f): """ - ω = 2*np.pi*np.array(f) + ω = 2 * np.pi * np.array(f) Rpore, Rct, Cdl, ε = p[0], p[1], p[2], p[3] - b1 = (1j*ω*Rpore*Cdl + Rpore/Rct)**(1/2) - b2 = (1j*2*ω*Rpore*Cdl + Rpore/Rct)**(1/2) + b1 = np.sqrt(1j * ω * Rpore * Cdl + Rpore / Rct) + b2 = np.sqrt(1j * 2 * ω * Rpore * Cdl + Rpore / Rct) - sinh1 = [] - for x in b1: - if x < 100: - sinh1.append(np.sinh(x)) - else: - sinh1.append(1e10) - sinh2 = [] - cosh2 = [] - for x in b1: - if x < 100: - sinh2.append(np.sinh(2*x)) - cosh2.append(np.cosh(2*x)) - else: - sinh2.append(1e10) - cosh2.append(1e10) - sinh3 = [] - cosh3 = [] - for x in b2: - if x < 100: - sinh3.append(np.sinh(x)) - cosh3.append(np.cosh(x)) - else: - sinh3.append(1e10) - cosh3.append(1e10) + sinh1 = np.where(b1 < 100, np.sinh(b1), 1e10) + sinh2 = np.where(b1 < 100, np.sinh(2 * b1), 1e10) + cosh2 = np.where(b1 < 100, np.cosh(2 * b1), 1e10) - mf = ((Rpore**3) / Rct)*ε*(F/(R*T)) / ((b1*np.array(sinh1))**2) - part1 = (b1/b2)*np.array(sinh2) / ((b2**2-4*b1**2)*np.tanh(b2)) - part2 = -np.array(cosh2) / (2*(b2**2-4*b1**2)) - 1/(2*b2**2) - Z = mf*(part1 + part2) + mf = ((Rpore ** 3) / Rct) * ε * (F / (R * T)) / ((b1 * sinh1) ** 2) + part1 = (b1 / b2) * sinh2 / ((b2 ** 2 - 4 * b1 ** 2) * np.tanh(b2)) + part2 = -cosh2 / (2 * (b2 ** 2 - 4 * b1 ** 2)) - 1 / (2 * b2 ** 2) + Z = mf * (part1 + part2) return Z @@ -623,14 +613,17 @@ def TDP(p, f): `_. """ - ω = 2*np.pi*np.array(f) + ω = 2 * np.pi * np.array(f) Rpore, Rct, Cdl, Aw, τd = p[0], p[1], p[2], p[3], p[4] - Zd = Aw / (np.sqrt(1j*ω*τd) * np.tanh(np.sqrt(1j*ω*τd))) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + Zd = Aw / (sqrt_1j_ω_τd * tanh_1j_ω_τd) - beta = (1j*ω*Rpore*Cdl + Rpore/(Zd+Rct))**(1/2) - Z = Rpore / (beta*np.tanh(beta)) + beta = np.sqrt(1j * ω * Rpore * Cdl + Rpore / (Zd + Rct)) + tanh_beta = np.tanh(beta) + Z = Rpore / (beta * tanh_beta) return Z @@ -706,38 +699,31 @@ def TDPn(p, f): """ - ω = 2*np.pi*np.array(f) + ω = 2 * np.pi * np.array(f) + Rpore, Rct, Cdl, Aw, τd, κ, ε = p[0], p[1], p[2], p[3], p[4], p[5], p[6] - Zd1 = Aw / (np.sqrt(1j*ω*τd) * np.tanh(np.sqrt(1j*ω*τd))) - Zd2 = Aw / (np.sqrt(1j*2*ω*τd) * np.tanh(np.sqrt(1j*2*ω*τd))) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + sqrt_1j_2ω_τd = np.sqrt(1j * 2 * ω * τd) + + Zd1 = Aw / (sqrt_1j_ω_τd * np.tanh(sqrt_1j_ω_τd)) + Zd2 = Aw / (sqrt_1j_2ω_τd * np.tanh(sqrt_1j_2ω_τd)) y1 = Rct / (Zd1 + Rct) y2 = Zd1 / (Zd1 + Rct) - b1 = (1j*ω*Rpore*Cdl + Rpore/(Zd1+Rct))**(1/2) - b2 = (1j*2*ω*Rpore*Cdl + Rpore/(Zd2+Rct))**(1/2) + b1 = np.sqrt(1j * ω * Rpore * Cdl + Rpore / (Zd1 + Rct)) + b2 = np.sqrt(1j * 2 * ω * Rpore * Cdl + Rpore / (Zd2 + Rct)) - sinh1 = [] - for x in b1: - if x < 100: - sinh1.append(np.sinh(x)) - else: - sinh1.append(1e10) - sinh2 = [] - cosh2 = [] - for x in b1: - if x < 100: - sinh2.append(np.sinh(2*x)) - cosh2.append(np.cosh(2*x)) - else: - sinh2.append(1e10) - cosh2.append(1e10) - const = -((Rct*κ*y2**2) - Rct*ε*(F/(R*T))*y1**2) / (Zd2 + Rct) - mf = ((Rpore**3)*const/Rct) / ((b1*np.array(sinh1))**2) - part1 = (b1/b2)*np.array(sinh2) / ((b2**2-4*b1**2)*np.tanh(b2)) - part2 = -np.array(cosh2) / (2*(b2**2-4*b1**2)) - 1/(2*b2**2) - Z = mf*(part1 + part2) + sinh1 = np.where(np.abs(b1) < 100, np.sinh(b1), 1e10) + sinh2 = np.where(np.abs(b1) < 100, np.sinh(2 * b1), 1e10) + cosh2 = np.where(np.abs(b1) < 100, np.cosh(2 * b1), 1e10) + + const = -((Rct * κ * y2**2) - Rct * ε * (F / (R * T)) * y1**2)/(Zd2+Rct) + mf = ((Rpore**3) * const / Rct) / ((b1 * sinh1)**2) + part1 = (b1 / b2) * sinh2 / ((b2**2 - 4 * b1**2) * np.tanh(b2)) + part2 = -cosh2 / (2 * (b2**2 - 4 * b1**2)) - 1 / (2 * b2**2) + Z = mf * (part1 + part2) return Z @@ -792,14 +778,16 @@ def TDS(p, f): """ - ω = 2*np.pi*np.array(f) + ω = 2 * np.pi * np.array(f) Rpore, Rct, Cdl, Aw, τd = p[0], p[1], p[2], p[3], p[4] - Zd = Aw*np.tanh(np.sqrt(1j*ω*τd)) / \ - (np.sqrt(1j*ω*τd) - np.tanh(np.sqrt(1j*ω*τd))) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + Zd = Aw * tanh_1j_ω_τd / (sqrt_1j_ω_τd - tanh_1j_ω_τd) + + beta = np.sqrt(1j * ω * Rpore * Cdl + Rpore / (Zd + Rct)) - beta = (1j*ω*Rpore*Cdl + Rpore/(Zd+Rct))**(1/2) - Z = Rpore / (beta*np.tanh(beta)) + Z = Rpore / (beta * np.tanh(beta)) return Z @@ -878,39 +866,32 @@ def TDSn(p, f): """ - ω = 2*np.pi*np.array(f) + ω = 2 * np.pi * np.array(f) Rpore, Rct, Cdl, Aw, τd, κ, ε = p[0], p[1], p[2], p[3], p[4], p[5], p[6] - Zd1 = Aw*np.tanh(np.sqrt(1j*ω*τd)) / \ - (np.sqrt(1j*ω*τd) - np.tanh(np.sqrt(1j*ω*τd))) - Zd2 = Aw*np.tanh(np.sqrt(1j*2*ω*τd)) / \ - (np.sqrt(1j*2*ω*τd) - np.tanh(np.sqrt(1j*2*ω*τd))) - y1 = Rct / (Zd1+Rct) - y2 = Zd1 / (Zd1+Rct) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + sqrt_1j_2ω_τd = np.sqrt(1j * 2 * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + tanh_1j_2ω_τd = np.tanh(sqrt_1j_2ω_τd) - b1 = (1j*ω*Rpore*Cdl + Rpore/(Zd1+Rct))**(1/2) - b2 = (1j*2*ω*Rpore*Cdl + Rpore/(Zd2+Rct))**(1/2) + Zd1 = Aw * tanh_1j_ω_τd / (sqrt_1j_ω_τd - tanh_1j_ω_τd) + Zd2 = Aw * tanh_1j_2ω_τd / (sqrt_1j_2ω_τd - tanh_1j_2ω_τd) - sinh1 = [] - for x in b1: - if x < 100: - sinh1.append(np.sinh(x)) - else: - sinh1.append(1e10) - sinh2 = [] - cosh2 = [] - for x in b1: - if x < 100: - sinh2.append(np.sinh(2*x)) - cosh2.append(np.cosh(2*x)) - else: - sinh2.append(1e10) - cosh2.append(1e10) - const = -((Rct*κ*y2**2) - Rct*ε*(F/(R*T))*y1**2) / (Zd2+Rct) - mf = ((Rpore**3)*const/Rct) / ((b1*np.array(sinh1))**2) - part1 = (b1/b2)*np.array(sinh2) / ((b2**2-4*b1**2)*np.tanh(b2)) - part2 = -np.array(cosh2) / (2*(b2**2-4*b1**2)) - 1/(2*b2**2) - Z = mf*(part1+part2) + y1 = Rct / (Zd1 + Rct) + y2 = Zd1 / (Zd1 + Rct) + + b1 = np.sqrt(1j * ω * Rpore * Cdl + Rpore / (Zd1 + Rct)) + b2 = np.sqrt(1j * 2 * ω * Rpore * Cdl + Rpore / (Zd2 + Rct)) + + sinh1 = np.where(b1 < 100, np.sinh(b1), 1e10) + sinh2 = np.where(b1 < 100, np.sinh(2 * b1), 1e10) + cosh2 = np.where(b1 < 100, np.cosh(2 * b1), 1e10) + + const = -((Rct * κ * y2**2) - Rct * ε * (F / (R * T)) * y1**2)/(Zd2+Rct) + mf = ((Rpore**3) * const / Rct) / ((b1 * sinh1)**2) + part1 = (b1 / b2) * sinh2 / ((b2**2 - 4 * b1**2) * np.tanh(b2)) + part2 = -cosh2 / (2 * (b2**2 - 4 * b1**2)) - 1 / (2 * b2**2) + Z = mf * (part1 + part2) return Z @@ -961,21 +942,21 @@ def TDC(p, f): `_. """ - ω = 2*np.pi*np.array(f) + ω = 2 * np.pi * np.array(f) Rpore, Rct, Cdl, Aw, τd = p[0], p[1], p[2], p[3], p[4] - i01 = [] - i11 = [] - for x in np.sqrt(1j*ω*τd): - if x < 100: - i01.append(iv(0, x)) - i11.append(iv(1, x)) - else: - i01.append(1e20) - i11.append(1e20) - Zd = Aw*np.array(i01) / (np.sqrt(1j*ω*τd)*np.array(i11)) - beta = (1j*ω*Rpore*Cdl + Rpore/(Zd+Rct))**(1/2) - Z = Rpore / (beta*np.tanh(beta)) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + + i01 = iv(0, sqrt_1j_ω_τd) + i11 = iv(1, sqrt_1j_ω_τd) + + i01 = np.where(sqrt_1j_ω_τd < 100, i01, 1e20) + i11 = np.where(sqrt_1j_ω_τd < 100, i11, 1e20) + + Zd = Aw * i01 / (sqrt_1j_ω_τd * i11) + + beta = np.sqrt(1j * ω * Rpore * Cdl + Rpore / (Zd + Rct)) + Z = Rpore / (beta * np.tanh(beta)) return Z @@ -1051,60 +1032,72 @@ def TDCn(p, f): """ - ω = 2*np.pi*np.array(f) + ω = 2 * np.pi * np.array(f) Rpore, Rct, Cdl, Aw, τd, κ, ε = p[0], p[1], p[2], p[3], p[4], p[5], p[6] - i01 = [] - i11 = [] - for x in np.sqrt(1j*ω*τd): - if x < 100: - i01.append(iv(0, x)) - i11.append(iv(1, x)) - else: - i01.append(1e20) - i11.append(1e20) - i02 = [] - i12 = [] - for x in np.sqrt(1j*2*ω*τd): - if x < 100: - i02.append(iv(0, x)) - i12.append(iv(1, x)) - else: - i02.append(1e20) - i12.append(1e20) - Zd1 = Aw*np.array(i01) / (np.sqrt(1j*ω*τd)*np.array(i11)) - Zd2 = Aw*np.array(i02) / (np.sqrt(1j*2*ω*τd)*np.array(i12)) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + sqrt_1j_2ω_τd = np.sqrt(1j * 2 * ω * τd) - y1 = Rct / (Zd1+Rct) - y2 = Zd1 / (Zd1+Rct) + i01 = np.where(sqrt_1j_ω_τd < 100, iv(0, sqrt_1j_ω_τd), 1e20) + i11 = np.where(sqrt_1j_ω_τd < 100, iv(1, sqrt_1j_ω_τd), 1e20) + i02 = np.where(sqrt_1j_2ω_τd < 100, iv(0, sqrt_1j_2ω_τd), 1e20) + i12 = np.where(sqrt_1j_2ω_τd < 100, iv(1, sqrt_1j_2ω_τd), 1e20) - b1 = (1j*ω*Rpore*Cdl + Rpore/(Zd1+Rct))**(1/2) - b2 = (1j*2*ω*Rpore*Cdl + Rpore/(Zd2+Rct))**(1/2) + Zd1 = Aw * i01 / (sqrt_1j_ω_τd * i11) + Zd2 = Aw * i02 / (sqrt_1j_2ω_τd * i12) - sinh1 = [] - for x in b1: - if x < 100: - sinh1.append(np.sinh(x)) - else: - sinh1.append(1e10) - sinh2 = [] - cosh2 = [] - for x in b1: - if x < 100: - sinh2.append(np.sinh(2*x)) - cosh2.append(np.cosh(2*x)) - else: - sinh2.append(1e10) - cosh2.append(1e10) - const = -((Rct*κ*y2**2) - Rct*ε*(F/(R*T))*y1**2) / (Zd2+Rct) - mf = ((Rpore**3)*const/Rct) / ((b1*np.array(sinh1))**2) - part1 = (b1/b2)*np.array(sinh2) / ((b2**2-4*b1**2)*np.tanh(b2)) - part2 = -np.array(cosh2) / (2*(b2**2-4*b1**2)) - 1/(2*b2**2) - Z = mf*(part1+part2) + y1 = Rct / (Zd1 + Rct) + y2 = Zd1 / (Zd1 + Rct) + + b1 = np.sqrt(1j * ω * Rpore * Cdl + Rpore / (Zd1 + Rct)) + b2 = np.sqrt(1j * 2 * ω * Rpore * Cdl + Rpore / (Zd2 + Rct)) + + sinh1 = np.where(b1 < 100, np.sinh(b1), 1e10) + sinh2 = np.where(b1 < 100, np.sinh(2 * b1), 1e10) + cosh2 = np.where(b1 < 100, np.cosh(2 * b1), 1e10) + + const = -((Rct * κ * y2**2) - Rct * ε * (F / (R * T)) * y1**2) / (Zd2+Rct) + mf = ((Rpore**3) * const / Rct) / ((b1 * sinh1)**2) + part1 = (b1 / b2) * sinh2 / ((b2**2 - 4 * b1**2) * np.tanh(b2)) + part2 = -cosh2 / (2 * (b2**2 - 4 * b1**2)) - 1 / (2 * b2**2) + Z = mf * (part1 + part2) return Z -# TLM Model # +################################################################## +# TLM Model +################################################################## + + +def A_matrices_TLMn(N, Rpore, Z12t): + """ + Construct the matrix `Ax` for the TLMn model + Parameters + ---------- + N : int + Number of circuit elements + Rpore : float + Pore electrolyte resistance + Z12t : np.complex128 + The single element impedance at 2ω + Returns + ------- + Ax : np.ndarray + The matrix `Ax` for the TLMn model + """ + + Ax = np.zeros((N, N), dtype=np.complex128) + # Construct matrix `A` + for i in range(N - 1): + for j in range(N - 1 - i): + # construct the Rpore term + Ax[i, j] = (N - 1 - i - j) * Rpore + # construct the Rpore Z12t term + Ax[i, 0] += Z12t + Ax[i, N - 1 - i] -= Z12t + # construct the last row of the matrix + Ax[-1, :] = 1 + return (Ax) @element(num_params=6, units=['Ohm', 'Ohm', 'F', 'Ohm', 'F', '-']) @@ -1117,7 +1110,6 @@ def TLM(p, f): ----- - **Parameters:** .. math:: @@ -1134,27 +1126,25 @@ def TLM(p, f): """ - - N = int(p[5]) frequencies = np.array(f) + N = int(p[5]) - Rct = p[1]*N - Cdl = p[2]/N - Rpore = p[0]/N - Rs = p[3]*N - Cs = p[4]/N + Rct = p[1] * N + Cdl = p[2] / N + Rpore = p[0] / N + Rs = p[3] * N + Cs = p[4] / N Z1b = RC([Rct, Cdl], frequencies) Z1s = RC([Rs, Cs], frequencies) Zran = Z1b + Z1s - Req = Zran - for i in range(1, N): - Req_inv = (1/(Req+Rpore)) + 1/Zran - Req = 1 / Req_inv + Req = np.copy(Zran) + inv_Zran = 1 / Zran + for _ in range(1, N): + Req = 1 / ((1 / (Req + Rpore)) + inv_Zran) - return (Req) -### + return Req @element(num_params=8, units=['Ohm', 'Ohm', 'F', '-', 'Ohm', 'F', '-', '-']) @@ -1196,16 +1186,18 @@ def TLMn(p, f): `_. """ - I1 = mTi(p[0:6], f) # calculate the current fraction (1st harmonic) + # calculate the current fraction (1st harmonic) + I1 = mTi(p[0:6], f) - N = int(p[5]) frequencies = np.array(f) - Rpore = p[0]/N - Rct = p[1]*N - Cdl = p[2]/N - Rs = p[3]*N - Cs = p[4]/N + N = int(p[5]) + + Rpore = p[0] / N + Rct = p[1] * N + Cdl = p[2] / N + Rs = p[3] * N + Cs = p[4] / N εb = p[6] εs = p[7] @@ -1214,52 +1206,39 @@ def TLMn(p, f): Z1s = RC([Rs, Cs], frequencies) Z2b = RCn([Rct, Cdl, εb], frequencies) Z2s = RCn([Rs, Cs, εs], frequencies) - Z1b2t = RC([Rct, Cdl], 2*frequencies) - Z1s2t = RC([Rs, Cs], 2*frequencies) + Z1b2t = RC([Rct, Cdl], 2 * frequencies) + Z1s2t = RC([Rs, Cs], 2 * frequencies) + Z1 = Z1b + Z1s Z2 = Z2b + Z2s Z12t = Z1b2t + Z1s2t + if N == 1: - return (Z2) + return Z2 if N == 2: - sum1 = Z1**2 / (2*Z1+Rpore)**2 - sum2 = (Z12t*Rpore+Rpore**2) / ((2*Z12t+Rpore)*(2*Z1+Rpore)) - Z = (sum1+sum2)*Z2 - return (Z) - Z = np.zeros((len(frequencies)), dtype=complex) - for freq in range(0, len(frequencies)): - Ii = I1[freq] - - A = np.arange(N-1, 0, -1) - A1 = np.arange(N-1, 0, -1) - - for i in range(0, N-2): - for j in range(0, N-1-i): - A1[j] = A1[j]-1 - A = np.vstack((A, A1)) - A = A*Rpore - A = np.append(A, np.zeros((N-1, 1)), axis=1) - A = np.append(A, np.zeros((1, N)), axis=0) - A2 = np.zeros((N-1, N)) - for i in range(0, N-1): - A2[i, 0] += 1 - A2[i, N-1-i] -= 1 - A2 = np.vstack((A2, np.zeros(N))) - A2 = A2*Z12t[freq] - - A3 = np.vstack((np.zeros((N-1, N)), np.ones(N))) - - Ax = A2+A+A3 - - b = np.zeros((N, 1), dtype=complex) - - for i in range(0, N-1): + sum1 = Z1**2 / (2 * Z1 + Rpore)**2 + sum2 = (Z12t * Rpore + Rpore**2) / \ + ((2 * Z12t + Rpore) * (2 * Z1 + Rpore)) + Z = (sum1 + sum2) * Z2 + return Z + len_freq = len(frequencies) + Z = np.zeros(len_freq, dtype=np.complex128) + for freq_idx in range(len_freq): + Ii = I1[freq_idx] + # initialize the Ax and b matrix + Ax = A_matrices_TLMn(N, Rpore, Z12t[freq_idx]) + + b = np.zeros((N, 1), dtype=np.complex128) + + # construct the b matrix + for i in range(N - 1): b[i] = Ii[-1]**2 - Ii[i]**2 - I2 = np.linalg.solve(Ax, -b*Z2[freq]) - Z[freq] = Z2[freq]*Ii[0]**2 + I2[-1]*Z12t[freq] - return (Z) + I2 = np.linalg.solve(Ax, -b * Z2[freq_idx]) + Z[freq_idx] = Z2[freq_idx] * Ii[0]**2 + I2[-1] * Z12t[freq_idx] + + return Z @element(num_params=6, units=['Ohm', 'Ohm', 'F', 'Ohm', 'F', '-']) @@ -1311,35 +1290,34 @@ def mTi(p, f): Z1b = RC([Rct, Cdl], frequencies) Z1s = RC([Rs, Cs], frequencies) Zran = Z1b + Z1s - Req = Zran + Req = np.copy(Zran) + inv_Zran = 1 / Zran for i in range(1, N): - Req_inv = (1/(Req+Rpore))+1/Zran - Req = 1/Req_inv + Req = 1 / ((1 / (Req + Rpore)) + inv_Zran) Req = Req+Rpore - I1 = np.zeros((len(frequencies), N), dtype=complex) - for freq in range(0, len(frequencies)): - b1 = np.ones(N)*Req[freq] - - A = np.identity(N)*Zran[freq] + len_freq = len(frequencies) + I1 = np.zeros((len_freq, N), dtype=np.complex128) - A1 = np.ones((N, N))*Rpore + for freq_idx in range(0, len_freq): + Req_freq = Req[freq_idx] + Zran_freq = Zran[freq_idx] - for i in range(0, N): + # initialize the matrix and fill the diagonal with Zran + Ax = np.eye(N) * Zran_freq + # Get lower triangular indices of Ax matrix + i_idx, j_idx = np.tril_indices(N, -1) + # Fill lower triangular part of Ax matrix + Ax[i_idx, j_idx] = -(i_idx - j_idx) * Rpore + # Add the scaled Rpore matrix directly + # with addition to each row + Ax += np.arange(1, N + 1).reshape(-1, 1) * Rpore - A1[i, :] = A1[i, :]*(i+1) + b = np.ones(N)*Req_freq - for j in range(0, i): - - A[i][j] = -(i-j)*Rpore - - A = A+A1 - - b = b1 - - I1[freq, :] = np.linalg.solve(A, b) + I1[freq_idx, :] = np.linalg.solve(Ax, b) return (I1) @@ -1393,16 +1371,16 @@ def TLMS(p, f): τd = p[4] Rs = p[5]*N Cs = p[6]/N - Z1b = RCS([Rct, Cdl, Aw, τd], frequencies) Z1s = RC([Rs, Cs], frequencies) + Zran = Z1b + Z1s - Req = Zran - for i in range(1, N): + Req = np.copy(Zran) + inv_Zran = 1 / Zran + for _ in range(1, N): - Req_inv = (1/(Req+Rpore)) + 1/Zran - Req = 1/Req_inv + Req = 1 / ((1 / (Req + Rpore)) + inv_Zran) return (Req) @@ -1447,11 +1425,12 @@ def TLMSn(p, f): """ - I1 = mTiS(p[0:8], f) # calculate the current fraction (1st harmonic) + frequencies = np.array(f) - N = int(p[7]) + # calculate the current fraction (1st harmonic) + I1 = mTiS(p[0:8], frequencies) - frequencies = np.array(f) + N = int(p[7]) Rpore = p[0]/N Rct = p[1]*N @@ -1484,38 +1463,25 @@ def TLMSn(p, f): Z = (sum1+sum2)*Z2 return (Z) - Z = np.zeros(len(frequencies), dtype=complex) - for freq in range(0, len(frequencies)): - Ii = I1[freq] - - A = np.arange(N-1, 0, -1) - A1 = np.arange(N-1, 0, -1) - - for i in range(0, N-2): - for j in range(0, N-1-i): - A1[j] = A1[j]-1 - A = np.vstack((A, A1)) - A = A*Rpore - A = np.append(A, np.zeros((N-1, 1)), axis=1) - A = np.append(A, np.zeros((1, N)), axis=0) - A2 = np.zeros((N-1, N)) - for i in range(0, N-1): - A2[i, 0] += 1 - A2[i, N-1-i] -= 1 - A2 = np.vstack((A2, np.zeros(N))) - A2 = A2*Z12t[freq] - - A3 = np.vstack((np.zeros((N-1, N)), np.ones(N))) + len_freq = len(frequencies) + Z = np.zeros(len_freq, dtype=np.complex128) - Ax = A2+A+A3 + for freq_idx in range(len_freq): + Ii = I1[freq_idx] + Z12t_freq = Z12t[freq_idx] + Z2_freq = Z2[freq_idx] - b = np.zeros((N, 1), dtype=complex) + # construct the Ax matrix + Ax = A_matrices_TLMn(N, Rpore, Z12t_freq) + # initialize the b matrix + b = np.zeros((N, 1), dtype=np.complex128) - for i in range(0, N-1): + # construct the b matrix + for i in range(N-1): b[i] = Ii[-1]**2 - Ii[i]**2 - I2 = np.linalg.solve(Ax, -b*Z2[freq]) - Z[freq] = Z2[freq]*Ii[0]**2 + I2[-1]*Z12t[freq] + I2 = np.linalg.solve(Ax, -b*Z2_freq) + Z[freq_idx] = Z2_freq*Ii[0]**2 + I2[-1]*Z12t_freq return (Z) @@ -1559,7 +1525,7 @@ def mTiS(p, f): """ N = int(p[7]) - frequencies = np.array(f) + frequencies = f Rpore = p[0]/N Rct = p[1]*N @@ -1573,36 +1539,33 @@ def mTiS(p, f): Z1s = RC([Rs, Cs], frequencies) Zran = Z1b + Z1s - Req = Zran - for i in range(1, N): + Req = np.copy(Zran) + inv_Zran = 1 / Zran + for _ in range(1, N): - Req_inv = (1/(Req+Rpore))+1/Zran - Req = 1/Req_inv + Req = 1 / ((1 / (Req + Rpore)) + inv_Zran) Req = Req+Rpore - - I1 = np.zeros((len(frequencies), N), dtype=complex) - for freq in range(0, len(frequencies)): - b1 = np.ones(N)*Req[freq] - - A = np.identity(N)*Zran[freq] - - A1 = np.ones((N, N))*Rpore - - for i in range(0, N): - - A1[i, :] = A1[i, :]*(i+1) - - for j in range(0, i): - - A[i][j] = -(i-j)*Rpore - - A = A+A1 - - b = b1 - - I1[freq, :] = np.linalg.solve(A, b) - + len_freq = len(frequencies) + I1 = np.zeros((len_freq, N), dtype=np.complex128) + for freq_idx in range(0, len_freq): + Req_freq = Req[freq_idx] + Zran_freq = Zran[freq_idx] + + # initialize the matrix and fill the diagonal with Zran + Ax = np.eye(N) * Zran_freq + # Get lower triangular indices of Ax matrix + i_idx, j_idx = np.tril_indices(N, -1) + # Fill lower triangular part of Ax matrix + Ax[i_idx, j_idx] = -(i_idx - j_idx) * Rpore + # Add the scaled Rpore matrix directly + # with addition to each row + Ax += np.arange(1, N + 1).reshape(-1, 1) * Rpore + + # construct the b matrix + b = np.ones(N)*Req_freq + + I1[freq_idx, :] = np.linalg.solve(Ax, b) return (I1) @@ -1647,12 +1610,12 @@ def mTiSn(p, f): """ - I1 = mTiS(p[0:8], f) # calculate the current fraction (1st harmonic) + frequencies = np.array(f) + # calculate the current fraction (1st harmonic) + I1 = mTiS(p[0:8], frequencies) N = int(p[7]) - frequencies = np.array(f) - Rpore = p[0]/N Rct = p[1]*N Cdl = p[2]/N @@ -1671,52 +1634,32 @@ def mTiSn(p, f): Z2 = Z2b + Z2s Z12t = Z1b2t + Z1s2t - + len_freq = len(frequencies) if N == 1: return (0) + I2 = np.zeros((len_freq, N), dtype=np.complex128) + if N == 2: - I2 = np.zeros((len(frequencies), N), dtype=complex) - I2[0] = Z2*Rpore / (2*Z12t+Rpore)**2 - I2[1] = -Z2*Rpore / (2*Z12t+Rpore)**2 + I2[:, 0] = Z2*Rpore / (2*Z12t+Rpore)**2 + I2[:, 1] = -Z2*Rpore / (2*Z12t+Rpore)**2 return (I2) - I2 = np.zeros((len(frequencies), N), dtype=complex) - - for freq in range(0, len(frequencies)): - Ii = I1[freq] - - A = np.arange(N-1, 0, -1) - A1 = np.arange(N-1, 0, -1) - - for i in range(0, N-2): - for j in range(0, N-1-i): - A1[j] = A1[j]-1 - A = np.vstack((A, A1)) - A = A*Rpore - A = np.append(A, np.zeros((N-1, 1)), axis=1) - A = np.append(A, np.zeros((1, N)), axis=0) - A2 = np.zeros((N-1, N)) - for i in range(0, N-1): - A2[i, 0] += 1 - A2[i, N-1-i] -= 1 - A2 = np.vstack((A2, np.zeros(N))) - A2 = A2*Z12t[freq] - - A3 = np.vstack((np.zeros((N-1, N)), np.ones(N))) - - Ax = A2+A+A3 - - b = np.zeros((N, 1), dtype=complex) - + for freq_idx in range(0, len_freq): + Ii = I1[freq_idx] + # construct the Ax matrix + Ax = A_matrices_TLMn(N, Rpore, Z12t[freq_idx]) + # initialize and b matrix + b = np.zeros((N, 1), dtype=np.complex128) + # construct the b matrix for i in range(0, N-1): b[i] = Ii[-1]**2-Ii[i]**2 # reverse the order to display # the correct result from small to larger N - I2[freq, :] = np.linalg.solve(Ax, -b*Z2[freq]).flatten()[::-1] + I2[freq_idx, :] = np.linalg.solve(Ax, -b*Z2[freq_idx]).flatten()[::-1] return (I2) @@ -1775,11 +1718,12 @@ def TLMD(p, f): Z1s = RC([Rs, Cs], frequencies) Zran = Z1b + Z1s - Req = Zran - for i in range(1, N): + Req = np.copy(Zran) + inv_Zran = 1 / Zran + for _ in range(1, N): + + Req = 1 / ((1 / (Req + Rpore)) + inv_Zran) - Req_inv = (1/(Req+Rpore))+1/Zran - Req = 1/Req_inv return (Req) @@ -1859,38 +1803,24 @@ def TLMDn(p, f): sum2 = (Z12t*Rpore+Rpore**2) / ((2*Z12t+Rpore)*(2*Z1+Rpore)) Z = (sum1+sum2)*Z2 return (Z) - Z = np.zeros((len(frequencies)), dtype=complex) - for freq in range(0, len(frequencies)): - Ii = I1[freq] - - A = np.arange(N-1, 0, -1) - A1 = np.arange(N-1, 0, -1) - - for i in range(0, N-2): - for j in range(0, N-1-i): - A1[j] = A1[j]-1 - A = np.vstack((A, A1)) - A = A*Rpore - A = np.append(A, np.zeros((N-1, 1)), axis=1) - A = np.append(A, np.zeros((1, N)), axis=0) - A2 = np.zeros((N-1, N)) - for i in range(0, N-1): - A2[i, 0] += 1 - A2[i, N-1-i] -= 1 - A2 = np.vstack((A2, np.zeros(N))) - A2 = A2*Z12t[freq] - - A3 = np.vstack((np.zeros((N-1, N)), np.ones(N))) - - Ax = A2+A+A3 - - b = np.zeros((N, 1), dtype=complex) - + len_freq = len(frequencies) + Z = np.zeros((len_freq), dtype=np.complex128) + for freq_idx in range(len_freq): + Ii = I1[freq_idx] + Z12t_freq = Z12t[freq_idx] + Z2_freq = Z2[freq_idx] + + # construct the Ax matrix + Ax = A_matrices_TLMn(N, Rpore, Z12t_freq) + # initialize the b matrix + b = np.zeros((N, 1), dtype=np.complex128) + + # construct the b matrix for i in range(0, N-1): b[i] = Ii[-1]**2-Ii[i]**2 - I2 = np.linalg.solve(Ax, -b*Z2[freq]) - Z[freq] = Z2[freq]*Ii[0]**2 + I2[-1]*Z12t[freq] + I2 = np.linalg.solve(Ax, -b*Z2_freq) + Z[freq_idx] = Z2_freq*Ii[0]**2 + I2[-1]*Z12t_freq return (Z) @@ -1934,7 +1864,7 @@ def mTiD(p, f): """ N = int(p[7]) - frequencies = np.array(f) + frequencies = f Rpore = p[0]/N Rct = p[1]*N @@ -1948,35 +1878,33 @@ def mTiD(p, f): Z1s = RC([Rs, Cs], frequencies) Zran = Z1b + Z1s - Req = Zran - for i in range(1, N): + Req = np.copy(Zran) + inv_Zran = 1 / Zran + for _ in range(1, N): - Req_inv = (1/(Req+Rpore))+1/Zran - Req = 1/Req_inv + Req = 1 / ((1 / (Req + Rpore)) + inv_Zran) Req = Req+Rpore - I1 = np.zeros((len(frequencies), N), dtype=complex) - for freq in range(0, len(frequencies)): - b1 = np.ones(N)*Req[freq] - - A = np.identity(N)*Zran[freq] + I1 = np.zeros((len(frequencies), N), dtype=np.complex128) + for freq_idx in range(0, len(frequencies)): + Req_freq = Req[freq_idx] + Zran_freq = Zran[freq_idx] - A1 = np.ones((N, N))*Rpore + # initialize the matrix and fill the diagonal with Zran + Ax = np.eye(N) * Zran_freq + # Get lower triangular indices of Ax matrix + i_idx, j_idx = np.tril_indices(N, -1) + # Fill lower triangular part of Ax matrix + Ax[i_idx, j_idx] = -(i_idx - j_idx) * Rpore + # Add the scaled Rpore matrix directly + # with addition to each row + Ax += np.arange(1, N + 1).reshape(-1, 1) * Rpore - for i in range(0, N): + b = np.ones(N)*Req_freq - A1[i, :] = A1[i, :]*(i+1) + I1[freq_idx, :] = np.linalg.solve(Ax, b) - for j in range(0, i): - - A[i][j] = -(i-j)*Rpore - - A = A+A1 - - b = b1 - - I1[freq, :] = np.linalg.solve(A, b) return (I1) @@ -2047,49 +1975,29 @@ def mTiDn(p, f): if N == 1: return (0) - + len_freq = len(frequencies) + I2 = np.zeros((len_freq, N), dtype=np.complex128) if N == 2: - I2 = np.zeros((len(frequencies), N), dtype=complex) - I2[0] = Z2*Rpore / (2*Z12t+Rpore)**2 - I2[1] = -Z2*Rpore / (2*Z12t+Rpore)**2 + I2[:, 0] = Z2*Rpore / (2*Z12t+Rpore)**2 + I2[:, 1] = -Z2*Rpore / (2*Z12t+Rpore)**2 return (I2) + for freq_idx in range(len_freq): + Ii = I1[freq_idx] - I2 = np.zeros((len(frequencies), N), dtype=complex) + # construct the Ax matrix + Ax = A_matrices_TLMn(N, Rpore, Z12t[freq_idx]) + # initialize the b matrix + b = np.zeros((N, 1), dtype=np.complex128) - for freq in range(0, len(frequencies)): - Ii = I1[freq] - - A = np.arange(N-1, 0, -1) - A1 = np.arange(N-1, 0, -1) - - for i in range(0, N-2): - for j in range(0, N-1-i): - A1[j] = A1[j]-1 - A = np.vstack((A, A1)) - A = A*Rpore - A = np.append(A, np.zeros((N-1, 1)), axis=1) - A = np.append(A, np.zeros((1, N)), axis=0) - A2 = np.zeros((N-1, N)) - for i in range(0, N-1): - A2[i, 0] += 1 - A2[i, N-1-i] -= 1 - A2 = np.vstack((A2, np.zeros(N))) - A2 = A2*Z12t[freq] - - A3 = np.vstack((np.zeros((N-1, N)), np.ones(N))) - - Ax = A2+A+A3 - - b = np.zeros((N, 1), dtype=complex) - - for i in range(0, N-1): - b[i] = Ii[-1]**2-Ii[i]**2 + # construct the b matrix + for i in range(N-1): + b[i] = Ii[-1]**2 - Ii[i]**2 # reverse the order to display # the correct result from small to larger N - I2[freq, :] = np.linalg.solve(Ax, -b*Z2[freq]).flatten()[::-1] + I2[freq_idx, :] = np.linalg.solve(Ax, -b*Z2[freq_idx]).flatten()[::-1] return (I2) @@ -2112,15 +2020,16 @@ def RCSQ(p, f): p[4] = τd ''' - ω = np.array(f)*2*np.pi + ω = 2 * np.pi * np.array(f) Rct, Qdl, alpha, Aw, τd = p[0], p[1], p[2], p[3], p[4] - Zd = Aw*np.tanh(np.sqrt(1j*ω*τd)) / \ - (np.sqrt(1j*ω*τd)-np.tanh(np.sqrt(1j*ω*τd))) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + Zd = Aw * tanh_1j_ω_τd / (sqrt_1j_ω_τd - tanh_1j_ω_τd) - tau = Rct*Qdl - Z = Rct/(Rct/(Rct+Zd) + tau*(1j*ω)**alpha) - return (Z) + tau = Rct * Qdl + Z = Rct / (Rct / (Rct + Zd) + tau * (1j * ω) ** alpha) + return Z @element(num_params=7, units=['Ohms', 'F', '-', 'Ohms', 's', '1/V', '-']) @@ -2143,24 +2052,27 @@ def RCSQn(p, f): p[6] = ε ''' - ω = np.array(f)*2*np.pi - Rct, Qdl, alpha, Aw, τd, κ, e = p[0], p[1], p[2], p[3], p[4], p[5], p[6] + ω = 2 * np.pi * np.array(f) + Rct, Qdl, alpha, Aw, τd, κ, ε = p[0], p[1], p[2], p[3], p[4], p[5], p[6] + + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + sqrt_1j_2ω_τd = np.sqrt(1j * 2 * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + tanh_1j_2ω_τd = np.tanh(sqrt_1j_2ω_τd) - Zd1 = Aw*np.tanh(np.sqrt(1j*ω*τd)) / \ - (np.sqrt(1j*ω*τd)-np.tanh(np.sqrt(1j*ω*τd))) - Zd2 = Aw*np.tanh(np.sqrt(1j*2*ω*τd)) / \ - (np.sqrt(1j*2*ω*τd)-np.tanh(np.sqrt(1j*2*ω*τd))) + Zd1 = Aw * tanh_1j_ω_τd / (sqrt_1j_ω_τd - tanh_1j_ω_τd) + Zd2 = Aw * tanh_1j_2ω_τd / (sqrt_1j_2ω_τd - tanh_1j_2ω_τd) - tau = Rct*Qdl - y1 = Rct/(Zd1+Rct) - y2 = (Zd1/(Zd1+Rct)) + tau = Rct * Qdl + y1 = Rct / (Zd1 + Rct) + y2 = Zd1 / (Zd1 + Rct) - Z1 = Rct/(y1+tau*(1j*ω)**alpha) - const = ((Rct*κ*y2**2)-Rct*e*F/(R*T)*y1**2)/(Zd2+Rct) + Z1 = Rct / (y1 + tau * (1j * ω) ** alpha) + const = ((Rct * κ * y2**2) - Rct * ε * F / (R * T) * y1**2) / (Zd2 + Rct) - Z2 = (const*Z1**2)/(tau*(1j*2*ω)**alpha+Rct/(Zd2+Rct)) + Z2 = (const * Z1**2) / (tau * (1j * 2 * ω) ** alpha + Rct / (Zd2 + Rct)) - return (Z2) + return Z2 @element(num_params=5, units=['Ohms', 'F', '-', 'Ohms', 's']) @@ -2181,14 +2093,16 @@ def RCDQ(p, f): p[4] = τd ''' - ω = np.array(f)*2*np.pi + ω = 2 * np.pi * np.array(f) Rct, Qdl, alpha, Aw, τd = p[0], p[1], p[2], p[3], p[4] - Zd = Aw / (np.sqrt(1j*ω*τd) * np.tanh(np.sqrt(1j*ω*τd))) + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + Zd = Aw / (sqrt_1j_ω_τd * tanh_1j_ω_τd) - tau = Rct*Qdl - Z = Rct/(Rct/(Rct+Zd) + tau*(1j*ω)**alpha) - return (Z) + tau = Rct * Qdl + Z = Rct / (Rct / (Rct + Zd) + tau * (1j * ω) ** alpha) + return Z @element(num_params=7, units=['Ohms', 'F', '-', 'Ohms', 's', '1/V', '-']) @@ -2211,22 +2125,27 @@ def RCDQn(p, f): p[6] = ε ''' - ω = np.array(f)*2*np.pi - Rct, Qdl, alpha, Aw, τd, κ, e = p[0], p[1], p[2], p[3], p[4], p[5], p[6] + ω = 2 * np.pi * np.array(f) + Rct, Qdl, alpha, Aw, τd, κ, ε = p[0], p[1], p[2], p[3], p[4], p[5], p[6] + + sqrt_1j_ω_τd = np.sqrt(1j * ω * τd) + sqrt_1j_2ω_τd = np.sqrt(1j * 2 * ω * τd) + tanh_1j_ω_τd = np.tanh(sqrt_1j_ω_τd) + tanh_1j_2ω_τd = np.tanh(sqrt_1j_2ω_τd) - Zd1 = Aw / (np.sqrt(1j*ω*τd) * np.tanh(np.sqrt(1j*ω*τd))) - Zd2 = Aw / (np.sqrt(1j*2*ω*τd) * np.tanh(np.sqrt(1j*2*ω*τd))) + Zd1 = Aw / (sqrt_1j_ω_τd * tanh_1j_ω_τd) + Zd2 = Aw / (sqrt_1j_2ω_τd * tanh_1j_2ω_τd) - tau = Rct*Qdl - y1 = Rct/(Zd1+Rct) - y2 = (Zd1/(Zd1+Rct)) + tau = Rct * Qdl + y1 = Rct / (Zd1 + Rct) + y2 = Zd1 / (Zd1 + Rct) - Z1 = Rct/(y1+tau*(1j*ω)**alpha) - const = ((Rct*κ*y2**2)-Rct*e*F/(R*T)*y1**2)/(Zd2+Rct) + Z1 = Rct / (y1 + tau * (1j * ω) ** alpha) + const = ((Rct * κ * y2**2) - Rct * ε * F / (R * T) * y1**2) / (Zd2 + Rct) - Z2 = (const*Z1**2)/(tau*(1j*2*ω)**alpha+Rct/(Zd2+Rct)) + Z2 = (const * Z1**2) / (tau * (1j * 2 * ω) ** alpha + Rct / (Zd2 + Rct)) - return (Z2) + return Z2 def get_element_from_name(name): diff --git a/nleis/nleis_tests/test_nleis.py b/nleis/nleis_tests/test_nleis.py index ff2c726..7ad3acc 100644 --- a/nleis/nleis_tests/test_nleis.py +++ b/nleis/nleis_tests/test_nleis.py @@ -71,8 +71,8 @@ def test_EISandNLEIS(): 'TDSn0_6', 'TDSn1_0', 'TDSn1_1', 'TDSn1_2', 'TDSn1_3', 'TDSn1_4', 'TDSn1_5', 'TDSn1_6'] assert all_units_NLEIS == ['Ohms', 'Ohms', 'F', 'Ohms', 's', - '1/V', '', 'Ohms', 'Ohms', 'F', 'Ohms', 's', - '1/V', ''] + '1/V', '-', 'Ohms', 'Ohms', 'F', 'Ohms', 's', + '1/V', '-'] # check _is_fit() assert not NLEIS_circuit._is_fit() @@ -194,8 +194,8 @@ def test_NLEISCustomCircuit(): 'TDSn0_6', 'TDSn1_0', 'TDSn1_1', 'TDSn1_2', 'TDSn1_3', 'TDSn1_4', 'TDSn1_5', 'TDSn1_6'] assert all_units_NLEIS == ['Ohms', 'Ohms', 'F', 'Ohms', 's', - '1/V', '', 'Ohms', 'Ohms', 'F', 'Ohms', 's', - '1/V', ''] + '1/V', '-', 'Ohms', 'Ohms', 'F', 'Ohms', 's', + '1/V', '-'] # check _is_fit() assert not NLEIS_circuit._is_fit() From f752d85f50a4b3eb7dafe22a57c6eceeea071000 Mon Sep 17 00:00:00 2001 From: Yuefan Ji Date: Mon, 13 Jan 2025 20:27:33 -0800 Subject: [PATCH 3/6] Code Improvement: * Enabled circuit processing and circuit evaluation using execution graph * Updated and passed all the corresponding tests * Updated requirements.txt and setup.py * Bumped up version number * Fixed some docstrings --- nleis/fitting.py | 185 +++++++++++++++++++++++- nleis/nleis.py | 89 +++++++++--- nleis/nleis_elements_pair.py | 7 + nleis/nleis_fitting.py | 43 +++--- nleis/nleis_tests/test_model_io.py | 35 ++++- nleis/nleis_tests/test_nleis_fitting.py | 18 ++- requirements.txt | 1 + setup.py | 1 + 8 files changed, 333 insertions(+), 46 deletions(-) diff --git a/nleis/fitting.py b/nleis/fitting.py index 4a12c2d..98eb850 100644 --- a/nleis/fitting.py +++ b/nleis/fitting.py @@ -7,6 +7,9 @@ from impedance.models.circuits.elements import get_element_from_name from impedance.models.circuits.fitting import check_and_eval, rmse +import networkx as nx +import re + # Note: a lot of codes are directly adopted from impedance.py., # which is designed to be enable a easy integration in the future, # but now we are keep them separated to ensure the stable performance @@ -209,6 +212,7 @@ def set_default_bounds(circuit, constants={}): def circuit_fit(frequencies, impedances, circuit, initial_guess, constants={}, bounds=None, weight_by_modulus=False, global_opt=False, + graph=False, **kwargs): """ Main function for fitting an equivalent circuit to data. @@ -279,6 +283,8 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants={}, if bounds is None: bounds = set_default_bounds(circuit, constants=constants) + cg = CircuitGraph(circuit, constants) + if not global_opt: if 'maxfev' not in kwargs: kwargs['maxfev'] = int(1e5) @@ -289,10 +295,17 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants={}, 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) + if graph: + popt, pcov = curve_fit(cg.compute_long, f, + np.hstack([Z.real, Z.imag]), + p0=initial_guess, + bounds=bounds, + **kwargs, + ) + else: + 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. @@ -320,6 +333,9 @@ def opt_function(x): return rmse(wrapCircuit(circuit, constants)(f, *x), np.hstack([Z.real, Z.imag])) + def opt_function_graph(x): + return rmse(cg(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 @@ -337,8 +353,12 @@ def __call__(self, **kwargs): basinhopping_bounds = BasinhoppingBounds(xmin=bounds[0], xmax=bounds[1]) - results = basinhopping(opt_function, x0=initial_guess, - accept_test=basinhopping_bounds, **kwargs) + if graph: + results = basinhopping(opt_function_graph, x0=initial_guess, + accept_test=basinhopping_bounds, **kwargs) + else: + results = basinhopping(opt_function, x0=initial_guess, + accept_test=basinhopping_bounds, **kwargs) popt = results.x # Calculate perror @@ -583,3 +603,156 @@ def extract_circuit_elements(circuit): current_element.append(char) extracted_elements.append(''.join(current_element)) return extracted_elements + +# Circuit Graph for computation optimization +# Special Thanks to Jake Anderson for the original code + + +class CircuitGraph: + # regular expression to find parallel and difference blocks + _parallel_difference_block_expression = re.compile(r'(?:p|d)\([^()]*\)') + # _parallel_difference_block_expression = re.compile(r'p\([^()]*\)') + + # regular expression to remove whitespace + _whitespce = re.compile(r"\s+") + + def __init__(self, circuit, constants=None): + # remove all whitespace from the circuit string + self.circuit = self._whitespce.sub("", circuit) + # parse the circuit string and initialize the graph + self.parse_circuit() + # compute the execution order of the graph + self.execution_order = list(nx.topological_sort(self.graph)) + # initialize the constants dictionary + self.constants = constants if constants is not None else dict() + + def parse_circuit(self): + # initialize the node counters for each type of block + self.snum = 1 + self.pnum = 1 + self.dnum = 1 + # initialize the circuit string to be parsed + 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 + ] + # graph initialization + self.graph = nx.DiGraph() + # add nodes to the graph + for e, f in zip(element_name, element_func): + self.graph.add_node(e, Z=f) + + # find unnested parallel and difference blocks + pd_blocks = self._parallel_difference_block_expression.findall( + parsing_circuit) + + while len(pd_blocks) > 0: + # add parallel or difference blocks to the graph + # unnesting each time around the loop + for pd in pd_blocks: + operator = pd[0] + pd_elem = pd[2:-1].split(",") + if operator == "p": + pnode = f"p{self.pnum}" + self.pnum += 1 + self.graph.add_node(pnode, Z=circuit_elements["p"]) + for elem in pd_elem: + elem = self.add_series_elements(elem) + self.graph.add_edge(elem, pnode) + parsing_circuit = parsing_circuit.replace(pd, pnode) + elif operator == "d": + dnode = f"d{self.dnum}" + self.dnum += 1 + self.graph.add_node(dnode, Z=circuit_elements["d"]) + for elem in pd_elem: + elem = self.add_series_elements(elem) + self.graph.add_edge(elem, dnode) + parsing_circuit = parsing_circuit.replace(pd, dnode) + pd_blocks = self._parallel_difference_block_expression.findall( + parsing_circuit) + + # pick up any top line series connections + self.add_series_elements(parsing_circuit) + + # assign layers to the nodes + for layer, nodes in enumerate(nx.topological_generations(self.graph)): + for n in nodes: + self.graph.nodes[n]["layer"] = layer + # function to add series elements to the graph + + 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] + + # function to visualize the graph + def visualize_graph(self, **kwargs): + pos = nx.multipartite_layout(self.graph, subset_key="layer") + nx.draw_networkx(self.graph, pos=pos, **kwargs) + + # function to compute the impedance of the circuit + 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) + + return np.squeeze(node_results[node]) + + # To enable comparision + + def __eq__(self, other): + if not isinstance(other, CircuitGraph): + return False + # Compare the internal graph attributes + return (self.graph.nodes == other.graph.nodes + and self.graph.edges == other.graph.edges) + + # To enable direct calling + + def __call__(self, f, *parameters): + Z = self.compute(f, *parameters) + return np.hstack([Z.real, Z.imag]) + + def compute_long(self, f, *parameters): + Z = self.compute(f, *parameters) + return np.hstack([Z.real, Z.imag]) + + 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) + + +def format_parameter_name(name, j, n_params): + return f"{name}_{j}" if n_params > 1 else f"{name}" diff --git a/nleis/nleis.py b/nleis/nleis.py index a67acd7..3d04c9a 100644 --- a/nleis/nleis.py +++ b/nleis/nleis.py @@ -10,6 +10,8 @@ from impedance.visualization import plot_bode from .visualization import plot_altair, plot_first, plot_second +from nleis.fitting import CircuitGraph + import json import matplotlib.pyplot as plt @@ -20,7 +22,7 @@ class EISandNLEIS: # ToDO: add SSO method and SO method def __init__(self, circuit_1='', circuit_2='', initial_guess=[], - constants=None, name=None, **kwargs): + constants=None, name=None, graph=False, **kwargs): """ Constructor for a customizable linear and nonlinear equivalent circuit model @@ -69,6 +71,7 @@ def __init__(self, circuit_1='', circuit_2='', initial_guess=[], NLEIS: circuit_2 = 'd(TDSn0-TDSn1)' """ + self.graph = graph # 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: @@ -206,6 +209,10 @@ def __init__(self, circuit_1='', circuit_2='', initial_guess=[], self.p1, self.p2 = individual_parameters( self.edited_circuit, self.initial_guess, self.constants_1, self.constants_2) + if self.circuit_1: + self.cg1 = CircuitGraph(self.circuit_1, self.constants_1) + if self.circuit_2: + self.cg2 = CircuitGraph(self.circuit_2, self.constants_2) def __eq__(self, other): if self.__class__ == other.__class__: @@ -306,7 +313,7 @@ def fit(self, frequencies, Z1, Z2, bounds=None, constants_1=self.constants_1, constants_2=self.constants_2, bounds=bounds, opt=opt, cost=cost, max_f=max_f, param_norm=param_norm, - positive=positive, + positive=positive, graph=self.graph, **kwargs) # self.parameters_ = list(parameters) self.parameters_ = parameters @@ -317,6 +324,11 @@ def fit(self, frequencies, Z1, Z2, bounds=None, self.conf1, self.conf2 = individual_parameters( self.edited_circuit, self.conf_, self.constants_1, self.constants_2) + # if cov is not None: + # self.cov_ = cov + # self.cov1, self.cov2 = individual_parameters( + # self.edited_circuit, self.cov_, self.constants_1, + # self.constants_2) self.p1, self.p2 = individual_parameters( self.edited_circuit, self.parameters_, self.constants_1, @@ -379,14 +391,14 @@ def predict(self, frequencies, max_f=10, use_initial=False): x1, x2 = wrappedImpedance(self.edited_circuit, self.circuit_1, self.constants_1, self.circuit_2, self.constants_2, f1, f2, - self.parameters_) + self.parameters_, graph=self.graph) return x1, x2 else: warnings.warn("Simulating circuit based on initial parameters") x1, x2 = wrappedImpedance(self.edited_circuit, self.circuit_1, self.constants_1, self.circuit_2, self.constants_2, f1, f2, - self.initial_guess) + self.initial_guess, graph=self.graph) return x1, x2 def get_param_names(self, circuit, constants): @@ -744,6 +756,7 @@ def save(self, filepath): if self._is_fit(): parameters_ = list(self.parameters_) model_conf_ = list(self.conf_) + # model_cov_ = list(self.cov_) # parameters_ = self.parameters_ # model_conf_ = self.conf_ data_dict = {"Name": model_name, @@ -756,6 +769,7 @@ def save(self, filepath): "Fit": True, "Parameters": parameters_, "Confidence": model_conf_, + # "Covariance": model_cov_, "Edited Circuit Str": edited_circuit_str } else: @@ -815,6 +829,9 @@ def load(self, filepath, fitted_as_initial=False): self.edited_circuit, self.initial_guess, self.constants_1, self.constants_2) + self.cg1 = CircuitGraph(self.circuit_1, self.constants_1) + self.cg2 = CircuitGraph(self.circuit_2, self.constants_2) + self.name = model_name if json_data["Fit"]: @@ -833,12 +850,16 @@ def load(self, filepath, fitted_as_initial=False): self.p1, self.p2 = individual_parameters( self.edited_circuit, self.parameters_, self.constants_1, self.constants_2) + # self.cov_ = np.array(json_data["Covariance"]) + # self.cov1, self.cov2 = individual_parameters( + # self.edited_circuit, self.cov_, + # self.constants_1, self.constants_2) class NLEISCustomCircuit(BaseCircuit): # this class can be fully integrated into CustomCircuit in the future # , but for the stable performance of nleis.py, we overwrite it here - def __init__(self, circuit='', **kwargs): + def __init__(self, circuit='', graph=False, **kwargs): """ Constructor for a customizable nonlinear equivalent circuit model for NLEIS analysis. @@ -871,7 +892,11 @@ def __init__(self, circuit='', **kwargs): """ super().__init__(**kwargs) + # self.cov_ = None self.circuit = circuit.replace(" ", "") + self.graph = graph + if self.circuit: + self.cg = CircuitGraph(self.circuit, self.constants) circuit_len = calculateCircuitLength(self.circuit) @@ -942,15 +967,20 @@ def fit(self, frequencies, impedance, bounds=None, impedance = impedance[mask] 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 = \ + circuit_fit(frequencies, impedance, + self.circuit, + self.initial_guess, + constants=self.constants, + bounds=bounds, + weight_by_modulus=weight_by_modulus, + graph=self.graph, + **kwargs) self.parameters_ = parameters if conf is not None: self.conf_ = conf + # if cov is not None: + # self.cov_ = cov else: raise ValueError('No initial guess supplied') @@ -982,20 +1012,31 @@ def predict(self, frequencies, max_f=10, use_initial=False): frequencies = np.array(frequencies, dtype=float) mask = np.array(frequencies) < max_f frequencies = frequencies[mask] + if self.graph: + self.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) + if self.graph: + return self.cg.compute(frequencies, *self.parameters_) + else: + return eval(buildCircuit(self.circuit, frequencies, + *self.parameters_, + constants=self.constants, + eval_string='', + index=0)[0], + circuit_elements) 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) + + if self.graph: + return self.cg.compute(frequencies, *self.initial_guess) + else: + return eval(buildCircuit(self.circuit, frequencies, + *self.initial_guess, + constants=self.constants, + eval_string='', + index=0)[0], + circuit_elements) def get_param_names(self): """ @@ -1211,3 +1252,9 @@ def plot(self, ax=None, f_data=None, Z2_data=None, else: raise ValueError("Kind must be one of 'altair'," + f"'nyquist', or 'bode' (received {kind})") + + # add on to the load function to create the graph + + def load(self, filepath, fitted_as_initial=False): + super().load(filepath, fitted_as_initial) + self.cg = CircuitGraph(self.circuit, self.constants) diff --git a/nleis/nleis_elements_pair.py b/nleis/nleis_elements_pair.py index 018e359..fb6503d 100644 --- a/nleis/nleis_elements_pair.py +++ b/nleis/nleis_elements_pair.py @@ -51,6 +51,13 @@ 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 return decorator diff --git a/nleis/nleis_fitting.py b/nleis/nleis_fitting.py index c3f403c..8d3e4c7 100644 --- a/nleis/nleis_fitting.py +++ b/nleis/nleis_fitting.py @@ -10,6 +10,7 @@ from .fitting import set_default_bounds, buildCircuit, extract_circuit_elements from scipy.optimize import minimize import warnings +from nleis.fitting import CircuitGraph # Customize warning format (here, simpler and just the message) warnings.formatwarning = lambda message, category, filename, lineno, \ @@ -76,7 +77,7 @@ def data_processing(f, Z1, Z2, max_f=10): def simul_fit(frequencies, Z1, Z2, circuit_1, circuit_2, edited_circuit, initial_guess, constants_1={}, constants_2={}, bounds=None, opt='max', cost=0.5, max_f=10, param_norm=True, - positive=True, **kwargs): + positive=True, graph=False, **kwargs): """ Main function for the simultaneous fitting of EIS and 2nd-NLEIS data. @@ -234,7 +235,7 @@ def simul_fit(frequencies, Z1, Z2, circuit_1, circuit_2, edited_circuit, popt, pcov = curve_fit( wrapCircuit_simul(edited_circuit, circuit_1, constants_1, circuit_2, constants_2, - ub, max_f), frequencies, Zstack, + ub, max_f, graph=graph), frequencies, Zstack, p0=initial_guess, bounds=bounds, **kwargs) # Calculate one standard deviation error estimates for fit parameters, @@ -259,7 +260,7 @@ def simul_fit(frequencies, Z1, Z2, circuit_1, circuit_2, edited_circuit, wrapNeg_log_likelihood(frequencies, Z1, Z2, edited_circuit, circuit_1, constants_1, circuit_2, constants_2, - ub, max_f, cost=cost), + ub, max_f, cost=cost, graph=graph), x0=initial_guess, bounds=bounds, **kwargs) return (res.x*ub, None) @@ -267,7 +268,8 @@ def simul_fit(frequencies, Z1, Z2, circuit_1, circuit_2, edited_circuit, def wrapNeg_log_likelihood(frequencies, Z1, Z2, edited_circuit, circuit_1, constants_1, - circuit_2, constants_2, ub, max_f=10, cost=0.5): + circuit_2, constants_2, ub, max_f=10, cost=0.5, + graph=False): ''' wraps function so we can pass the circuit string for negtive log likelihood optimization''' @@ -298,7 +300,9 @@ def wrappedNeg_log_likelihood(parameters): f2 = frequencies[mask] x1, x2 = wrappedImpedance(edited_circuit, circuit_1, constants_1, circuit_2, - constants_2, f1, f2, parameters*ub) + constants_2, f1, f2, parameters*ub, + graph=graph) + # No normalization in currently applied # Z1max = max(np.abs(Z1)) # Z2max = max(np.abs(Z2)) @@ -314,7 +318,7 @@ def wrappedNeg_log_likelihood(parameters): def wrapCircuit_simul(edited_circuit, circuit_1, constants_1, circuit_2, - constants_2, ub, max_f=10): + constants_2, ub, max_f=10, graph=False): """ wraps function so we can pass the circuit string for simultaneous fitting """ def wrappedCircuit_simul(frequencies, *parameters): @@ -343,7 +347,7 @@ def wrappedCircuit_simul(frequencies, *parameters): x1, x2 = wrappedImpedance(edited_circuit, circuit_1, constants_1, circuit_2, constants_2, - f1, f2, parameters*ub) + f1, f2, parameters*ub, graph=graph) y1_real = np.real(x1) y1_imag = np.imag(x1) @@ -357,7 +361,7 @@ def wrappedCircuit_simul(frequencies, *parameters): def wrappedImpedance(edited_circuit, circuit_1, constants_1, circuit_2, - constants_2, f1, f2, parameters): + constants_2, f1, f2, parameters, graph=False): """ Calculate EIS and 2nd-NLEIS impedances using the provided circuits. @@ -402,15 +406,20 @@ def wrappedImpedance(edited_circuit, circuit_1, constants_1, circuit_2, p1, p2 = individual_parameters( edited_circuit, parameters, constants_1, constants_2) - - x1 = eval(buildCircuit(circuit_1, f1, *p1, - constants=constants_1, eval_string='', - index=0)[0], - circuit_elements) - x2 = eval(buildCircuit(circuit_2, f2, *p2, - constants=constants_2, eval_string='', - index=0)[0], - circuit_elements) + if graph: + graph_EIS = CircuitGraph(circuit_1, constants_1) + graph_NLEIS = CircuitGraph(circuit_2, constants_2) + x1 = graph_EIS.compute(f1, *p1) + x2 = graph_NLEIS.compute(f2, *p2) + else: + x1 = eval(buildCircuit(circuit_1, f1, *p1, + constants=constants_1, eval_string='', + index=0)[0], + circuit_elements) + x2 = eval(buildCircuit(circuit_2, f2, *p2, + constants=constants_2, eval_string='', + index=0)[0], + circuit_elements) return (x1, x2) diff --git a/nleis/nleis_tests/test_model_io.py b/nleis/nleis_tests/test_model_io.py index 2be413c..39832a4 100644 --- a/nleis/nleis_tests/test_model_io.py +++ b/nleis/nleis_tests/test_model_io.py @@ -1,5 +1,5 @@ import numpy as np -from nleis.nleis import EISandNLEIS +from nleis.nleis import EISandNLEIS, NLEISCustomCircuit import os import os.path @@ -20,6 +20,7 @@ def test_model_io(): circ_str_1 = 'L0-R0-TDS0-TDS1' circ_str_2 = 'd(TDSn0,TDSn1)' + # test for EISandNLEIS initial_guess = [1e-7, 1e-3, # L0,RO 5e-3, 1e-3, 10, 1e-2, 100, 10, 0.1, # TDS0 + additioal nonlinear parameters @@ -53,3 +54,35 @@ def test_model_io(): circuit_1 = EISandNLEIS(circ_str_1, circ_str_2, initial_guess=p_fit) assert str(circuit_1) == str(fitted_template) assert circuit_1 == fitted_template + + # test for NLIESCustomCircuit + initial_guess = [ + 5e-3, 1e-3, 10, 1e-2, 100, 10, 0.1, + # TDS0 + additioal nonlinear parameters + 1e-3, 1e-3, 1e-3, 1e-2, 1000, 0, 0 + # TDS1 + additioal nonlinear parameters + ] + + circuit_1 = NLEISCustomCircuit(circ_str_2, + initial_guess=initial_guess) + + circuit_1.save(os.path.join(data_dir, 'test_io.json')) + + circuit_2 = NLEISCustomCircuit() + circuit_2.load(os.path.join(data_dir, 'test_io.json')) + + assert circuit_1 == circuit_2 + + circuit_1.fit(frequencies, Z2) + p_fit = list(circuit_1.parameters_) + circuit_1.save(os.path.join(data_dir, 'test_io.json')) + circuit_2 = NLEISCustomCircuit() + circuit_2.load(os.path.join(data_dir, 'test_io.json')) + + assert str(circuit_1) == str(circuit_2) + assert circuit_1 == circuit_2 + + fitted_template = NLEISCustomCircuit() + fitted_template.load(os.path.join( + data_dir, 'test_io.json'), fitted_as_initial=True) + circuit_1 = NLEISCustomCircuit(circ_str_2, initial_guess=p_fit) diff --git a/nleis/nleis_tests/test_nleis_fitting.py b/nleis/nleis_tests/test_nleis_fitting.py index 31ea4f9..c09e26e 100644 --- a/nleis/nleis_tests/test_nleis_fitting.py +++ b/nleis/nleis_tests/test_nleis_fitting.py @@ -132,9 +132,15 @@ def test_set_default_bounds(): default_bounds[1][6] = 0.5 default_bounds[1][7] = 0.5 assert np.allclose(default_bounds, bounds_from_func) -# This test remain unchanged as we are adopting it from impedance.py + circuit = 'RCSQ' + bounds_from_func = set_default_bounds(circuit, constants=constants) + default_bounds = (np.zeros(5), np.inf*np.ones(5)) + default_bounds[1][2] = 1 + assert np.allclose(default_bounds, bounds_from_func) + +# This test remain unchanged as we are adopting it from impedance.py def test_circuit_fit(): # Test trivial model (10 Ohm resistor) @@ -298,6 +304,16 @@ def test_buildCircuit(): 's([R([100],[1000.0,5.0,0.01]),' +\ 'R([1],[1000.0,5.0,0.01])])' + # Test constants in circuit + circuit = 'Wo1' + params = [100] + frequencies = [1000.0, 5.0, 0.01] + constants = {'Wo1_1': 1} + + assert buildCircuit(circuit, frequencies, *params, + constants=constants)[0].replace(' ', '') == \ + 'Wo([100,1],[1000.0,5.0,0.01])' + def test_mae(): a = np.array([2 + 4*1j, 3 + 2*1j]) diff --git a/requirements.txt b/requirements.txt index 1674d92..66ded7e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,3 +10,4 @@ scipy>=1.0 pandas >= 2.0.2 sphinx_rtd_theme impedance >= 1.7.1 +networkx>=2.6.3 diff --git a/setup.py b/setup.py index 6dc3927..1a1079f 100644 --- a/setup.py +++ b/setup.py @@ -31,5 +31,6 @@ python_requires=">=3.8", install_requires=['altair>=3.0', 'matplotlib>=3.5', 'numpy>=1.14', 'scipy>=1.0', + 'networkx>=2.6.3', 'impedance>=1.7.1', 'pandas >= 2.0.2'], ) From 2429eff3419ae0a0eed05a46b01661d4f0e83b13 Mon Sep 17 00:00:00 2001 From: Yuefan Ji Date: Mon, 13 Jan 2025 20:30:48 -0800 Subject: [PATCH 4/6] Added more test and fixed docstrings --- nleis/__init__.py | 2 +- nleis/fitting.py | 5 +- nleis/nleis.py | 13 ++- nleis/nleis_fitting.py | 28 +++++- nleis/nleis_tests/__init__.py | 2 +- nleis/nleis_tests/test_graph.py | 152 ++++++++++++++++++++++++++++++++ nleis/nleis_tests/test_nleis.py | 49 +--------- 7 files changed, 197 insertions(+), 54 deletions(-) create mode 100644 nleis/nleis_tests/test_graph.py diff --git a/nleis/__init__.py b/nleis/__init__.py index 381f6f1..17ea3e6 100644 --- a/nleis/__init__.py +++ b/nleis/__init__.py @@ -1,4 +1,4 @@ -__version__ = "0.1.1" +__version__ = "0.2" try: from .nleis import * # noqa: F401, F403 diff --git a/nleis/fitting.py b/nleis/fitting.py index 98eb850..43cb0a6 100644 --- a/nleis/fitting.py +++ b/nleis/fitting.py @@ -257,6 +257,10 @@ def circuit_fit(frequencies, impedances, circuit, initial_guess, constants={}, If global optimization should be used (uses the basinhopping algorithm). Defaults to False + graph : bool, optional + Whether to use execution graph to process the circuit. + Defaults to False, which uses eval based code + kwargs : Keyword arguments passed to scipy.optimize.curve_fit or scipy.optimize.basinhopping @@ -611,7 +615,6 @@ def extract_circuit_elements(circuit): class CircuitGraph: # regular expression to find parallel and difference blocks _parallel_difference_block_expression = re.compile(r'(?:p|d)\([^()]*\)') - # _parallel_difference_block_expression = re.compile(r'p\([^()]*\)') # regular expression to remove whitespace _whitespce = re.compile(r"\s+") diff --git a/nleis/nleis.py b/nleis/nleis.py index 3d04c9a..8ef5f0f 100644 --- a/nleis/nleis.py +++ b/nleis/nleis.py @@ -48,6 +48,9 @@ def __init__(self, circuit_1='', circuit_2='', initial_guess=[], name : str, optional A name for the model. + graph : bool, optional + Whether to use execution graph to process the circuit. + Defaults to False, which uses eval based code Notes ----- @@ -872,6 +875,10 @@ def __init__(self, circuit='', graph=False, **kwargs): circuit : str A string representing the nonlinear equivalent circuit for NLEIS. + graph : bool, optional + Whether to use execution graph to process the circuit. + Defaults to False, which uses eval based code + Notes ----- A custom NLEIS circuit is defined as a string comprised of elements @@ -909,7 +916,7 @@ def __init__(self, circuit='', graph=False, **kwargs): f'the circuit length ({circuit_len})') def fit(self, frequencies, impedance, bounds=None, - weight_by_modulus=False, max_f=10, **kwargs): + weight_by_modulus=False, max_f=np.inf, **kwargs): """ Fit the nonlinear equivalent circuit model to NLEIS data. @@ -986,7 +993,7 @@ def fit(self, frequencies, impedance, bounds=None, return self - def predict(self, frequencies, max_f=10, use_initial=False): + def predict(self, frequencies, max_f=np.inf, use_initial=False): """ Predict impedance using an nonlinear equivalent circuit model @@ -1097,7 +1104,7 @@ def extract(self): return dict def plot(self, ax=None, f_data=None, Z2_data=None, - kind='nyquist', max_f=10, **kwargs): + kind='nyquist', max_f=np.inf, **kwargs): """ Visualizes the model and optional data as Nyquist, Bode, or Altair (interactive) plots. diff --git a/nleis/nleis_fitting.py b/nleis/nleis_fitting.py index 8d3e4c7..582740d 100644 --- a/nleis/nleis_fitting.py +++ b/nleis/nleis_fitting.py @@ -160,6 +160,10 @@ def simul_fit(frequencies, Z1, Z2, circuit_1, circuit_2, edited_circuit, positive : bool, optional If True, high-frequency inductance is eliminated. Defaults to True. + graph : bool, optional + Whether to use execution graph to process the circuit. + Defaults to False, which uses eval based code + kwargs : Additional keyword arguments passed to `scipy.optimize.curve_fit`. @@ -282,12 +286,21 @@ def wrappedNeg_log_likelihood(parameters): frequencies: list of floats Z1: EIS data Z2: NLEIS data + edited_circuit : string circuit_1 : string constants_1 : dict circuit_2 : string constants_2 : dict ub : list of floats upper bound if bounds are provided max_f: int + cost : float, optional + Weighting between EIS and 2nd-NLEIS data. A value greater than 0.5 + puts more weight on EIS, + while a value less than 0.5 puts more weight on 2nd-NLEIS. + Default is 0.5. + graph : bool, optional + Whether to use execution graph to process the circuit. + Defaults to False, which uses eval based code parameters : list of floats Returns @@ -327,14 +340,22 @@ def wrappedCircuit_simul(frequencies, *parameters): Parameters ---------- + edited_circuit : string circuit_1 : string constants_1 : dict circuit_2 : string constants_2 : dict - max_f: int - parameters : list of floats + ub : list of floats upper bound if bounds are provided + max_f: float + graph : bool, optional + Whether to use execution graph to process the circuit. + Defaults to False, which uses eval based code + frequencies : list of floats + parameters : list of floats + + Returns ------- array of floats @@ -394,6 +415,9 @@ def wrappedImpedance(edited_circuit, circuit_1, constants_1, circuit_2, parameters : list of float Full set of parameters derived from the edited circuit string. + graph : bool, optional + Whether to use execution graph to process the circuit. + Defaults to False, which uses eval based code Returns ------- diff --git a/nleis/nleis_tests/__init__.py b/nleis/nleis_tests/__init__.py index 485f44a..0988857 100644 --- a/nleis/nleis_tests/__init__.py +++ b/nleis/nleis_tests/__init__.py @@ -1 +1 @@ -__version__ = "0.1.1" +__version__ = "0.2" diff --git a/nleis/nleis_tests/test_graph.py b/nleis/nleis_tests/test_graph.py new file mode 100644 index 0000000..ce2b5a3 --- /dev/null +++ b/nleis/nleis_tests/test_graph.py @@ -0,0 +1,152 @@ +import numpy as np +import matplotlib.pyplot as plt # noqa: F401 + +import timeit + +from nleis.nleis import EISandNLEIS, NLEISCustomCircuit # noqa: F401 +from nleis.fitting import CircuitGraph + + +def test_graph_NLEISCustomCircuit(): + circ_str = 'd(TDSn0,TDSn1)' + initial_guess = [ + 5e-3, 1e-3, 10, 1e-2, 100, 10, 0.1, + # TDS0 + additioal nonlinear parameters + 1e-3, 1e-3, 1e-3, 1e-2, 1000, 1, 0.01, + # TDS1 + additioal nonlinear parameters + ] + # initialize + eval_circuit = NLEISCustomCircuit( + circ_str, initial_guess=initial_guess, graph=False) + graph_circuit = NLEISCustomCircuit( + circ_str, initial_guess=initial_guess, graph=True) + + frequencies = np.geomspace(1e-3, 1e1, 100) + + def time_predict_eval(): + return eval_circuit.predict(frequencies, max_f=np.inf) + + def time_predict_graph(): + return graph_circuit.predict(frequencies, max_f=np.inf) + eval_time = timeit.timeit(time_predict_eval, number=10) + graph_time = timeit.timeit(time_predict_graph, number=10) + Z2_eval = eval_circuit.predict(frequencies, max_f=np.inf) + Z2_graph = graph_circuit.predict(frequencies, max_f=np.inf) + assert np.allclose(Z2_eval, Z2_graph) + assert graph_time < eval_time + print(f'eval_time: {eval_time}, graph_time: {graph_time}') + + # test the fitting of the graph circuit using curve_fit + graph_circuit.fit(frequencies, Z2_graph, max_f=np.inf) + assert np.allclose(graph_circuit.predict( + frequencies, max_f=np.inf), Z2_graph) + + # test the fitting of the graph circuit using with global_opt = True + graph_circuit.fit(frequencies, Z2_graph, max_f=np.inf, global_opt=True) + assert np.allclose(graph_circuit.predict( + frequencies, max_f=np.inf), Z2_graph) + + +def test_graph_EISandNLEIS(): + circ_str_1 = 'L0-R0-TDS0-TDS1' + circ_str_2 = 'd(TDSn0,TDSn1)' + initial_guess = [1e-7, 1e-3, # L0,RO + 5e-3, 1e-3, 10, 1e-2, 100, 10, 0.1, + # TDS0 + additioal nonlinear parameters + 1e-3, 1e-3, 1e-3, 1e-2, 1000, 1, 0.01, + # TDS1 + additioal nonlinear parameters + ] + + eval_circuit = EISandNLEIS( + circ_str_1, circ_str_2, initial_guess=initial_guess, graph=False) + graph_circuit = EISandNLEIS(circ_str_1, circ_str_2, + initial_guess=initial_guess, graph=True) + + frequencies = np.geomspace(1e-3, 1e3, 100) + + def time_predict_eval(): + return eval_circuit.predict(frequencies) + + def time_predict_graph(): + return graph_circuit.predict(frequencies) + eval_time = timeit.timeit(time_predict_eval, number=10) + graph_time = timeit.timeit(time_predict_graph, number=10) + + Z1_eval, Z2_eval = eval_circuit.predict(frequencies, max_f=np.inf) + Z1_graph, Z2_graph = graph_circuit.predict(frequencies, max_f=np.inf) + assert np.allclose(Z1_eval, Z1_graph) + assert np.allclose(Z2_eval, Z2_graph) + + assert graph_time < eval_time + print(f'eval_time: {eval_time}, graph_time: {graph_time}') + + # test the fitting of the graph circuit using curve_fit and opt = 'max' + graph_circuit.fit(frequencies, Z1_graph, Z2_graph, max_f=np.inf) + assert np.allclose(graph_circuit.predict( + frequencies, max_f=np.inf), (Z1_graph, Z2_graph)) + + # test the fitting of the graph circuit using curve_fit and opt = 'neg' + graph_circuit.fit(frequencies, Z1_graph, Z2_graph, opt='neg', max_f=np.inf) + assert np.allclose(graph_circuit.predict( + frequencies, max_f=np.inf), (Z1_graph, Z2_graph)) + + +def test_CircuitGraph(): + # Test multiple parallel elements + circuit = "R0-p(C1,R1,R2)" + + assert len(CircuitGraph(circuit).execution_order) == 6 + + # Test nested parallel groups + circuit = "R0-p(p(R1, C1)-R2, C2)" + + assert len(CircuitGraph(circuit).execution_order) == 9 + + # Test parallel elements at beginning and end + circuit = "p(C1,R1)-p(C2,R2)" + + assert len(CircuitGraph(circuit).execution_order) == 7 + + # Test single element circuit + circuit = "R1" + + assert len(CircuitGraph(circuit).execution_order) == 1 + + # Test complex circuit with d and p + circuit = "d(p(R1,C1),p(R2,C2))" + params = [0.1, 0.01, 1, 1000] + frequencies = [1000.0, 5.0, 0.01] + + cg = CircuitGraph(circuit) + + # test for execution_order + assert len(cg.execution_order) == 7 + + # test for compute + assert len(cg.compute(frequencies, *params)) == len(frequencies) + + # test for calculate_circuit_length + assert cg.calculate_circuit_length() == 4 + + # test for __call__ + assert np.allclose(cg.compute_long( + frequencies, *params), cg(frequencies, *params)) + + # test for __eq__ + assert not cg == 1 + cg2 = CircuitGraph(circuit) + assert cg == cg2 + + # test for graph visualization + f, ax = plt.subplots() + cg.visualize_graph(ax=ax) + plt.close(f) + + # test for constants inputs + circuit = "d(p(R1,C1),p(R2,C2))" + params1 = [0.1, 0.01, 1] + constants = {'C2': 1000} + + cg = CircuitGraph(circuit, constants) + + assert np.allclose(cg(frequencies, *params1), cg2(frequencies, *params)) diff --git a/nleis/nleis_tests/test_nleis.py b/nleis/nleis_tests/test_nleis.py index 15186d0..da164cc 100644 --- a/nleis/nleis_tests/test_nleis.py +++ b/nleis/nleis_tests/test_nleis.py @@ -190,49 +190,6 @@ def test_eq(): simul_circuit == NLEIS_circuit -def test_NLEISCustomCircuit(): - circ_str = 'd(TDSn0,TDSn1)' - initial_guess = [ - 5e-3, 1e-3, 10, 1e-2, 100, 10, 0.1, - # TDS0 + additioal nonlinear parameters - 1e-3, 1e-3, 1e-3, 1e-2, 1000, 0, 0, - # TDS1 + additioal nonlinear parameters - ] - - NLEIS_circuit = NLEISCustomCircuit( - circ_str, initial_guess=initial_guess) - - assert not NLEIS_circuit._is_fit() - - # check get_param_names() - full_names_NLEIS, all_units_NLEIS = NLEIS_circuit.get_param_names() - assert full_names_NLEIS == ['TDSn0_0', 'TDSn0_1', 'TDSn0_2', 'TDSn0_3', - 'TDSn0_4', 'TDSn0_5', - 'TDSn0_6', 'TDSn1_0', 'TDSn1_1', 'TDSn1_2', - 'TDSn1_3', 'TDSn1_4', 'TDSn1_5', 'TDSn1_6'] - assert all_units_NLEIS == ['Ohms', 'Ohms', 'F', 'Ohms', 's', - '1/V', '-', 'Ohms', 'Ohms', 'F', 'Ohms', 's', - '1/V', '-'] - - # check _is_fit() - assert not NLEIS_circuit._is_fit() - - # check complex frequencies raise TypeError - with pytest.raises(TypeError): - NLEIS_circuit.predict([0.42, 42 + 42j]) - - # test is_fit method - NLEIS_circuit.fit(f, Z2) - assert NLEIS_circuit._is_fit() - - # test extract method - dict = NLEIS_circuit.extract() - assert list(dict.keys()) == ['TDSn0_0', 'TDSn0_1', 'TDSn0_2', 'TDSn0_3', - 'TDSn0_4', 'TDSn0_5', - 'TDSn0_6', 'TDSn1_0', 'TDSn1_1', 'TDSn1_2', - 'TDSn1_3', 'TDSn1_4', 'TDSn1_5', 'TDSn1_6'] - - def test_EISandNLEIS_fitting(): circ_str_1 = 'L0-R0-TDS0-TDS1' circ_str_2 = 'd(TDSn0,TDSn1)' @@ -353,8 +310,8 @@ def test_NLEISCustomCircuit(): 'TDSn0_6', 'TDSn1_0', 'TDSn1_1', 'TDSn1_2', 'TDSn1_3', 'TDSn1_4', 'TDSn1_5', 'TDSn1_6'] assert all_units_NLEIS == ['Ohms', 'Ohms', 'F', 'Ohms', 's', - '1/V', '', 'Ohms', 'Ohms', 'F', 'Ohms', 's', - '1/V', ''] + '1/V', '-', 'Ohms', 'Ohms', 'F', 'Ohms', 's', + '1/V', '-'] # check _is_fit() assert not NLEIS_circuit._is_fit() @@ -388,7 +345,7 @@ def test_NLEISCustomCircuit(): None, f, Z2, kind='bode')[0], type(ax)) # check altair plotting with a fit circuit - chart = NLEIS_circuit.plot(f_data=f, Z2_data=Z2, + chart = NLEIS_circuit.plot(f_data=f, Z2_data=Z2, max_f=10, kind='altair') datasets = json.loads(chart.to_json())['datasets'] From 2457aa165771b151431c71d78022af3fe7744d3a Mon Sep 17 00:00:00 2001 From: Yuefan Ji Date: Tue, 14 Jan 2025 10:33:39 -0800 Subject: [PATCH 5/6] Documentation and docstring improvements: * Added speed benchmark for graph-based function evaluation to examples * Added release notes * Improved CircuitGraph class docstrings --- docs/source/conf.py | 2 + docs/source/examples.rst | 3 +- docs/source/examples/graph_example.ipynb | 632 +++++++++++++++++++++++ docs/source/index.rst | 3 +- docs/source/release-notes.rst | 45 ++ nleis/fitting.py | 33 ++ nleis/nleis_elements_pair.py | 3 + 7 files changed, 719 insertions(+), 2 deletions(-) create mode 100644 docs/source/examples/graph_example.ipynb create mode 100644 docs/source/release-notes.rst diff --git a/docs/source/conf.py b/docs/source/conf.py index bd07311..f469097 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -95,6 +95,8 @@ # a list of builtin themes. # html_theme = 'sphinx_rtd_theme' +# Todo: change to pydata_sphinx_theme +# html_theme = 'pydata_sphinx_theme' # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the diff --git a/docs/source/examples.rst b/docs/source/examples.rst index 98cf9b2..737d601 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -6,4 +6,5 @@ Examples :maxdepth: 1 :glob: - examples/nleis_example \ No newline at end of file + examples/nleis_example + examples/graph_example \ No newline at end of file diff --git a/docs/source/examples/graph_example.ipynb b/docs/source/examples/graph_example.ipynb new file mode 100644 index 0000000..e68d6e7 --- /dev/null +++ b/docs/source/examples/graph_example.ipynb @@ -0,0 +1,632 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Graph Based Evaluation and Benchmark**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The original function evaluation in `impedance.py` and `nleis.py` relies on string parsing and the use of `eval()`, which can become slow for complex circuits and large inputs of parameters and frequencies. To address this, we have introduced a graph-based function evaluation that improves performance significantly, achieving at least a 3X speedup in certain applications.\n", + "\n", + "The following example demonstrates the implementation and provides benchmark results comparing the graph-based approach to the eval()-based evaluation." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "from nleis import NLEISCustomCircuit, EISandNLEIS\n", + "import timeit\n", + "import numpy as np\n", + "import warnings\n", + "import matplotlib.pyplot as plt\n", + "from tabulate import tabulate\n", + "\n", + "warnings.filterwarnings(\"ignore\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Enable Graph-based Calculation**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The graph-based evalution has been seamleassly built into `EISandNELIS` and `NLEISCustomCircuit`. The user can easily enable it by set `graph = True` in the circuit initialization. Beside this, all other method remains the same and can be used as usual." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Example for NLEISCustomCircuit" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "circ_str = 'd(TDSn0,TDSn1)'\n", + "initial_guess = [\n", + " 5e-3, 1e-3, 10, 1e-2, 100, 10, 0.1,\n", + " # TDS0 + additioal nonlinear parameters\n", + " 1e-3, 1e-3, 1e-3, 1e-2, 1000, 1, 0.01,\n", + " # TDS1 + additioal nonlinear parameters\n", + " ]\n", + "circuit = NLEISCustomCircuit(circ_str, initial_guess=initial_guess,graph=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once the circuit has been initialized, it is possible to retreieve the graph by calling `circuit.cg`. Once nice thing about the graph is that it is visualizable, which can be achieved using the `.visualize_graph` method. " + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAADTCAYAAADNnRQhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAriUlEQVR4nO3de3Bc1Z0n8G+/pW633k/LLbUe/bjXmOBQLsCBEPPyGoLBEGIcA6FwPFl2YAonpJbATCDMhJAKu6ECk0yyQx6Ay6YygcTLTBLIRmE9ITiwhBnjllqtR6v1flpSq1vq1737h6KLxb0ysi35qlvfTxVlqe/t2+cIkL91zvmdY5BlWQYRERER6cKodwOIiIiI1jKGMSIiIiIdMYwRERER6YhhjIiIiEhHDGNEREREOmIYIyIiItIRwxgRERGRjsxLuUmSJPT398PpdMJgMKx0m4iIiIiynizLiEajWL9+PYzGxce/lhTG+vv74XK5lq1xRERERGtFT08PNmzYsOj1JYUxp9OpPKygoGB5WkZERESUw6ampuByuZQctZglhbH5qcmCggKGMSIiIqIz8FFLvLiAn4iIiEhHSxoZyxWxRBrhsRiSaQlWsxHuUgcctjX1IyAiIqJVJueTSGgoioPHImgODiMyHod8yjUDgNoSO7b5KrD3klp4Kk8/p0tERES03AyyLMsfddPU1BQKCwsxOTmZNWvGesbjePiV4zjaPgqT0YCMtHg3569f0VSGJ3ZtgqvEfh5bSkRERLloqfkpJ9eMHX47gmu+8wbe7BwDgNMGsVOvv9k5hmu+8wYOvx1Z8TYSERERATk4TflscwhPvdZ2Vu/NSDIykoyHXj6O0ekE7tvmWebWERERES2UUyNjh9+OnHUQ+7CnXmvDSxwhIyIiohW2akbGlnrMUnNzM9xuN+rr65XXzGYznAWFmLFXwOq6AM6LdsBcWKF6b3piCBN/OIREz/tIR8dgzHPAUlyDvLoLUXTFXtX9XztyAlsby067huwb3/gGjh07hmPHjmF4eBiPPvooHnvssSX1hYiIiGjVhLEXXnhhwffPP/88Xn/9ddXrgiBgZmYGALBnzx5cf/31kCQJz/zqPbQe/zOibx9B9J0jKN3xN3CIVyrvS53sx+BPDsBgscFx4bUwF1YgMz2O5GAHJt/6F80wlpZkPPzKcbyw75JF2/23f/u3qKqqwubNm/Gb3/zmXH4EREREtAatmjB2xx13LPj+rbfewuuvv656HQDC4TAA4OMf/zjuuOMOhIai+FpLKUrdV6Pwk8MYeunvMPqv34Gl1AVrZQMAYOrtX0JKzaLmnmdUo2aZ2IRmmzKSjKPto2gfjqKpQnvbi66uLrjdboyOjqK8vPwMe01ERERrXU6sGTt4LAKTcW6a01xYgbIbHgAyaUwe+7lyT/rkAEzOMs3pS5OjaMH3vd+7B8M/+zpme05g8Kdfgn9DGRoaGvD888+r3ut2u5ezK0RERLTG5EQYaw4OL9i+wlYjwFxUjdnwe8pr5sIKZKZGMBP+jyU9M3VyACO/+CZs7otQf8N/RXFxMe6++26cOHFiuZtPREREa1jWh7HpRBqR8bjqdUt5HaT4JKTE3DXnxTfCYLJg+PAj6P/R32D8tz9EvO0tSKlZzeemx3tRfvNDKL7yLqT91+Hnv3wVVqsVP/7xj1e0P0RERLS2rJo1Y2ereywGrS1djdY8AICUjMNos8NaXofqe76LyT8cxkz7nxB9pxPRd47AYM1H8VX74Lzovyx4v6WsFnmuCwAAMoCY0Q6fz4fOzs4V7hERERGtJVkfxpJpSfN1KTk34mW0frAthaWkBmU3fhmylEFqtAczHX/C1Fs/x/ivn4W5qAr57ouUe00FCxfjJ9MSiouLcfLkyeXvBBEREa1ZWT9NaTVrdyE10g2jvQhGm3qPMIPRBGuFG4WXfRbltzwCAIid+P3CewwLn/vkN/4BPT09GBsbQzgchiRph0AiIiKiM5H1I2PuUgcMwIKpykRfC9ITA3Bs3PaR77dWzx15lJkeX/wmWUboz28qIay+vh52ux1+vx+iKEIQBGzYsAEAGNKIiIjojGR9GHPYzKgtsaP7L4v405PDGP3XpwGTGQWX3KLcN9vzPmzr/TCYFnZ5puMdAHNTmIupK3PgjT/9EZ/61KcwOzuLRx99FIFAAC0tLWhpacGrr76KiYkJAHM78r/88stKSBMEAaIowuv1Ii8vb3k7T0RERFkvq8PYu+++ixdffBElfRGcCEQw29+GePBNwACUffrLsFZ8cGTS1Fs/R3KwHfm+rbCWuwEAyaEOxN7/HYx5Tji33KT5GSajAdu8H+xNlpeXhx07dmDHjh0A5k4OCIfDGBkZwTPPPIOmpiY4HA688847aG5uxujoKADAaDSioaFBCWfzf/r9fjid2hvKEhERUe7L6jB26NAhHDp0CGazGZI5H+aS9XBu2al5NmXBZZ9FPPB7zPa8j9iJ30NOJWBaVwy78EkUfuJ2WIqqND8jI8m449LaRdvw3HPP4Y033lC+b2v74KDy5uZmXHjhhcoI2vxo2qFDhxCJfHAIucvlWhDS5r8uLS092x8NERERZQmDLMtaO0MsMDU1hcLCQkxOTqKgoOB8tOuM3fncMbzZObZg89dzZTIasLWh9LRnU56t6elpBINBJaDN/9ne3q6sOysvL18wijb/Z3V19ZIPViciIiJ9LDU/5UwY6xmP45rvvIHEIltdnA2b2YjfHrgSrhJ1ReZKSSQSCIVCqpAWDAaRTCYBAAUFBZohra6uDkZj1hfIEhER5YQ1F8YA4PDbETz08vFle963btmE3VsWn6I8n9LpNLq6uhYUDsx/HYvFAAD5+fnw+XyqoNbU1ASLxaJzD4iIiNaWNRnGAODZ5hCeeq3to2/8CF+5zoe/3ta0DC1aWZIkobe3d0E4CwQCCAQCyga1ZrMZHo9HFdJ8Ph/y8/N17gEREVFuWrNhDJgbIXv0yAmkJfmM1pCZjAaYjQY8vnPjqhkRO1uyLGN4eFgV0lpaWjAwMAAAMBgMqK+vVxUOCIKQFf+eiYiIVrM1HcaAuTVkD79yHEfbR2EyGk4byuavX9FUhid2bTqva8T0MDExoZrqDAQCCIfDyj01NTWqNWmCIKC8vHzxBxMREZFizYexeaGhKA4ei6C5bRiRsfiCnfoNAGpL7djmrcAdl9aiqWJt7/cVi8UQDAZVIa29vR2ZTAYAUFZWphnSampqWOFJRER0CoYxDbFEGuGxGJJpCVazEe5SBxy2rN5q7bxIJpMIhUKq0bTW1lYkEgkAgNPpVE11iqIIt9sNk8mkcw+IiIjOP4YxWnGZTAbhcFi1Ji0QCGB6ehrA3IkFPp9PFdKamppgtVp17gEREdHKYRgj3ciyjL6+Ps2QNjY2BgAwmUzweDyqkwf8fj/s9txes0dERGsDwxitSiMjI5oVnn19fQDmKjzr6upUa9IEQUBRUZG+jSciIjoDDGOUVSYnJ9Ha2qoKaV1dXZj/T7S6ulrz5IHy8nIWDxAR0arDMEY5YWZmRrPCMxQKIZ1OAwBKSko0KzxdLhdDGhER6YZhjHJaKpVCe3u7KqS1trZidnYWALBu3Tr4/X5VSGtoaGCFJxERrTiGMVqTMpkMuru7NdelTU1NAQBsNhu8Xq8qpHk8HthsNp17QEREuYJhjOgUsixjYGBAs8JzZGQEwFyFZ2Njo2rK0+/3w+Fw6NwDIiLKNgxjREs0OjqqeTxUb2+vck9dXZ3murTi4mIdW05ERKsZwxjROZqamkJra6sqqHV2dkKSJABAZWWlZoVnZWUliweIiNY4hjGiFTI7O4u2tjbVlGdbWxtSqRQAoKioSDOkuVwuGI1GnXtARETnA8MY0XmWSqXQ2dmphLNTKzzj8TgAwG63K5vYnhrUGhsbYTbznFQiolzCMEa0SkiShEgkolnhOTExAQCwWq3weDyq0TSv14u8vDx9O0BERGeFYYxolZNlGYODg5ohbWhoCABgNBrR0NCw4Fio+QpPp9Opcw+IiOh0GMaIstj4+LhmhWckElHucblcmhWepaWlOraciIjmMYwR5aDp6WmlwvPUkNbR0aFUeFZUVGiGtOrqalZ4EhGdRwxjRGtIIpFAKBRSTXcGg0Ekk0kAQGFhoapwQBRF1NXVscKTiGgFMIwREdLpNLq6ulQhraWlBbFYDACQn58Pv9+vCmmNjY2wWCw694CIKHsxjBHRoiRJQm9vr+bxUCdPngQAWCwWeDyeBSFNEAT4fD7k5+fr3AMiotWPYYyIzpgsyxgeHtas8BwYGAAAGAwG1NfXq9akCYLA3w9ERKdgGCOiZTUxMaEZ0sLhsHJPTU2NZkgrLy/Xr+FERDphGCOi8yIWiyEYDKqCWnt7OzKZDACgrKxMs8KzpqaGFZ5ElLMYxohIV8lkEqFQSBXSgsEgEokEAMDpdGqGNLfbDZPJpHMPiIjODcMYEa1KmUwGXV1dqurOQCCA6elpAEBeXh58Pp8qqDU1NcFqtercAyKipWEYI6KsIssy+vr6NCs8x8bGAABmsxlNTU2qkObz+WC323XuARHRQgxjRJQzRkZGNENaf38/gLkKT7fbrTnlWVhYqHPriWitYhgjopw3OTmpTHOeGtS6urow/6uturpaFdBEUUR5eTmLB4hoRTGMEdGaFY/H0dbWphpNC4VCSKfTAICSkhLNkLZhwwaGNCJaFgxjREQfkkql0N7eriocaG1txezsLABg3bp1yv5opwa1hoYGVngS0RlhGCMiWqJMJoPu7m7NTW2npqYAADabDV6vVzWa5vF4YLPZdO4BEa1GDGNEROdIlmUMDAxoFg+MjIwAAEwmExobG1Uhze/3w+Fw6NwDItITwxgR0QoaHR1VFQ4EAgH09vYq99TV1WlWeBYXF+vYciI6XxjGiIh0MDU1hdbWVlVQ6+zshCRJAICqqirNkFZZWcniAaIcwjBGRLSKzM7OalZ4trW1IZVKAQCKi4tVhQOiKMLlcsFoNOrcAyI6UwxjRERZIJVKobOzU7PCMx6PAwDsdrtmSGtoaIDZbNa5B0S0GIYxIqIsJkkSIpGIZoXnxMQEAMBqtcLr9aqCmtfrRV5enr4dICKGMSKiXCTLMgYHBzVD2tDQEADAaDSioaFBs8LT6XTq3AOitYNhjIhojRkfH9es8IxEIso9LpdLFdIEQUBpaamOLSfKTQxjREQEAJienlYqPE8NaR0dHUqFZ0VFhWaFZ3V1NSs8ic4SwxgREZ1WIpFAKBRSTXcGg0Ekk0kAQGFhoWZIq6urY4Un0UdgGCMiorOSTqfR1dWlCmktLS2IxWIAgPz8fPj9flVQa2xshMVi0bkHRKsDwxgRES0rSZLQ29urmu4MBAI4efIkAMBiscDj8ahCmtfrRX5+vs49IDq/GMaIiOi8kGUZw8PDmhWeAwMDAACDwYCGhgbVNhyCIPDvFcpZDGNERKS7iYkJzQrPcDis3FNTU6NakyaKIsrKyvRrONEyYBgjIqJVKxaLIRgMqkJae3s7MpkMAKCsrEwzpK1fv54VnpQVGMaIiCjrJJNJhEIh1Whaa2srEokEAMDpdGpWeLrdbphMJp17QPQBhjEiIsoZmUwG4XBYs8IzGo0CAPLy8uDz+VQhrampCVarVece0FrEMEZERDlPlmX09fVpVniOjY0BAMxmM5qamlSjaT6fD3a7XeceUC5jGCMiojVtZGREs8Kzr68PwFyFp9vt1pzyLCws1Ln1lAsYxoiIiDRMTk5qHg/V1dWF+b8S169fr9qGQxRFlJeXs3iAloxhjIiI6AzMzMxoVniGQiGk02kAQElJiWaF54YNGxjSSIVhjIiIaBmkUim0t7drVnjOzMwAANatW6e5oW1DQwMrPNcwhjEiIqIVJEkSuru7VWvSAoEApqamAAA2mw1er1c1mubxeGCz2XTuweoSS6QRHoshmZZgNRvhLnXAYTPr3axzwjBGRESkA1mWMTAwoFk8MDw8DAAwmUxobGxUhTS/3w+Hw6FzD86f0FAUB49F0BwcRmQ8jlMDiQFAbYkd23wV2HtJLTyVTr2aedYYxoiIiFaZsbExzZDW09Oj3FNXV6dZ4VlcXKxjy5dXz3gcD79yHEfbR2EyGpCRFo8i89evaCrDE7s2wVWSPduRMIwRERFliWg0qlnh2dnZCUmSAABVVVWaIa2ysjKrigcOvx3Bo0dOIC3Jpw1hH2YyGmA2GvD1nRtx+5baFWzh8mEYIyIiynKzs7Noa2tThbS2tjakUikAQHFxsapwQBRFuFwuGI1GnXuw0LPNITz1Wts5P+fB67y4b5tnGVq0shjGiIiIclQ6nUZHR4fqaKiWlhbE43EAgN1u19wrraGhAWbz+V8Yf/jtCB56+fii1yeOHsTkHw6h7qFXl/S8b92yCbtX+QjZUvNTdpcpEBERrUFmsxk+nw8+nw8333yz8rokSejp6VGtSXv11VcxMTEBALBarfB6vaqQ5vV6z7rCs7OzEw6HA5WVlZrXe8bjePTIiTN65kzXu4i1HEWyP4jUWC9MzjJs+G8/Uq5/7cgJbG0sy6o1ZIthGCMiIsoRRqMRdXV1qKurw44dO5TXZVnG0NCQKqT94Ac/wODgoPLehoYGzQpPp/P0lYzXX389+vv78eMf/xi33nqr6vrDrxxH+gzWhwFA7MQbiLcehbWyEaZ1JarraUnGw68cxwv7Ljmj565GDGNEREQ5zmAwoKqqClVVVbjqqqsWXDt58qRqTdrhw4fR3d2t3ONyuVQhTRAElJaWIpFIIBQKQZIkfOYzn8EXvvAFPP3008oWHaGhKI62j55xm4uuvAulO+6HwWTG8M++juRI94LrGUnG0fZRtA9H0VSRfdtenIphjIiIaA0rLi7G1q1bsXXr1gWvT09PIxgMLghp//Zv/4bvfve7SoVnRUUFamtrle8B4Ec/+hGam5vxs5/9DJs3b8bBYxHV9hWzPSdw8v/8M5IjYZidpSi4RD2aZnaWfmTbTUYDXnwrgsd2bjzb7q8KDGNERESksm7dOlx88cW4+OKLF7w+PxI2H9B+85vfLLguSRI6Ojpw8cUX4+///u/RbN26IIglh8MYfulrMNoLUHT55yBLGUz8+0GY7EVn3MaMJKO5bRiPgWGMiIiI1gibzYYLLrgAF1xwgfLan/70J2QyGQBza88kSYLdbgcseYiMxxe8f+LoiwBkVO39FsyFFQAAh+8T6H/ur8+qPZGxOGKJdFYfnbS6NiAhIiKirNLS0qIEsc2bN+Pxxx/H+++/j2g0ipvv+MKCI45kKYPZrj8j33OpEsQAwFLmQn7Dx8/q82UA4bHYOfRAfwxjREREdNai0SgA4GMf+xhuv/123Hbbbdi4cSMMBgOSaWnBvVJ8CnI6AUvxetVzzCU1Z92GD39OtmEYIyIiorPm9/thMBjwn//5n/jqV78Kn88Ht9uN3bt3oyfcdV7aYDVnd5zJ3glWIiIi0k0ymUQoFEI8Hsf8YT7zf3Z3d6O7uxvTiRQMwj5lqtJoL4DBbEPqZL/qeenxvrNqhwGAu9RxVu9dLRjGiIiIaFHxeFzzEPP29nZlrdip5g8t/9KXvoRvfvObuObpf0f3XxbxG4wm5NVvxkzoLaQnh5V1Y6nRHsx0vntW7asttWf14n2AYYyIiIgATExMqAJXS0sLwuGwck9NTQ1EUcT27dtx4MABCIKA6upqeL1eAIDJZEJVVRVeeuklfOITnwAAbPNV4IVj3cr2FkVX7MVg17sYPPjf4fz4DYCUwdT/+9+wlNUiNfLBZyWHuxAPHQMApE4OQE7EMPGHwwAAa0U97J5LYDIasM37QSFAtmIYIyIiWiNkWcbw8LBm6BoYGAAwN7JVX18PURRx2223LTgWqbCwUPO5ZWVlGB0dxS233IIf/vCHKCoqUq7tvaQWP/ljWPneWlGPis8+jpO/+2dMHH0RZmcZii7fi8z0OCZPDWODHZg8+uKCz5n/3nHB1bB7LkFGknHHpav7sPClMMjzE7ynsdRTx4mIiEh/sixrHhje0tKC8fFxAHOHjc8fGH7qMUc+nw/5+fln9Hk//elPYbFYsGfPHmWa8lR3PncMb3aOLdj89VyZjAZsbShd1WdTLjU/MYwRERFlqXQ6ja6uLs3QFYvN7b2Vn58Pv9+/IHCJoojGxkZYLJbz0s6e8Tiu+c4bSCzjFhQ2sxG/PXAlXCX2ZXvmcltqfuI0JRER0SqXSCTQ1tamBK350BUMBpFMJgEABQUFEEURmzZtwu7du5XQVVdXB6NR360fXCV2fH3nRjz08vFle+bjOzeu6iB2JhjGiIiIVonp6WnNysXOzk6lcrGiogKCIODyyy/H/v37ldGu6upqzSnC1eL2LbUYnU7gqdfazvlZX7nOh91bsn+t2DyGMSIiovNsfHxcNcoVCAQQiUSUe1wuFwRBwA033LBgirG0tFTHlp+b+7Z5ULbOhkePnEBaks9oDZnJaIDZaMDjOzfmVBADuGaMiIhoRciyjMHBQc3KxaGhIQBzh2o3NDQsWMs1X7nodDp17sHK6RmP4+FXjuNo+yhMRsNpQ9n89SuayvDErk1ZNTXJBfxERETngSRJiEQimqFrYmICAGCxWOD1elWhy+v1Ii8vT98O6Cg0FMXBYxE0tw0jMhZfcKi4AXMbum7zVuCOS2vRVJF94ZRhjIiIaBml02l0dHSoAldrayvi8bkd5u12u7JVxKmhq7GxEWYzVwadTiyRRngshmRagtVshLvUkfU767OakoiI6CzMzs6ira1NFbra2tqQSqUAAEVFRRBFEZs3b8bevXuV4OVyuXSvXMxWDpsZG9drbyqb6xjGiIhoTYpGo8oi+lNDV2dnJyRpbj+sqqoqCIKAK6+8Evfee68SuiorK1d15SJlF4YxIiLKaWNjY6pRrkAggN7eXuWeuro6CIKAnTt3LpheLC4u1rHltFYwjBERUdaTZRn9/f2a20WMjIwAmKtcbGpqgiAIuPPOO5XQ5fP5sG7dOp17QGsZwxgREWUNSZIQDoc1KxenpqYAAFarFT6fD6IoYtu2bcool8fjgc1m07kHRGoMY0REtOqkUim0t7erQlcwGMTMzAwAwOFwKKNbN998sxK66uvrWblIWYX/tRIRkW5mZmYQDAZVo1yhUAjpdBoAUFJSAlEUsWXLFtx1111K6HK5XFxETzmBYYyIiFbc1NSU5tRiV1cX5re7rK6uhiiKuPrqq3H//fcro17l5eUMXZTTGMaIiGjZjIyMaIauvr4+5R632w1RFLFr1y5llEsQBBQVFenXcCIdMYwREdEZkWUZfX19mttFjI2NAQBMJhM8Hg8EQcDdd9+9oHLRbs+eswWJzgeGMSIi0pTJZBAOh1Whq6WlBdFoFABgs9ng9/shCAKuvfZaJXQ1NTXBarXq3AOi7MAwRkS0xiWTSYRCIc3KxUQiAQBwOp1K0Lr11luVr91uN0wmk849IMpuDGNERGtEPB5Ha2urKnS1t7cjk8kAAEpLSyGKIi677DLcc889ypqumpoaLqInWiEMY0REOWZiYkI1rRgIBNDd3a1ULtbU1EAURWzfvh0PPPCAErrKy8t1bj3R2sMwRkSUhWRZxvDwsGbl4sDAAADAYDCgvr4eoijitttuUwKX3+9HYWGhzj0gonkMY0REq5gsy+jp6VGNcrW0tGB8fBwAYDab4fF4IIoi9u3bp4Qun8+H/Px8nXtARB+FYYyIaBVIp9Po6urSrFyMxWIAgPz8fKVycceOHcoi+sbGRlgsFp17QERni2GMiOg8SiQSCIVCqtAVDAaRTCYBAAUFBRBFEZs2bcLu3buV0FVXVwej0ahzD4houTGMERGtgOnpac3Kxc7OTqVysaKiAoIg4PLLL8f+/fuV6cXq6mpWLhKtIQxjRETnYHx8XLWWKxAIIBKJKPe4XC4IgoAbbrhBGeUSBAGlpaU6tpyIVguGMSKijyDLMgYHBzW3ixgaGgIAGI1GNDQ0QBAE7NmzRwldfr8fTqdT5x4Q0WrGMEZE9BeSJCESiWhuFzExMQEAsFgs8Hq9EEURX/ziF5VRLq/Xi7y8PH07QERZiWGMiNacdDqNjo4OVehqbW1FPB4HANjtdvj9foiiiE9/+tNK6GpsbITZzF+dRLR8+BuFiHLW7Ows2traVKNcbW1tSKVSAICioiKIoojNmzdj7969EAQBgiCgtraWlYtEdF4wjBFR1otGo2htbVWFrs7OTkiSBACorKyEKIq48sorce+99ypruiorK1m5SES6YhgjoqwxNjamuSlqT0+Pck9dXR0EQcDOnTuVUS5BEFBSUqJjy4mIFscwRkSriizLGBgYUIWuQCCAkZERAHOVi01NTRAEAXv37oUoihBFET6fD+vWrdO5B0REZ4ZhjIh0IUkSwuGw5nYRU1NTAACr1QqfzwdRFLFt2zZlatHj8cBms+ncAyKi5cEwRkQrKpVKob29XTXKFQwGMTMzAwBwOBxK0LrpppuUysX6+npWLhJRzuNvOSJaFjMzMwgGg6pRrlAohHQ6DQAoKSmBIAjYsmUL7rrrLiV0bdiwgZWLRLRmMYwR0RmZmprS3BS1q6sLsiwDAKqrqyGKIq6++mrcf//9yqhXeXk5KxeJiD6EYYyINI2MjGiGrr6+PuUet9sNURSxa9cuZZRLEAQUFRXp13AioizDMEa0hsmyjL6+Ps3KxbGxMQCAyWSCx+OBIAi4++67lVEun88Hu92ucw+IiLIfwxjRGpDJZNDV1bVgLdf819FoFABgs9ng9/shCAKuvfZaJXQ1NTXBarXq3AMiotzFMEaUQ5LJJEKhkGblYiKRAAA4nU4laN16663K1263GyaTSeceEBGtPQxjRFkoHo+jtbVVFbra29uRyWQAAGVlZRAEAZdddhnuueceZU1XTU0NF9ETEa0iDGNEq9jExITmIvpwOKzcU1NTA1EUsX37dhw4cEBZRF9eXq5fw4mIaMkYxoh0JssyhoeHNUPXwMAAAMBgMKC+vh6iKOK2225TRrn8fj8KCwt17gEREZ0LhjGi80SWZfT09GhWLp48eRIAYDab4fV6IQgC9u3bp4Qun8+H/Px8nXtAREQrgWGMaJml02l0dXWpQldLSwtisRgAID8/X6lc3LFjh7KIvrGxERaLReceEBHR+cQwRnSWEokE2traVIErGAwimUwCAAoLCyEIAi688ELs3r1bCV11dXU8/oeIiAAwjBF9pOnpac3KxY6ODkiSBACoqKiAIAi4/PLLsX//fmV6sbq6mpWLRER0WgxjRH8xPj6uGuUKBAKIRCLKPS6XC4Ig4IYbblBGuQRBQGlpqY4tJyKibMYwRmuKLMsYHBzUrFwcGhoCABiNRjQ0NEAURezZs2dB5aLT6dS5B0RElGsYxignSZKESCSiGuVqaWnBxMQEAMBiscDr9UIURXzyk59UQpfX60VeXp6+HSAiojWDYYyyWjqdRkdHh2qUq7W1FfF4HABgt9uVjVBvvPFGZXqxoaEBZjP/FyAiIn3xbyLKCrOzswgGg6o1XW1tbUilUgCAoqIiiKKIzZs3Y+/evUrocrlcrFwkIqJVi2GMVpVoNKqaVgwEAujq6lIqF6uqqiAIAj71qU/h3nvvVUJXZWUlKxeJiCjrMIyRLkZHRzVDV29vr3JPXV0dBEHATTfdtKBysbi4WMeWExERLS+GMVoxsiyjv79fs3JxZGQEAGAymdDY2AhRFHHnnXcuOP5n3bp1OveAiIho5TGM0TmTJAnhcFgzdE1NTQEArFYrfD4fRFHEVVddpYQuj8cDm82mcw+IiIj0wzBGS5ZKpdDe3q5ZuTg7OwsAcDgcStDatWuXUsVYX1/PykUiIiINa+pvx1gijfBYDMm0BKvZCHepAw7bmvoRLMnMzAyCwaAqdIVCIaTTaQBASUkJRFHEli1b8PnPf15Z07VhwwYuoiciIjoDOZ9EQkNRHDwWQXNwGJHxOORTrhkA1JbYsc1Xgb2X1MJTubZ2V5+cnFQtom9paUFXVxdkee4ntX79egiCgGuuuQb333+/ErrKy8sZuoiIiJaBQZ7/W/c0pqamUFhYiMnJSRQUFJyPdp2znvE4Hn7lOI62j8JkNCAjLd7N+etXNJXhiV2b4Cqxn8eWrryRkRHVKFcgEEB/fz8AwGAwwO12L6hYFEURfr8fRUVF+jaeiIgoSy01P+VkGDv8dgSPHjmBtCSfNoR9mMlogNlowNd3bsTtW2pXsIXLT5Zl9Pb2ah50PTY2BgAwm81oampShS6fzwe7PbcCKBERkd6Wmp9ybpry2eYQnnqt7azem/lLeHvo5eMYnU7gvm2eZW7ductkMujq6lKNcrW2tiIajQIA8vLylMrFa6+9VgleTU1NsFqtOveAiIiITpVTYezw25GzDmIf9tRrbShfZ8NunUbIkskkQqGQapQrGAwikUgAAJxOJwRBwMaNG/GZz3xGCV1utxsmk0mXdhMREdGZWTVhbKmLwZubm+F2u1FfX6+8Zjab4SwoxIy9AlbXBXBetAPmwgrVe9MTQ5j4wyEket5HOjoGY54DluIa5NVdiKIr9qru/9qRE9jaWHbaNWSSJOGpp57C97//fQwMDMDr9eKrX/0q9uzZs6T+xGIxzcrF9vZ2ZDIZAEBZWRlEUcRll12Gffv2KdOL69ev5yJ6IiKiLLdqwtgLL7yw4Pvnn38er7/+uup1QRAwMzMDANizZw+uv/56SJKEZ371HlqP/xnRt48g+s4RlO74GzjEK5X3pU72Y/AnB2Cw2OC48FqYCyuQmR5HcrADk2/9i2YYS0syHn7lOF7Yd8mi7X7kkUfw5JNPYv/+/diyZQt++ctf4nOf+xwMBgNuv/125b6JiQnNTVHD4bByT01NDURRxPbt23HgwAFlj67y8vIz+lkSERFR9li1C/jvu+8+/OM//iO0mhcOh1FfX49vf/vbePDBBxEaiuLap/8vACA9OYyhl/4O6ckhVN/1P2GtbAAAjL32fUy/92vUfPF/qUbNMrEJmBxFi7bltwc+iaYK9bYXfX19qK+vx1/91V/h2WefhSzLGBoawnXXXYeenh7s2bNHGfUaHBwEMDcC2NDQsGAR/fw/2VAcQUREREuzphbwHzwWUbanMBdWoOyGBzD4wlcweeznKN/5FQBA+uQATM4yzenLDwex3u/dA2t5HQou/Qwmfvcc/P8jjNoNNXjsscdw1113QZIk9PT04Nvf/jZSqRQGBwdx+eWXIxAI4OTJk8pzfvWrX+Hiiy/G/v37lfDl9XqRn5+/oj8PIiIiyh45Ecaag8MLtrCw1QgwF1VjNvye8pq5sAKz4fcwE/4P5Ls/9pHPTJ0cwMgvvol1F16His1XIfHeEXz+85/Ht771LXR3dyMWiyn3dnR0YOPGjdixYwdEUYTD4cD27dvxpS99Cffff/+y9pWIiIhyS9aHselEGpHxuOp1S3kdZkJvQUrEYbTZ4bz4RsTeb8bw4UdgqWhAXu0FyKu9EHn1F8FoyVO9Pz3ei8q9TyLPdQFkWUb0j/8Cg8EAo9GIxx57DKIo4qmnnkJ3dzf+/Oc/L3hvPD7XnvlNVYmIiIgWY9S7AeeqeywGrUVvRutcwJKSc8HIWl6H6nu+C8fGbchMDiH6zhGMvPwP6H3mTkTf+7Xq/ZayWuS5LgAwt87rj8fbsWnTJng8Hjz44IO4/vrrYTAYYLPZVO/Ny5v77PlCAyIiIqLFZP3IWDItab4uJWcBAEbrB9tSWEpqUHbjlyFLGaRGezDT8SdMvfVzjP/6WZiLqpDvvki511SwsIIxmZFQXFy8YE1Yfn6+sufXqWZnZ5XrRERERKeT9SNjVrN2F1Ij3TDai2C0qfcIMxhNsFa4UXjZZ1F+yyMAgNiJ3y+8x7DwufOfc2p1Z3V1NQYHB1UVnwMDAwDmDtkmIiIiOp2sD2PuUgc+vO1poq8F6YkB5Ndv/sj3W6vnjjzKTI8veo/hL5/zYRdddBHi8ThaWloWvH7s2DHlOhEREdHpZH0Yc9jMqD1lh/z05DBG//VpwGRGwSW3KK/P9rwPOZNWvX+m4x0Ac1OYi6kttcNhU8/o3nTTTbBYLPje976nvCbLMv7pn/4JNTU12Lp169l0iYiIiNaQrF4z9u677+LFF19ESV8EJwIRzPa3IR58EzAAZZ/+MqwVHxyZNPXWz5EcbEe+byus5W4AQHKoA7H3fwdjnhPOLTdpfobJaMA2r3pvMgDYsGEDHnjgAWW/sS1btuAXv/gFjh49ioMHD/J8SCIiIvpIWR3GDh06hEOHDsFsNkMy58Ncsh7OLTs1z6YsuOyziAd+j9me9xE78XvIqQRM64phFz6Jwk/cDktRleZnZCQZd1y6+GHhTz75JIqLi/GDH/wAP/nJT+DxePDiiy/ic5/73LL2lYiIiHLTqj0O6Uzd+dwxvNk5tmDz13NlMhqwtaH0tGdTEhEREWlZan7K+jVj857YtQlm44eX8p8bs9GAJ3ZtWtZnEhEREZ0qZ8KYq8SOr+/cuKzPfHznRrhK1FtjEBERES2XnAljAHD7llo8eJ13WZ71let82L1l8bViRERERMshqxfwa7lvmwdl62x49MgJpCX5jNaQmYwGmI0GPL5zI4MYERERnRc5NTI27/YttfjtgSuxtaEUwFzIOp3561sbSvHbA1cyiBEREdF5k3MjY/NcJXa8sO8ShIaiOHgsgua2YUTG4gsOFTdgbkPXbd4K3HFpLZoqnHo1l4iIiNaonNnaYiliiTTCYzEk0xKsZiPcpQ7NnfWJiIiIztVS89OaSiIOmxkb1xfq3QwiIiIiRU6uGSMiIiLKFksaGZufyZyamlrRxhARERHlivnc9FErwpYUxqLRKADA5XKdY7OIiIiI1pZoNIrCwsWXSS1pAb8kSejv74fT6YTBsLxHDhERERHlIlmWEY1GsX79ehiNi68MW1IYIyIiIqKVwQX8RERERDpiGCMiIiLSEcMYERERkY4YxoiIiIh0xDBGREREpCOGMSIiIiIdMYwRERER6ej/A5AD0aKABDiyAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "f, ax = plt.subplots(figsize=(6,2), layout=\"constrained\")\n", + "circuit.cg.visualize_graph(ax=ax)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Example for EISandNLEIS" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The graph initialization follows the same logic for `EISandNLEIS` class. The only difference is that graph are initialized separately for EIS and NLEIS circuit." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "circ_str_1 = 'L0-R0-TDS0-TDS1'\n", + "circ_str_2 = 'd(TDSn0,TDSn1)'\n", + "initial_guess = [1e-7, 1e-3, # L0,RO\n", + " 5e-3, 1e-3, 10, 1e-2, 100, 10, 0.1,\n", + " # TDS0 + additioal nonlinear parameters\n", + " 1e-3, 1e-3, 1e-3, 1e-2, 1000, 1, 0.01,\n", + " # TDS1 + additioal nonlinear parameters\n", + " ]\n", + "simul_circuit = EISandNLEIS(circ_str_1, circ_str_2,\n", + " initial_guess=initial_guess, graph=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "EIS and NLEIS graph can be extracted separately by calling `simul_circuit.cg1` and `simul_circuit.cg2` respectively" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAADTCAYAAADNnRQhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4bElEQVR4nO3de3hb5Z0n8K+ObpZkWb47jiNbvkmWQgKhhEAoUFMCZCihoaXABmhLykxnl3aadtil2enQdjp05hmm8HTb3afddso0sEC7wAxL6SzQMZSWEqBLS4hvUmxHduK7Y+tmS5bO2T+EDhbn2JFj2bLk7+d5/BAfHUnnOMH++ve+v/fVSJIkgYiIiIhyQsj1BRARERFtZAxjRERERDnEMEZERESUQwxjRERERDnEMEZERESUQwxjRERERDnEMEZERESUQ7pMThJFEadPn4bVaoVGo1ntayIiIiLKe5IkIRgMYvPmzRCExetfGYWx06dPw263Z+3iiIiIiDaKwcFBbNmyZdHHMwpjVqtVfrGSkpLsXBkRERFRAQsEArDb7XKOWkxGYSw1NFlSUsIwRkRERLQMZ5vixQn8RERERDmUUWWsUISjcQxMhhGLizDoBDgqLLAYN9SXgIiIiNaZgk8i3tEgHjvqR0fPGPxTEUgLHtMAqC83o91VjQO76tFas/SYLhEREVG2aSRJks52UiAQgM1mw8zMTN7MGRuciuDwM8fwqm8CWkGDhLj4baYev7ylEg/s3wZ7uXkNr5SIiIgKUab5qSDnjD3xph9XP/QKXuubBIAlg9jCx1/rm8TVD72CJ970r/o1EhEREQEFOEz5vQ4vHnyh95yemxAlJEQJ9z19DBOhKO5pb83y1RERERGlK6jK2BNv+s85iH3Qgy/04klWyIiIiGiVFUxlbHAqgvufPY7QOy9h8vmH339AI0BrKUWRYwdKr7wDOmtl2vPmJwYx9av/iehQJzRaHUzNO1H20c9Ba7bhr589jt3NlZxDRkRERKumYMLY4WeOIb5gbpjt8gPQ2TZBSsQQO9WD0LGXEB3qxObPfR8anQEAEA9MYOSx/wLBaEHplXdCis0h8MbTiI0PoPbT30FcMODwM8dw5OCuXN0WERERFbiCCGPe0SBe9U2kHTM1XQRj7Xtzvs6/FoK5BIHX/zci3qOwuC8HAMz87meQ5qOo+czD0NmqAQCGzU6MPfFXCB37FawXXIdXfRPwjQXRUs1lL4iIiCj7CmLO2GNH/dAKS281YNyyFQAQnx6Wj0V6XoOpZaccxADA5LgAuvI6RLpeBZBc9uLR1zl3jIiIiFZHQYSxjp6xsy9fMTMKABCKigEA8eAExMg0DJtaFOcaa52IjfYlnydK6Ogdy/IVExERESXl/TBlKBqHfyqiOC5Gw0hEZiDF5xE93YPp3z4OaPUwNV8MAEiEzgAAtMXliudqi8sgzgUhxeeh0enhn4wgHI1z6yQiIiLKurxPFycnw1CriY098Vdpn2ttNai+4SvQlSS7KaV4FACg0eoVz9VoDfI5Gp0eEoCByTC2brZl9dqJiIiI8j6MxeKi6vHya/4cuvI6iHNhhI+9iLnB42nBS6MzAgCkxLziuVIilnbOUu9DREREtBJ5H8YMOvVpb4Zap9xNaXZegpFH/zMmnv0HbP7TH0AwmKAtLgMAJEJTiucmQmcgFFmh0b0f3v7ub7+Fi52b4fF44Ha7UV9fD0EoiCl3RERElEN5H8YcFRZoANWhyhSNoEXZlZ/G6OOHEfz9c7BdejN01koIZhtiIz7F+dHhXhhqGt8/IEnwvv0annv0HUQiyflpZrMZbW1tcjhL/be5uRk6Xd5/WYmIiGiN5H1qsBh1qC8346TKJP6Fihq2w1DrROCtf0XJzhuh0Rlgdu1G+Ni/Ix4Yh66kCgAwO/AHxKdOoWTnjfJzGyoteOWN30EURfj9fnR1daGzsxNdXV3o6urCc889h+npaQCAXq+H0+mUw1kqqDmdThQVFa3a14GIiIjyU96HMQBod1XjyNGTZ13eomTXTZj4l79D6NhLsO74E9gu/RQi3b/F6P86DOtF+yDNzyJw9Gnoqxwo3rYHQHKdsXZnch0yQRDgcDjgcDiwd+9e+XUlScLo6Kgc0FL//cEPfoDR0VH5uU1NTWlVNI/Hg7a2NlitXFCWiIhooyqIMHZgVz0e+d3AWc8zu3ZDV1qLwNFnUHz+tdCVVKHmP3wbZ/79R5h+5RFoBB1MLTtRdtVBeb5YQpRw+yX1S76uRqPBpk2bsGnTJlx11VVpj01NTckVtFRIe/zxx+H3v7+QrN1uTwtpqT9XVFQs/4tBREREeUUjSdLS5SQAgUAANpsNMzMzKCkpWYvrWrY7fnwUr/VNnrU6thxaQYPdTRWrsjdlKBRCT0+Poprm8/kgisnOzaqqKsWcNI/Hg9raWmg0S+84QERERLmVaX4qmDA2OBXB1Q+9gmgWl6Aw6gS8dOhK2MvNWXvNs4lGo/B6vYqQ1tPTg1gsueRGSUmJakhraGhghycREdE6seHCGAA88aYf9z19LGuv9/c3bcMtO5ceolwr8Xgc/f39aY0DqT+Hw2EAgMlkgsvlUgS1lpYW6PXKxW2JiIho9WzIMAYA3+vw4sEXelf8Ovde48J/alfuW7neiKKIoaGhtHDW2dmJzs5OnDmT3PJJp9OhtbVVEdJcLhdMJlOO74CIiKgwbdgwBiQrZPc/exxxUVrWHDKtoIFO0OCb+7aum4rYuZIkCWNjY4qQ1tXVheHhYQDJxoPGxkZF44Db7c6Lv2ciIqL1bEOHMSA5h+zwM8fwqm8CWkGzZChLPX55SyUe2L9tTeeI5cL09LRiqLOzsxMDAwPyOXV1dYo5aW63G1VVVbm7cCIiojyy4cNYinc0iMeO+tHROwb/ZCRtpX4NgPoKM9qd1bj9knq0VG/s9b7C4TB6enoUIc3n8yGRSAAAKisrVUNaXV0dOzyJiIgWYBhTEY7GMTAZRiwuwqAT4KiwwGIsiKXWVlUsFoPX61VU07q7uxGNRgEAVqtVMdTp8XjgcDig1WpzfAdERERrj2GMVl0ikcDAwIBiTlpnZydCoRAAoKioCC6XSxHSWlpaYDAYcnwHREREq4dhjHJGkiScOnVKNaRNTk4CALRaLVpbWxU7D7S1tcFsLuw5e0REtDEwjNG6ND4+rtrheerUKQDJDs+GhgbFnDS3243S0tLcXjwREdEyMIxRXpmZmUF3d7cipPX39yP1T7S2tlZ154Gqqio2DxAR0brDMEYFYXZ2VrXD0+v1Ih6PAwDKy8tVOzztdjtDGhER5QzDGBW0+fl5+Hw+RUjr7u7G3NwcAKC4uBhtbW2KkNbU1MQOTyIiWnUMY7QhJRIJnDx5UnVeWiAQAAAYjUY4nU5FSGttbYXRaMzxHRARUaFgGCNaQJIkDA8Pq3Z4jo+PA0h2eDY3NyuGPNva2mCxWHJ8B0RElG8YxogyNDExobo91NDQkHxOQ0OD6ry0srKyHF45ERGtZwxjRCsUCATQ3d2tCGp9fX0QRREAUFNTo9rhWVNTw+YBIqINjmGMaJXMzc2ht7dXMeTZ29uL+fl5AEBpaalqSLPb7RAEIcd3QEREa4FhjGiNzc/Po6+vTw5nCzs8I5EIAMBsNsuL2C4Mas3NzdDpuE8qEVEhYRgjWidEUYTf71ft8JyengYAGAwGtLa2KqppTqcTRUVFub0BIiI6JwxjROucJEkYGRlRDWmjo6MAAEEQ0NTUlLYtVKrD02q15vgOiIhoKQxjRHlsampKtcPT7/fL59jtdtUOz4qKihxeORERpTCMERWgUCgkd3guDGknTpyQOzyrq6tVQ1ptbS07PImI1hDDGNEGEo1G4fV6FcOdPT09iMViAACbzaZoHPB4PGhoaGCHJxHRKmAYIyLE43H09/crQlpXVxfC4TAAwGQyoa2tTRHSmpubodfrc3wHRET5i2GMiBYliiKGhoZUt4c6c+YMAECv16O1tTUtpLndbrhcLphMphzfARHR+scwRkTLJkkSxsbGVDs8h4eHAQAajQaNjY2KOWlut5vfH4iIFmAYI6Ksmp6eVg1pAwMD8jl1dXWqIa2qqip3F05ElCMMY0S0JsLhMHp6ehRBzefzIZFIAAAqKytVOzzr6urY4UlEBYthjIhyKhaLwev1KkJaT08PotEoAMBqtaqGNIfDAa1Wm+M7ICJaGYYxIlqXEokE+vv7Fd2dnZ2dCIVCAICioiK4XC5FUGtpaYHBYMjxHRARZYZhjIjyiiRJOHXqlGqH5+TkJABAp9OhpaVFEdJcLhfMZnOO74CIKB3DGBEVjPHxcdWQdvr0aQDJDk+Hw6E65Gmz2XJ89US0UTGMqQhH4xiYDCMWF2HQCXBUWGAx6nJ9WUR0jmZmZuRhzoVBrb+/H6lvbbW1tYqA5vF4UFVVxeYBIlpVDGPv8Y4G8dhRPzp6xuCfimDhzWoA1Jeb0e6qxoFd9WitsebqMokoiyKRCHp7exXVNK/Xi3g8DgAoLy9XDWlbtmxhSCOirNjwYWxwKoLDzxzDq74JaAUNEuLit5l6/PKWSjywfxvs5Zx7QlSI5ufn4fP5FI0D3d3dmJubAwAUFxfL66MtDGpNTU3s8CSiZdnQYeyJN/24/9njiIvSkiHsg7SCBjpBg2/s24pbd9av4hUS0XqSSCRw8uRJ1UVtA4EAAMBoNMLpdCqqaa2trTAajTm+AyJajzZsGPtehxcPvtC74tf5y2ucuKe9NQtXRET5SpIkDA8PqzYPjI+PAwC0Wi2am5sVIa2trQ0WiyXHd0BEuVQQYSzTeRsdHR1wOBxobGx8/6CghWC0QF++Gcb6bbBesBc6W7XiufHpUUz/9nFEB99FPDgJocgCfVkdihq24wcP/R1uWVAh6+rqwqFDh/Cb3/wGBoMB119/Pb7zne9wqxeiDWhiYkLRONDZ2YmhoSH5nIaGBtUOz7KyshxeORGtlYIIY48++mja5z/96U/x4osv4siRI2nH9+zZgxPDU7hshwdmz5UwNV0ESCLEuRBiI15Een4HaICKvV+ExXOl/Lz5M6cx8sghaPRGWLbvgc5WjURoCrGRE5jt+z2cX/1XvHToStjLzRgaGsKOHTtgs9nwxS9+EaFQCA8++CDq6+vxxhtvcCFKIgKQ/H7Z3d2tCGp9fX0QRREAsGnTJtWQVlNTw+YBogKSaX5a1+s63H777Wmfv/7663jxxRcVxwHg84+8BgAw1DSj+Lz2tMfiV4xh9MmvYeIXD0FfYYehpgkAEHjzXyHOz6Hurv+mqJolwtOIixIOP3MMRw7uwgMPPIBwOIzf//73qK9PVssuvvhi7NmzB4888gj+9E//NGv3TUT5q6SkBBdffDEuvvjitONzc3OKDs+XX34ZP/zhDzE/Pw8AKCsrUzQOeDwe2O12CIKQi9shojWwrsNYpryjQbw5cGbRx3W2alRe/yWMHLkXM0efQtW+ewEA8TPD0ForVYcvtZZSJEQJr/om4BsL4qmnnsLHPvYxOYgBwNVXXw2n04mf/exnDGNEtKSioiJs374d27dvTzs+Pz+Pvr6+tA7Pt99+G48//jgikQgAwGw2q4a0pqYm6HQF8W2caEMriP+LHzvqhyAsXdo31rmhK63F3MAf5GM6WzXmBv6A2YE/wuQ4X/V5WkGD//6LNzE2NoaLLrpI8fjFF1+M559/fkXXT0Qbl16vh8vlgsvlwv79++XjoijC7/cr5qQ999xzmJ6eBgAYDAY4nU5FUHM6nSgqKsrRHRHRchVEGOvoGYOYwRIW+qoGzHpfhxiNQDCaYf3QDQi/24GxJ/4r9NVNKKo/D0X121HUeAEEffIbWUKU8PIfkt2ZtbW1itesra3F1NQUotEo29uJKGsEQYDD4YDD4cDevXvl45IkYWRkRBHSfv3rX2N0dFR+blNTk2qHp9XKxa2J1pu8D2OhaBz+qUhG5wqGZMASY8kwZqhqQO1d38XMb5/ArO8NBN/qQ/CtZ6ExmFB21UFYL7gOADA8OQMAqmEr9dvn7OwswxgRrTqNRoPa2lrU1tbiqquuSntsampK0Tjw+OOPw+/3y+fY7XZFSHO73aioqFjrWyGi9+R9GDs5GUamy7qKseQK24Lh/RX29eV1qLzhK5DEBOYnBjF74g0EXn8KU//2PehKN8HkuADQJUNWX18fBgcHUVJSguLiYmi1WnnVbpPJlM3bIiJatvLyclx22WW47LLL0o6HQiG5wzMV0n7xi1/gu9/9rtzhWV1drdrhWVtbyw5PolWW92EsFhczPnd+/CQEcykEo3K7I42ghaHaAUO1A8bNbRh9/DDCx1+GyXEBtMXJNYG++tWv4qtf/ar8nOLiYiQSCQiCgMsvvxwlJSVpH1arVXFM7bjRaOQ3OyJaNcXFxbjooosU816j0Si8Xm/acOdvfvMb/PjHP0YsFgMA2Gw21ZDW0NDADk+iLMn7MGbQZfbNIHqqC/HpYVi2tp/1XENtcuX9RGgKAKCzVsJWWoYPXbgD9957L4LBIAKBAAKBAP72b/8WlZWVOP/88+XjY2Nj8uOpj0Qisej76fX6jIPbUsdT1ToiokwYjUacd955OO+889KOx+Nx9Pf3p4W0d955B08++STC4TCA5GhAW1ubIqg1NzdDr9fn4naI8lbehzFHhQVnqynFZ8Yw8YuHAa0OJbtuko/PDb4L4+Y2aLTpX4bZE28BSA5hAoAGwCdvvhn/69Ej2Lp1K+x2OwDgV7/6FSYnJ/Gtb30Ln//85xd9f0mSMDc3lxbOFga6xY6PjY3B5/OlHU99I1yMxWJZUaBLHSsqKmK1jmiD0ul0aG1tRWtrK2688Ub5uCiKGBoaUjQP/PKXv8SZM8nlhfR6PVpbWxUhzel0cjoH0SLyPoxZjDrUl5txYjr5eWz0BELvdiRX4I+GERv2ItLzGqABKj/2FRiq398yKfD6U4iN+GBy7YahyiE/P/zuv0MossK6M/lNqL7CjPv/41/hX55+Cu3t7fiLv/gLhEIh/MM//AO2bduGz372s0teo0ajgclkgslkQk1NzYruNx6PIxQKZRzoUh99fX2KY/F4fNH30el0GQe3pY5brVZW64gKhCAIqK+vR319Pa699lr5uCRJGBsbU4S0H/3oRxgeHgaQ/D7Y1NSkWIbD7Xav+z2PiVZb3ocxAGh3VaN/YAAAEOl8BZHOV5J7UxrM0JVvhnXnPtW9KUsu/RQinS9jbvBdhI+/DGk+Cm1xGczuK2C77FboSzdBK2jQ7qyG3W7HK6+8gi9/+cu477775L0p//Ef/3FNuyh1Oh1KS0tRWlq6otdJVeuWE+iCwSAmJiYUwS7Tat25BrrUB6t1ROuTRqNBTU0Nampq8JGPfCTtsenpaUWH589//nMMvPc9GwDq6uoUc9I8Hg8qKyvX9kaIcmRd702ZKe9oEHse/vWqvf5Lh65ASzXX5llMIpGQq3XLGYb94PGZmZmzVutWGuisViusVitXLSfKsXA4jJ6eHkU1zefzyXNsKysrVUPa5s2b+YsZ5YWC2Ch8Oe748VG81jeJRAaLv2ZKK2iwu6kCRw7uytpr0uIkSUI0Gl1RoEv9ORQKLfleZrN5xfPqSkpKYDKZ+EOBKItisRi8Xq+imtbd3Y1oNAoAsFqtqh2eDoeD0yJoXdlwYWxwKoKrH3oF0WUsdXE2Rp2Alw5dCXu5cikMWt9S1brlDsOqnZvaxFmNVqtdcaBLfc5qHdHiEokEBgYG0qpoqcAWDAYBJBfhdrlcipDW0tICg8GQ4zugjWjDhTEAeOJNP+57+ljWXu/vb9qGW3bWn/1EKmgLq3XnGuiCwaD8A2MxZrN5RYGO1TraiCRJwqlTpxTDnZ2dnZicnASQnOLQ0tKiqKa5XC6Yzfxlm1bPhgxjAPC9Di8efKF3xa9z7zUu/Kf2lixcEVGSKIoZza07WzVvZmZmyWqdIAgrmle3cG4d14uifDY+Pq4IaV1dXTh16hSAZOOBw+FQHfK02Ww5vnoqBBs2jAHJCtn9zx5HXJSWNYdMK2igEzT45r6trIjRupaq1q1kXl3q86WYTKYVBbrUn81mM6t1tG7MzMwotofq7OxEf38/Uj8SN2/erFiGw+PxoKqqiv+WKWMbOowByTlkh585hld9E9AKmiVDWerxy1sq8cD+bZwjRhuGKIoIh8NZGYZNTa5Wk6rWrXQYltU6Wk2zs7OqHZ5er1fu9C4vL1ft8NyyZQtDGils+DCW4h0N4rGjfnT0jsE/GUnbVFyD5IKu7c5q3H5JPZevIFqBaDSqWnlbbtUuGAxiqW9LRUVFK55XZ7VaYbFY+MOTMjI/Pw+fz6fa4Tk7Owsguf+n2oK2TU1N7PDcwBjGVISjcQxMhhGLizDoBDgqLLAY2cFGtJ58sFp3tkA3MzOTFuQWnne2at0Hd4s412FYVus2JlEUcfLkScWctM7OTgQCAQDJ/T+dTqeimtba2rqmC4bng0L8Gc0wRkQbXiwWW/G8uuVU61Y6DMtqXWGQJAnDw8OqzQNjY2MAksviNDc3K0JaW1sbLBZLju9g7cijVz1j8E+pjF6Vm9HuqsaBXfVorcm/0SuGMSKiLBFFEZFIZMXz6gKBAObm5hZ9H41Gs+JdJlLHua7W+jQ5Oaka0gYHB+VzGhoaVDs8y8rKcnjl2bVR5nUzjBERrUMfrNadS6BLfSz17dtoNJ5zoFt4zGKxQBCENfwKbUzBYFC1w7Ovrw+imFzMfNOmTaohraamJq8qqitd8eAb+7bi1jxZ8YBhjIiogEmShHA4nJVdJjKt1q10GJbVuuWbm5tDb2+vIqT19vbK6w2WlZUpGgc8Hg/sdvu6C9LZWgv0L69x4p721ixc0epiGCMioozMz89n1AmbyfFUFUeN0WjMyjAsq3VAPB7HiRMnFFtDdXV1IRKJAEju6qG2VlpTU1NOtl/biLvkMIwREdGakiQpbW7dSponUktGLEZtf9dzqdoVWkejKIoYHBxU7fCcnp4GABgMBjidTkVIczqd5/z16Ovrg8ViQU1Njerjy90/en5yCMG3f4nYcA+iIyeAxDzqPv9j6Erff/182D+aYYyIiPJWqlqXjWHYpap1BoPhnNeqW/h5cXHxuq7WSZKE0dFR1Y3WR0ZGACSXe2lqalLt8LRal+5kbGtrw+nTp/GTn/wEn/jEJxSP3/Hjo3itbzLjOWKhd17C5C+/C32lHdBoMT/WpwhjWkGD3U0VOHJw1zK+EmuLYYyIiDY8SZIwOzu7oqVNUsdTw3+LSe3pupJhWKvVCqPRuKYT8s+cOaPa4Xny5En5HLvdrghpbrcbFRUVCIfDsFqtckPJ5z73OTz88MPyEh3e0SD2PPzrZV1TYjYIjaCFYDRj5ujTmO74J0UYS3np0BXrdtH2TPNTfq+mRkREtASNRgOz2Qyz2YxNmzat6LXi8fhZO2HVjo+MjCiOJRKJRd9Hr9evONClqnWZrP5fVlaG3bt3Y/fu3WnHQ6EQenp60kLa888/j+9+97tytbG6ujotiAHAP/3TP6GjowM///nPsWPHDjx21K9YvkKMRjD96qOI9L6ORHgKgtECQ3UjSj/yGRg3tUBryixcaQUNHn3dj6/v25rR+esVwxgREVEGdDodysrKVrze18Jq3XKGW8fGxuDz+dKOna1aV1xcvKJhWLvdjq1bt6ZV66LRKLxerxzQfvKTn6S9pyiKOHHiBD70oQ/hb/7mb9Bh2K0Ynpz8v99HpOe3sF74Megr6yHOBhAd6sT85CCMm1oy/lomRAkdvWP4OhjGiIiIKEPZrtaFQqFlD8OOjY2lHZuZmTlrtW6x4DYxMaH6HEEQMJcA/FPKwDh74i1Yz78W5R/93IruHwD8kxGEo/G83jopf6+ciIhog9PpdCgtLUVpaemKXkeSJMzNzS17Xt3Y2Nii69QlEglMzAFqE9MFowXR072IByehs1as7NoBDEyGsXWzbUWvk0sMY0RERBucRqOByWSCyWRSLE8RjUaXDGmvvfaaalVNo9Ggzt4A+JXvV9b+WUz+4iGc+u+fhWFTM0xNF8Gy7aPQl55bpTCW4ZIZ6xXDGBERUYERRTHj4cuzVcBisdii76PVahddOkSn00GrUV+wweK+HEb7Vsz2/g6z/W8j8MbTCBx9ClX7D8PUfNGy79egW7/LimSCYYyIiGgdkCQJ0Wh0RYvlpj5CodCS72U2m1Un8Dscjown+xcVFWFoaAgHDhzAO++8o3iP+fl5/ObfnoHGfVB1qFJXXA7rhdfDeuH1SISnMfzIX2Dmdz9bdhjTAHBUWJb1nPWGYYyIiGgFEolEWhVqJZu/p/abVKPValXDUXl5uRyiMumQtFqty9oOKRKJyJuY/+pXv5K7KH0+36KT/jUaDb785S/j29/+Nq5++Dc4uWASvyQmIMXmIBS9H6C0llJoi8shxRe//8XUV5jzevI+wDBGREQbUKoKtZLqU+pYplWoD340NjYua6kJk8m0qovBTk9Pqy7+OjAwIJ9TV1cHj8eDa6+9FocOHYLb7YZOp8Nll10GIBkYN23ahCeffFI+1u6qxpGjJ+XlLaTYLIa+/xmYXZfBUN0IjaEIcwN/RGzYi7KrDgIAxLkwAr//PwCA6KkuAEDw/z0HjdECociCkg/dkHw/QYN2Z/WqfU3WCsMYERHljUQiseg2ScsNUmerQtlsNkU4qqysRFNT07K2ScrFptyLkSQJY2NjqqFreHgYQLKq1djYCI/Hg5tvvjltWySbTb1jsbKyEhMTE7jpppvwwx/+MK2788CuejzyuwH5c43eCOuFf4LZ/rcR6X0NkCToympRfs1/hPXCPwEAiHMhzLz6aNp7BN54BgCgLamWw1hClHD7Jet7s/BMcDskIiJaVR9cNmElQSocDi/5XhaLZVkbhS92vKioaE23JMo2SZJUNwzv6urC1NQUgOQE+9SG4Qu3OXK5XDCZTMt6v3/+53+GXq/Hbbfdpvp1W+7elJng3pRERFTwUlWolQ7jBQIBxOPxRd9Hp9NlvNXPUgHLarVmtP1PIYnH4+jv71cNXangajKZ0NbWlha4PB4Pmpubodfr1+Q6B6ciuPqhVxDN4hIURp2Alw5dCXu5OWuvmW2Z5qf1UzslIqIVS221k41hvLNVoYqLi1WDUWq/wkyD1FpvjJ2PotEoent75aCVCl09PT3y0hMlJSXweDzYtm0bbrnlFjl0NTQ0QBByu/SDvdyMb+zbivuePpa11/zmvq3rOogtB8MYEdE6sHAT6nOtPqU+ltrWRqfTwWazKcJRdXU1WlpaMq5IZboJNS1PKBSSOxcXVrv6+vrkv9fq6mq43W58+MMfxt133y1Xu2pra9d1qL11Zz0mQlE8+ELvil/r3mtcuGVn/s8VS2EYIyI6Rws3fD6X6tPC48vZ8HlhOKqpqcl4GI9VqPVjampKUeXq7OyE3//+cvV2ux1utxvXX3992hBjRcXKtg/KpXvaW1FZbMT9zx5HXJSWNYdMK2igEzT45r6tBRXEAM4ZI6INaGEVaqVB6mybK6t15C13PlRxcXHOh5lo+SRJwsjIiGrn4ujoKIDkZtpNTU1pc7lSnYtWqzXHd7B6BqciOPzMMbzqm4BW0CwZylKPX95SiQf2b8uroUlO4CeigiJJEiKRSFaG8WZnZ5d8r4WBaLmTyD84F4oKnyiK8Pv9qqFrenoaQDKYO51ORehyOp0oKirK7Q3kkHc0iMeO+tHROwb/ZCRtpX4Nkgu6tjurcfsl9Wipzr9wyjBGROvC/Px8WiBayVYvi+2BBwAGg2HF3XisQtFS4vE4Tpw4oQhc3d3d8jCz2WyWl4pYGLqam5vX1Xpj61E4GsfAZBixuAiDToCjwpL3K+uzm1JFIf5FE62GhVWoc60+pY4vVYXSaDTykgQfDEe1tbUZBamFc6GIsmFubg69vb2K0NXb2ysvFFtaWgqPx4MdO3bgwIEDcvCy2+0M8+fIYtRh62b1RWULXcFXxuQSaM8Y/FMqJdByM9pd1Tiwqx6tNflXAiVaaH5+PivDeMFgcMkqlNFozMownsVi4Q8uyplgMChPol8Yuvr6+uR//5s2bVJUuTweD2pqatgIQWe14YcpN8rkQMp/kiQhHA6f81pQCz+fm5tb9H1SVaiVDuNZrVZWoSivTE5OKqpcnZ2dGBoaks9paGhQBC63242ysrIcXjnluw0dxp5407+ittlv7NuKWwusbZayLxaLrWh18tTxs1WhioqKVtyNV1JSArPZzCoUFSxJknD69GnV5SLGx8cBJDsXW1paFKHL5XKhuLg4x3dAhWjDzhn7Xof3nBeUS7wX3u57+hgmQlHc096a5aujXBNFEeFwOCvLGkSj0UXfR6PRLBqOtmzZknGQslqtMBgMa/gVIlrfRFHEwMCAaudiIBAAkGzmcLlc8Hg8aG9vl4NXa2srq7q0Lq3rMJbpeHxHRwccDgcaGxvfPyhoIRgt0JdvhrF+G6wX7IXOVq14bnx6FNO/fRzRwXcRD05CKLJAX1aHv3p1O6oe+jt5Ybk33ngDjzzyCI4ePYp33nkH8XgcGRQVKUtisVhWhvGCweCSf29FRUWq4WjLli3LqkhZLBbOJyFagfn5efh8PkXo6unpkZtCLBaLXN36+Mc/LoeuxsZGdi5SXlnX/1qPHDmS9vlPf/pTvPjii4rjbrcbJ4aTu9CbPVfC1HQRIIkQ50KIjXgRfPNZBN96FhV7vwiL50r5efNnTmPkkUPQ6I2wbN8Dna0aidAUYiMnMPP6/8ZfP3sHdjdXwl5uxvPPP48f/ehH2L59O5qamtDbu/LtHApdqgp1rtWnhccyrUIt/LDZbLDb7RnPh2IVimjtzc7OoqenR1Hl8nq98ubi5eXl8Hg82LlzJ+688045dNntdv7SQwUhr+aM3XPPPfj+97+vWtnY//fP4F/uuwml7XfBtuumtMfiM2MYffJriM+MovbO78BQ0wQAmHzhfyD0h39D3Z/9T0XVLBGehsFaht1NFThycBdGR0dRUlICk8m05HUUgmg0mpVhvFAotOTXyGQyLXvek9pxs9nMb8hE61wgEFAdWuzv75e/T9TW1iom0Hs8HlRVVfH/ccpLG2rOmHc0iDcHziz6uM5Wjcrrv4SRI/di5uhTqNp3LwAgfmYYWmul6vCl1lKKhCjhVd8EfGNBtNTUrNr1Z4MoigiFQue0DtQHj8VisUXfRxAE1XBUWlqK+vr6jIOU1WqFXq9fw68QEa2F8fFx1dB16tQp+RyHwwGPx4P9+/fLocvtdqO0tDR3F06UQwURxh476ocgLP1bk7HODV1pLeYG/iAf09mqMTfwB8wO/BEmx/mqz9MKGjz6uh9f37c1m5csi0ajWRnGCwaDS76PyWRSDUepAJVpRYpVKCKSJAmnTp1SXS5icnISAKDVatHa2gq3243PfOYzaZ2LZjOXDyJaqCDCWEfPGMQMlrDQVzVg1vs6xGgEgtEM64duQPjdDow98V+hr25CUf15KKrfjqLGCyDok3uFJUQJHb1j+DreD2OpkvrQ0NCKg9RSVSitVqsajsrLy9HQ0LCsoT1OZiWi5UokEhgYGFCErq6uLvkXQKPRiLa2NrjdbuzZs0cOXS0tLZyDSZShvP8JHYrG4Z+KZHSuYEgGLDGWDGOGqgbU3vVdzPz2Ccz63kDwrT4E33oWGoMJZVcdhPWC6wAAJyfCaDvvfATPTMhzoQDAbrervo/ZbFYNRw6H46wLai78MJlMrEIR0aqLxWLwer2qnYup5hmr1SoHrU984hPynx0OB7RabY7vgCi/5X0YOzkZRqbT6MVYcnVywfB+iVxfXofKG74CSUxgfmIQsyfeQOD1pzD1b9+DrnQTTI4LAI0Gl133cdRbBVitVjz33HPo6OjAL3/5S9XVyVmFIqL1KBKJoLu7WxG6fD4fEokEAKCiogIejweXXnop7rrrLnlOV11dHX85JFoleZ8aYvHFVy7/oPnxkxDMpRCMyvkKGkELQ7UDhmoHjJvbMPr4YYSPv5wMYwDu+eKXsKM+uS1GX18fOjo6cN1112XlHoiIsml6eloxrNjZ2YmTJ0/K0yzq6urg8Xhw7bXX4ktf+pIcuqqqqnJ89UQbT96HMYMus+1doqe6EJ8ehmVr+9lfsza58n4iNLXs9yEiWguSJGFsbEy1c3F4eBhAcg2+xsZGeDwe3HzzzXLgamtrg81my/EdEFFK3ocxR4UFZyucx2fGMPGLhwGtDiUL1iCbG3wXxs1t0GjTvwyzJ94CkBzCBADNe+9DRLTWJEnC4OCgosrV1dWFqankL4w6nQ6tra3weDw4ePCgHLpcLhdMJlOO74CIzibvw5jFqEN9uRknppOfx0ZPIPRuR3IF/mgYsWEvIj2vARqg8mNfgaH6/S2TAq8/hdiIDybXbhiqHPLzw+/+O4QiK6w7bwQA1FeYMTFyCg+9t/L/W28lw9q3vvUtAEBDQwPuuOOOtblhIipI8Xgc/f39qp2L4XAYQHKJmlTn4t69e+VJ9M3NzVy3jyiP5X0YA4B2VzX6BwYAAJHOVxDpfCW5N6XBDF35Zlh37lPdm7Lk0k8h0vky5gbfRfj4y5Dmo9AWl8HsvgK2y26FvnQTtIIG7c5q9Pf342tf+1ra81OfX3nllQxjRJSRaDQKr9erCF09PT3yUjclJSXweDzYtm0bbrnlFjl0NTQ0QBA4ZYKo0OTVdkiL8Y4GsefhX6/a67906Aq0VFtX7fWJqPCEQiHVzsW+vj65c7G6ulqx9Y/b7UZtbS07F4kKwIbaDqm1xorLWyrxWt8kEhks/popraDB7qYKBjEiWtTU1JRiLldnZyf8fr98jt1uh9vtxvXXX58WuioqKnJ45US0XhREGAOAB/Zvw9UPvZLVMKYTNHhg/7asvR4R5SdJkjAyMqK6XMTo6CiA5L6tTU1NcLvduO222+TQ1dbWBquVv9AR0eIKJozZy834xr6tuO/pY1l7zW/u2wp7OfdQI9ooRFGE3+9XXS5ienoaAKDX6+F0OuHxePBnf/ZncpXL6XSiqKgotzdARHmpYMIYANy6sx4ToSgefKF3xa917zUu3LKzPgtXRUTrTTwex4kTJxShq7u7G5FIcns1s9mMtrY2eDwefOxjH5NDV3NzM3fZIKKsKrjvKPe0t6Ky2Ij7nz2OuCgta9hSK2igEzT45r6tDGJEBWBubg69vb2KKldvby/m5+cBAKWlpfB4PNixYwcOHDgAt9sNt9uN+vp6di4S0ZooiG5KNYNTERx+5hhe9U1AK2iWDGWpxy9vqcQD+7dxaJIozwSDQXR3dytCV19fH0QxuWVaTU2NomvR4/GgpqaGnYtEtCoyzU8FG8ZSvKNBPHbUj47eMfgnI2mbimuQXNC13VmN2y+pZ9ck0To3OTmpuijq4OCgfE5DQ0Na4Ep9lJeX5/DKiWgjYhhTEY7GMTAZRiwuwqAT4KiwwGIsuJFaorwmSRKGh4cVoauzsxPj4+MAkp2LLS0tctDyeDzweDxwuVwoLi7O8R0QESVtqHXGMmUx6rB1MzfHJVoPRFHEwMCA6nIRgUAAAGAwGOByueDxeNDe3i4Hr9bWVhiNxhzfARFRdmyoMEZEa29+fh4+n09R5erp6cHs7CwAwGKxyEHrxhtvlIcYGxsb2blIRAWP3+WIKCtmZ2fR09OjqHJ5vV7E43EAQHl5OdxuN3bu3Ik777xTDl1btmxh5yIRbVgMY0S0LIFAQHVR1P7+fqSmoNbW1sLj8eCjH/0ovvCFL8hVr6qqKnYuEhF9AMMYEakaHx9XDV2nTp2Sz3E4HPB4PNi/f39a92JpaWnuLpyIKM8wjBFtYJIk4dSpU6qdi5OTkwAArVaL1tZWuN1ufOYzn5GrXC6XC2Yz1+QjIlophjGiDSCRSKC/vz9tLlfqz8FgEABgNBrR1tYGt9uNPXv2yKGrpaUFBoMhx3dARFS4GMaICkgsFoPX61XtXIxGowAAq9UqB61PfOIT8p8dDge0Wm2O74CIaONhGCPKQ5FIBN3d3YrQ5fP5kEgkAACVlZVwu9249NJLcdddd8lzuurq6jiJnohoHWEYI1rHpqenVSfRDwwMyOfU1dXB4/Hg2muvxaFDh+RJ9FVVVbm7cCIiyhjDGFGOSZKEsbEx1dA1PDwMANBoNGhsbITH48HNN98sV7na2tpgs3FXCSKifMYwRrRGJEnC4OCgaufimTNnAAA6nQ5OpxNutxsHDx6UQ5fL5YLJZMrxHRAR0WpgGCPKsng8jv7+fkXo6urqQjgcBgCYTCa5c3Hv3r3yJPrm5mbo9foc3wEREa0lhjGicxSNRtHb26sIXD09PYjFYgAAm80Gt9uN7du345ZbbpFDV0NDA7f/ISIiAAxjRGcVCoVUOxdPnDgBURQBANXV1XC73fjwhz+Mu+++Wx5erK2tZeciEREtiWGM6D1TU1OKKldnZyf8fr98jt1uh9vtxvXXXy9XudxuNyoqKnJ45URElM8YxmhDkSQJIyMjqp2Lo6OjAABBENDU1ASPx4PbbrstrXPRarXm+A6IiKjQMIxRQRJFEX6/X1Hl6urqwvT0NABAr9fD6XTC4/HgiiuukEOX0+lEUVFRbm+AiIg2DIYxymvxeBwnTpxQVLm6u7sRiUQAAGazWV4I9YYbbpCHF5uamqDT8X8BIiLKLf4korwwNzeHnp4exZyu3t5ezM/PAwBKS0vh8XiwY8cOHDhwQA5ddrudnYtERLRuMYzRuhIMBhXDip2dnejv75c7Fzdt2gS3242PfOQj+PM//3M5dNXU1LBzkYiI8g7DGOXExMSEaugaGhqSz2loaIDb7caNN96Y1rlYVlaWwysnIiLKLoYxWjWSJOH06dOqnYvj4+MAAK1Wi+bmZng8Htxxxx1p2/8UFxfn+A6IiIhWH8MYrZgoihgYGFANXYFAAABgMBjgcrng8Xhw1VVXyaGrtbUVRqMxx3dARESUOwxjlLH5+Xn4fD7VzsW5uTkAgMVikYPW/v375S7GxsZGdi4SERGp2FA/HcPROAYmw4jFRRh0AhwVFliMG+pLkJHZ2Vn09PQoQpfX60U8HgcAlJeXw+PxYOfOnfj0pz8tz+nasmULJ9ETEREtQ8EnEe9oEI8d9aOjZwz+qQikBY9pANSXm9HuqsaBXfVordlYq6vPzMwoJtF3dXWhv78fkpT8Sm3evBlutxtXX301vvCFL8ihq6qqiqGLiIgoCzRS6qfuEgKBAGw2G2ZmZlBSUrIW17Vig1MRHH7mGF71TUAraJAQF7/N1OOXt1Tigf3bYC83r+GVrr7x8XFFlauzsxOnT58GAGg0GjgcjrSORY/Hg7a2NpSWlub24omIiPJUpvmpIMPYE2/6cf+zxxEXpSVD2AdpBQ10ggbf2LcVt+6sX8UrzD5JkjA0NKS60fXk5CQAQKfToaWlRRG6XC4XzObCCqBERES5lml+Krhhyu91ePHgC73n9NzEe+HtvqePYSIUxT3trVm+upVLJBLo7+9XVLm6u7sRDAYBAEVFRXLn4p49e+Tg1dLSAoPBkOM7ICIiooUKKow98ab/nIPYBz34Qi+qio24JUcVslgsBq/Xq6hy9fT0IBqNAgCsVivcbje2bt2KT37yk3Locjgc0Gq1ObluIiIiWp6CCWODUxHc/+xxAEDonZcw+fzD2PTph2CsXby6FfzjCwi88TTi06PQlVTC+qF9KLnoBvnxv372OHY3V67qHLJwOKzauejz+ZBIJAAAlZWV8Hg8uPTSS3Hw4EF5eHHz5s2cRE9ERJTnCiaMHX7mGOLLmB8WfPuXmPq/34fZtRslOz+O6OBxnHnpB5DiUdgu+SQAIC5KOPzMMRw5uGvF1zc9Pa26KOrAwIB8Tl1dHTweD6699locOnRIXqOrqqpqxe9PRERE61NBhDHvaBCv+iYyPl+cj2L610dgat6Jqv2HAQDWC64DIGHmt0+g+ILroC0qRkKU8KpvAr6xIFqqz77shSRJGBsbU+1cHBkZAZDsXGxqaoLb7canPvUpOXC53e68aI4gIiKi7CqIMPbYUf9Zl69YKOp/B+JsANYL/yTtuPXC6xE+/jJmfW+i+Lx2AMkOy0df9+Pr+7bK54miiMHBQdWNrs+cOQMA0Ov1aG1thdvtxt133y0PLTqdTphMpizdOREREeW7gghjHT1jy1rCIjbaBwAwbEqfT2bY1AJoBMRGTwDvhbGEKOG5t/thPP5sWudiOBwGAJhMJrS1tcHj8WDv3r3yJPrm5mbo9fos3SEREREVqrwPY6FoHP6pyLKekwhNARoBWktp2nGNVg/BZE0+vsD4LPD3/+1huFubcP755+PWW2+VQ1dDQwMEQVjpbRAREdEGlfdh7ORkGJnXxJLEeAwarXrVSqMzQIrH0o9pNPjtH3uxtc52jldJREREpC7vSzqxuLjs5wg6A6TEvOpjUjwGjU65MGossfz3ISIiIjqbvA9jBt3yb0FbXA5IIhLh6bTjUmIe4mww+XgW3oeIiIjobPI+YTgqLFjusqf66iYAQGzEm3Y8OuwDJBGGmqa045r33oeIiIgo2/I+jFmMOtQvc4X8oobtEIqsCP6/59OOh95+Hhq9EabmnWnH6yvMsBjzfnodERERrUMFkTDaXdU4cvSkYnmL0DsvYrbv94rzSy7ah9IrbsfUC/8D4898G0VNFyI6eBzh4x0oveJOaE3vL/CqFTRod1av+j0QERHRxlQQYezArno88rsBxfHQ288rTwZQvO1qWC+8HhC0CLzxL4j4jkJnrULZR++G9aJ9aecmRAm3X5KbzcKJiIio8BVEGGutseLylkq81jeJhCihePvVKN5+9VmfZ73guve2QVKnFTTY3VSR0VZIREREROci7+eMpTywfxt0wnKn8i9NJ2jwwP5tWX1NIiIiooUKJozZy834xoL9I7Phm/u2wr7M5gAiIiKi5SiYMAYAt+6sx19e48zKa917jQu37ORcMSIiIlpdBTFnbKF72ltRWWzE/c8eR1yUlrWBuFbQQCdo8M19WxnEiIiIaE0UVGUs5dad9Xjp0JXY3VQBIBmylpJ6fHdTBV46dCWDGBEREa2ZgquMpdjLzThycBe8o0E8dtSPjt4x+CcjaZuKa5Bc0LXdWY3bL6ln1yQRERGtOY0kSWcdxwsEArDZbJiZmUFJSclaXNeqCEfjGJgMIxYXYdAJcFRYuLI+ERERrYpM89OGSiIWow5bN9tyfRlEREREsoKcM0ZERESULzKqjKVGMgOBwKpeDBEREVGhSOWms80IyyiMBYNBAIDdbl/hZRERERFtLMFgEDbb4tOkMprAL4oiTp8+DavVCo0mu1sOERERERUiSZIQDAaxefNmCMLiM8MyCmNEREREtDo4gZ+IiIgohxjGiIiIiHKIYYyIiIgohxjGiIiIiHKIYYyIiIgohxjGiIiIiHKIYYyIiIgoh/4/fIbvJn3p4FMAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "f, ax = plt.subplots(figsize=(6,2), layout=\"constrained\")\n", + "simul_circuit.cg1.visualize_graph(ax=ax)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAADTCAYAAADNnRQhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAriUlEQVR4nO3de3Bc1Z0n8G+/pW633k/LLbUe/bjXmOBQLsCBEPPyGoLBEGIcA6FwPFl2YAonpJbATCDMhJAKu6ECk0yyQx6Ay6YygcTLTBLIRmE9ITiwhBnjllqtR6v1flpSq1vq1737h6KLxb0ysi35qlvfTxVlqe/t2+cIkL91zvmdY5BlWQYRERER6cKodwOIiIiI1jKGMSIiIiIdMYwRERER6YhhjIiIiEhHDGNEREREOmIYIyIiItIRwxgRERGRjsxLuUmSJPT398PpdMJgMKx0m4iIiIiynizLiEajWL9+PYzGxce/lhTG+vv74XK5lq1xRERERGtFT08PNmzYsOj1JYUxp9OpPKygoGB5WkZERESUw6ampuByuZQctZglhbH5qcmCggKGMSIiIqIz8FFLvLiAn4iIiEhHSxoZyxWxRBrhsRiSaQlWsxHuUgcctjX1IyAiIqJVJueTSGgoioPHImgODiMyHod8yjUDgNoSO7b5KrD3klp4Kk8/p0tERES03AyyLMsfddPU1BQKCwsxOTmZNWvGesbjePiV4zjaPgqT0YCMtHg3569f0VSGJ3ZtgqvEfh5bSkRERLloqfkpJ9eMHX47gmu+8wbe7BwDgNMGsVOvv9k5hmu+8wYOvx1Z8TYSERERATk4TflscwhPvdZ2Vu/NSDIykoyHXj6O0ekE7tvmWebWERERES2UUyNjh9+OnHUQ+7CnXmvDSxwhIyIiohW2akbGlnrMUnNzM9xuN+rr65XXzGYznAWFmLFXwOq6AM6LdsBcWKF6b3piCBN/OIREz/tIR8dgzHPAUlyDvLoLUXTFXtX9XztyAlsby067huwb3/gGjh07hmPHjmF4eBiPPvooHnvssSX1hYiIiGjVhLEXXnhhwffPP/88Xn/9ddXrgiBgZmYGALBnzx5cf/31kCQJz/zqPbQe/zOibx9B9J0jKN3xN3CIVyrvS53sx+BPDsBgscFx4bUwF1YgMz2O5GAHJt/6F80wlpZkPPzKcbyw75JF2/23f/u3qKqqwubNm/Gb3/zmXH4EREREtAatmjB2xx13LPj+rbfewuuvv656HQDC4TAA4OMf/zjuuOMOhIai+FpLKUrdV6Pwk8MYeunvMPqv34Gl1AVrZQMAYOrtX0JKzaLmnmdUo2aZ2IRmmzKSjKPto2gfjqKpQnvbi66uLrjdboyOjqK8vPwMe01ERERrXU6sGTt4LAKTcW6a01xYgbIbHgAyaUwe+7lyT/rkAEzOMs3pS5OjaMH3vd+7B8M/+zpme05g8Kdfgn9DGRoaGvD888+r3ut2u5ezK0RERLTG5EQYaw4OL9i+wlYjwFxUjdnwe8pr5sIKZKZGMBP+jyU9M3VyACO/+CZs7otQf8N/RXFxMe6++26cOHFiuZtPREREa1jWh7HpRBqR8bjqdUt5HaT4JKTE3DXnxTfCYLJg+PAj6P/R32D8tz9EvO0tSKlZzeemx3tRfvNDKL7yLqT91+Hnv3wVVqsVP/7xj1e0P0RERLS2rJo1Y2ereywGrS1djdY8AICUjMNos8NaXofqe76LyT8cxkz7nxB9pxPRd47AYM1H8VX74Lzovyx4v6WsFnmuCwAAMoCY0Q6fz4fOzs4V7hERERGtJVkfxpJpSfN1KTk34mW0frAthaWkBmU3fhmylEFqtAczHX/C1Fs/x/ivn4W5qAr57ouUe00FCxfjJ9MSiouLcfLkyeXvBBEREa1ZWT9NaTVrdyE10g2jvQhGm3qPMIPRBGuFG4WXfRbltzwCAIid+P3CewwLn/vkN/4BPT09GBsbQzgchiRph0AiIiKiM5H1I2PuUgcMwIKpykRfC9ITA3Bs3PaR77dWzx15lJkeX/wmWUboz28qIay+vh52ux1+vx+iKEIQBGzYsAEAGNKIiIjojGR9GHPYzKgtsaP7L4v405PDGP3XpwGTGQWX3KLcN9vzPmzr/TCYFnZ5puMdAHNTmIupK3PgjT/9EZ/61KcwOzuLRx99FIFAAC0tLWhpacGrr76KiYkJAHM78r/88stKSBMEAaIowuv1Ii8vb3k7T0RERFkvq8PYu+++ixdffBElfRGcCEQw29+GePBNwACUffrLsFZ8cGTS1Fs/R3KwHfm+rbCWuwEAyaEOxN7/HYx5Tji33KT5GSajAdu8H+xNlpeXhx07dmDHjh0A5k4OCIfDGBkZwTPPPIOmpiY4HA688847aG5uxujoKADAaDSioaFBCWfzf/r9fjid2hvKEhERUe7L6jB26NAhHDp0CGazGZI5H+aS9XBu2al5NmXBZZ9FPPB7zPa8j9iJ30NOJWBaVwy78EkUfuJ2WIqqND8jI8m449LaRdvw3HPP4Y033lC+b2v74KDy5uZmXHjhhcoI2vxo2qFDhxCJfHAIucvlWhDS5r8uLS092x8NERERZQmDLMtaO0MsMDU1hcLCQkxOTqKgoOB8tOuM3fncMbzZObZg89dzZTIasLWh9LRnU56t6elpBINBJaDN/9ne3q6sOysvL18wijb/Z3V19ZIPViciIiJ9LDU/5UwY6xmP45rvvIHEIltdnA2b2YjfHrgSrhJ1ReZKSSQSCIVCqpAWDAaRTCYBAAUFBZohra6uDkZj1hfIEhER5YQ1F8YA4PDbETz08vFle963btmE3VsWn6I8n9LpNLq6uhYUDsx/HYvFAAD5+fnw+XyqoNbU1ASLxaJzD4iIiNaWNRnGAODZ5hCeeq3to2/8CF+5zoe/3ta0DC1aWZIkobe3d0E4CwQCCAQCyga1ZrMZHo9HFdJ8Ph/y8/N17gEREVFuWrNhDJgbIXv0yAmkJfmM1pCZjAaYjQY8vnPjqhkRO1uyLGN4eFgV0lpaWjAwMAAAMBgMqK+vVxUOCIKQFf+eiYiIVrM1HcaAuTVkD79yHEfbR2EyGk4byuavX9FUhid2bTqva8T0MDExoZrqDAQCCIfDyj01NTWqNWmCIKC8vHzxBxMREZFizYexeaGhKA4ei6C5bRiRsfiCnfoNAGpL7djmrcAdl9aiqWJt7/cVi8UQDAZVIa29vR2ZTAYAUFZWphnSampqWOFJRER0CoYxDbFEGuGxGJJpCVazEe5SBxy2rN5q7bxIJpMIhUKq0bTW1lYkEgkAgNPpVE11iqIIt9sNk8mkcw+IiIjOP4YxWnGZTAbhcFi1Ji0QCGB6ehrA3IkFPp9PFdKamppgtVp17gEREdHKYRgj3ciyjL6+Ps2QNjY2BgAwmUzweDyqkwf8fj/s9txes0dERGsDwxitSiMjI5oVnn19fQDmKjzr6upUa9IEQUBRUZG+jSciIjoDDGOUVSYnJ9Ha2qoKaV1dXZj/T7S6ulrz5IHy8nIWDxAR0arDMEY5YWZmRrPCMxQKIZ1OAwBKSko0KzxdLhdDGhER6YZhjHJaKpVCe3u7KqS1trZidnYWALBu3Tr4/X5VSGtoaGCFJxERrTiGMVqTMpkMuru7NdelTU1NAQBsNhu8Xq8qpHk8HthsNp17QEREuYJhjOgUsixjYGBAs8JzZGQEwFyFZ2Njo2rK0+/3w+Fw6NwDIiLKNgxjREs0OjqqeTxUb2+vck9dXZ3murTi4mIdW05ERKsZwxjROZqamkJra6sqqHV2dkKSJABAZWWlZoVnZWUliweIiNY4hjGiFTI7O4u2tjbVlGdbWxtSqRQAoKioSDOkuVwuGI1GnXtARETnA8MY0XmWSqXQ2dmphLNTKzzj8TgAwG63K5vYnhrUGhsbYTbznFQiolzCMEa0SkiShEgkolnhOTExAQCwWq3weDyq0TSv14u8vDx9O0BERGeFYYxolZNlGYODg5ohbWhoCABgNBrR0NCw4Fio+QpPp9Opcw+IiOh0GMaIstj4+LhmhWckElHucblcmhWepaWlOraciIjmMYwR5aDp6WmlwvPUkNbR0aFUeFZUVGiGtOrqalZ4EhGdRwxjRGtIIpFAKBRSTXcGg0Ekk0kAQGFhoapwQBRF1NXVscKTiGgFMIwREdLpNLq6ulQhraWlBbFYDACQn58Pv9+vCmmNjY2wWCw694CIKHsxjBHRoiRJQm9vr+bxUCdPngQAWCwWeDyeBSFNEAT4fD7k5+fr3AMiotWPYYyIzpgsyxgeHtas8BwYGAAAGAwG1NfXq9akCYLA3w9ERKdgGCOiZTUxMaEZ0sLhsHJPTU2NZkgrLy/Xr+FERDphGCOi8yIWiyEYDKqCWnt7OzKZDACgrKxMs8KzpqaGFZ5ElLMYxohIV8lkEqFQSBXSgsEgEokEAMDpdGqGNLfbDZPJpHMPiIjODcMYEa1KmUwGXV1dqurOQCCA6elpAEBeXh58Pp8qqDU1NcFqtercAyKipWEYI6KsIssy+vr6NCs8x8bGAABmsxlNTU2qkObz+WC323XuARHRQgxjRJQzRkZGNENaf38/gLkKT7fbrTnlWVhYqHPriWitYhgjopw3OTmpTHOeGtS6urow/6uturpaFdBEUUR5eTmLB4hoRTGMEdGaFY/H0dbWphpNC4VCSKfTAICSkhLNkLZhwwaGNCJaFgxjREQfkkql0N7eriocaG1txezsLABg3bp1yv5opwa1hoYGVngS0RlhGCMiWqJMJoPu7m7NTW2npqYAADabDV6vVzWa5vF4YLPZdO4BEa1GDGNEROdIlmUMDAxoFg+MjIwAAEwmExobG1Uhze/3w+Fw6NwDItITwxgR0QoaHR1VFQ4EAgH09vYq99TV1WlWeBYXF+vYciI6XxjGiIh0MDU1hdbWVlVQ6+zshCRJAICqqirNkFZZWcniAaIcwjBGRLSKzM7OalZ4trW1IZVKAQCKi4tVhQOiKMLlcsFoNOrcAyI6UwxjRERZIJVKobOzU7PCMx6PAwDsdrtmSGtoaIDZbNa5B0S0GIYxIqIsJkkSIpGIZoXnxMQEAMBqtcLr9aqCmtfrRV5enr4dICKGMSKiXCTLMgYHBzVD2tDQEADAaDSioaFBs8LT6XTq3AOitYNhjIhojRkfH9es8IxEIso9LpdLFdIEQUBpaamOLSfKTQxjREQEAJienlYqPE8NaR0dHUqFZ0VFhWaFZ3V1NSs8ic4SwxgREZ1WIpFAKBRSTXcGg0Ekk0kAQGFhoWZIq6urY4Un0UdgGCMiorOSTqfR1dWlCmktLS2IxWIAgPz8fPj9flVQa2xshMVi0bkHRKsDwxgRES0rSZLQ29urmu4MBAI4efIkAMBiscDj8ahCmtfrRX5+vs49IDq/GMaIiOi8kGUZw8PDmhWeAwMDAACDwYCGhgbVNhyCIPDvFcpZDGNERKS7iYkJzQrPcDis3FNTU6NakyaKIsrKyvRrONEyYBgjIqJVKxaLIRgMqkJae3s7MpkMAKCsrEwzpK1fv54VnpQVGMaIiCjrJJNJhEIh1Whaa2srEokEAMDpdGpWeLrdbphMJp17QPQBhjEiIsoZmUwG4XBYs8IzGo0CAPLy8uDz+VQhrampCVarVece0FrEMEZERDlPlmX09fVpVniOjY0BAMxmM5qamlSjaT6fD3a7XeceUC5jGCMiojVtZGREs8Kzr68PwFyFp9vt1pzyLCws1Ln1lAsYxoiIiDRMTk5qHg/V1dWF+b8S169fr9qGQxRFlJeXs3iAloxhjIiI6AzMzMxoVniGQiGk02kAQElJiWaF54YNGxjSSIVhjIiIaBmkUim0t7drVnjOzMwAANatW6e5oW1DQwMrPNcwhjEiIqIVJEkSuru7VWvSAoEApqamAAA2mw1er1c1mubxeGCz2XTuweoSS6QRHoshmZZgNRvhLnXAYTPr3axzwjBGRESkA1mWMTAwoFk8MDw8DAAwmUxobGxUhTS/3w+Hw6FzD86f0FAUB49F0BwcRmQ8jlMDiQFAbYkd23wV2HtJLTyVTr2aedYYxoiIiFaZsbExzZDW09Oj3FNXV6dZ4VlcXKxjy5dXz3gcD79yHEfbR2EyGpCRFo8i89evaCrDE7s2wVWSPduRMIwRERFliWg0qlnh2dnZCUmSAABVVVWaIa2ysjKrigcOvx3Bo0dOIC3Jpw1hH2YyGmA2GvD1nRtx+5baFWzh8mEYIyIiynKzs7Noa2tThbS2tjakUikAQHFxsapwQBRFuFwuGI1GnXuw0LPNITz1Wts5P+fB67y4b5tnGVq0shjGiIiIclQ6nUZHR4fqaKiWlhbE43EAgN1u19wrraGhAWbz+V8Yf/jtCB56+fii1yeOHsTkHw6h7qFXl/S8b92yCbtX+QjZUvNTdpcpEBERrUFmsxk+nw8+nw8333yz8rokSejp6VGtSXv11VcxMTEBALBarfB6vaqQ5vV6z7rCs7OzEw6HA5WVlZrXe8bjePTIiTN65kzXu4i1HEWyP4jUWC9MzjJs+G8/Uq5/7cgJbG0sy6o1ZIthGCMiIsoRRqMRdXV1qKurw44dO5TXZVnG0NCQKqT94Ac/wODgoPLehoYGzQpPp/P0lYzXX389+vv78eMf/xi33nqr6vrDrxxH+gzWhwFA7MQbiLcehbWyEaZ1JarraUnGw68cxwv7Ljmj565GDGNEREQ5zmAwoKqqClVVVbjqqqsWXDt58qRqTdrhw4fR3d2t3ONyuVQhTRAElJaWIpFIIBQKQZIkfOYzn8EXvvAFPP3008oWHaGhKI62j55xm4uuvAulO+6HwWTG8M++juRI94LrGUnG0fZRtA9H0VSRfdtenIphjIiIaA0rLi7G1q1bsXXr1gWvT09PIxgMLghp//Zv/4bvfve7SoVnRUUFamtrle8B4Ec/+hGam5vxs5/9DJs3b8bBYxHV9hWzPSdw8v/8M5IjYZidpSi4RD2aZnaWfmTbTUYDXnwrgsd2bjzb7q8KDGNERESksm7dOlx88cW4+OKLF7w+PxI2H9B+85vfLLguSRI6Ojpw8cUX4+///u/RbN26IIglh8MYfulrMNoLUHT55yBLGUz8+0GY7EVn3MaMJKO5bRiPgWGMiIiI1gibzYYLLrgAF1xwgfLan/70J2QyGQBza88kSYLdbgcseYiMxxe8f+LoiwBkVO39FsyFFQAAh+8T6H/ur8+qPZGxOGKJdFYfnbS6NiAhIiKirNLS0qIEsc2bN+Pxxx/H+++/j2g0ipvv+MKCI45kKYPZrj8j33OpEsQAwFLmQn7Dx8/q82UA4bHYOfRAfwxjREREdNai0SgA4GMf+xhuv/123Hbbbdi4cSMMBgOSaWnBvVJ8CnI6AUvxetVzzCU1Z92GD39OtmEYIyIiorPm9/thMBjwn//5n/jqV78Kn88Ht9uN3bt3oyfcdV7aYDVnd5zJ3glWIiIi0k0ymUQoFEI8Hsf8YT7zf3Z3d6O7uxvTiRQMwj5lqtJoL4DBbEPqZL/qeenxvrNqhwGAu9RxVu9dLRjGiIiIaFHxeFzzEPP29nZlrdip5g8t/9KXvoRvfvObuObpf0f3XxbxG4wm5NVvxkzoLaQnh5V1Y6nRHsx0vntW7asttWf14n2AYYyIiIgATExMqAJXS0sLwuGwck9NTQ1EUcT27dtx4MABCIKA6upqeL1eAIDJZEJVVRVeeuklfOITnwAAbPNV4IVj3cr2FkVX7MVg17sYPPjf4fz4DYCUwdT/+9+wlNUiNfLBZyWHuxAPHQMApE4OQE7EMPGHwwAAa0U97J5LYDIasM37QSFAtmIYIyIiWiNkWcbw8LBm6BoYGAAwN7JVX18PURRx2223LTgWqbCwUPO5ZWVlGB0dxS233IIf/vCHKCoqUq7tvaQWP/ljWPneWlGPis8+jpO/+2dMHH0RZmcZii7fi8z0OCZPDWODHZg8+uKCz5n/3nHB1bB7LkFGknHHpav7sPClMMjzE7ynsdRTx4mIiEh/sixrHhje0tKC8fFxAHOHjc8fGH7qMUc+nw/5+fln9Hk//elPYbFYsGfPHmWa8lR3PncMb3aOLdj89VyZjAZsbShd1WdTLjU/MYwRERFlqXQ6ja6uLs3QFYvN7b2Vn58Pv9+/IHCJoojGxkZYLJbz0s6e8Tiu+c4bSCzjFhQ2sxG/PXAlXCX2ZXvmcltqfuI0JRER0SqXSCTQ1tamBK350BUMBpFMJgEABQUFEEURmzZtwu7du5XQVVdXB6NR360fXCV2fH3nRjz08vFle+bjOzeu6iB2JhjGiIiIVonp6WnNysXOzk6lcrGiogKCIODyyy/H/v37ldGu6upqzSnC1eL2LbUYnU7gqdfazvlZX7nOh91bsn+t2DyGMSIiovNsfHxcNcoVCAQQiUSUe1wuFwRBwA033LBgirG0tFTHlp+b+7Z5ULbOhkePnEBaks9oDZnJaIDZaMDjOzfmVBADuGaMiIhoRciyjMHBQc3KxaGhIQBzh2o3NDQsWMs1X7nodDp17sHK6RmP4+FXjuNo+yhMRsNpQ9n89SuayvDErk1ZNTXJBfxERETngSRJiEQimqFrYmICAGCxWOD1elWhy+v1Ii8vT98O6Cg0FMXBYxE0tw0jMhZfcKi4AXMbum7zVuCOS2vRVJF94ZRhjIiIaBml02l0dHSoAldrayvi8bkd5u12u7JVxKmhq7GxEWYzVwadTiyRRngshmRagtVshLvUkfU767OakoiI6CzMzs6ira1NFbra2tqQSqUAAEVFRRBFEZs3b8bevXuV4OVyuXSvXMxWDpsZG9drbyqb6xjGiIhoTYpGo8oi+lNDV2dnJyRpbj+sqqoqCIKAK6+8Evfee68SuiorK1d15SJlF4YxIiLKaWNjY6pRrkAggN7eXuWeuro6CIKAnTt3LpheLC4u1rHltFYwjBERUdaTZRn9/f2a20WMjIwAmKtcbGpqgiAIuPPOO5XQ5fP5sG7dOp17QGsZwxgREWUNSZIQDoc1KxenpqYAAFarFT6fD6IoYtu2bcool8fjgc1m07kHRGoMY0REtOqkUim0t7erQlcwGMTMzAwAwOFwKKNbN998sxK66uvrWblIWYX/tRIRkW5mZmYQDAZVo1yhUAjpdBoAUFJSAlEUsWXLFtx1111K6HK5XFxETzmBYYyIiFbc1NSU5tRiV1cX5re7rK6uhiiKuPrqq3H//fcro17l5eUMXZTTGMaIiGjZjIyMaIauvr4+5R632w1RFLFr1y5llEsQBBQVFenXcCIdMYwREdEZkWUZfX19mttFjI2NAQBMJhM8Hg8EQcDdd9+9oHLRbs+eswWJzgeGMSIi0pTJZBAOh1Whq6WlBdFoFABgs9ng9/shCAKuvfZaJXQ1NTXBarXq3AOi7MAwRkS0xiWTSYRCIc3KxUQiAQBwOp1K0Lr11luVr91uN0wmk849IMpuDGNERGtEPB5Ha2urKnS1t7cjk8kAAEpLSyGKIi677DLcc889ypqumpoaLqInWiEMY0REOWZiYkI1rRgIBNDd3a1ULtbU1EAURWzfvh0PPPCAErrKy8t1bj3R2sMwRkSUhWRZxvDwsGbl4sDAAADAYDCgvr4eoijitttuUwKX3+9HYWGhzj0gonkMY0REq5gsy+jp6VGNcrW0tGB8fBwAYDab4fF4IIoi9u3bp4Qun8+H/Px8nXtARB+FYYyIaBVIp9Po6urSrFyMxWIAgPz8fKVycceOHcoi+sbGRlgsFp17QERni2GMiOg8SiQSCIVCqtAVDAaRTCYBAAUFBRBFEZs2bcLu3buV0FVXVwej0ahzD4houTGMERGtgOnpac3Kxc7OTqVysaKiAoIg4PLLL8f+/fuV6cXq6mpWLhKtIQxjRETnYHx8XLWWKxAIIBKJKPe4XC4IgoAbbrhBGeUSBAGlpaU6tpyIVguGMSKijyDLMgYHBzW3ixgaGgIAGI1GNDQ0QBAE7NmzRwldfr8fTqdT5x4Q0WrGMEZE9BeSJCESiWhuFzExMQEAsFgs8Hq9EEURX/ziF5VRLq/Xi7y8PH07QERZiWGMiNacdDqNjo4OVehqbW1FPB4HANjtdvj9foiiiE9/+tNK6GpsbITZzF+dRLR8+BuFiHLW7Ows2traVKNcbW1tSKVSAICioiKIoojNmzdj7969EAQBgiCgtraWlYtEdF4wjBFR1otGo2htbVWFrs7OTkiSBACorKyEKIq48sorce+99ypruiorK1m5SES6YhgjoqwxNjamuSlqT0+Pck9dXR0EQcDOnTuVUS5BEFBSUqJjy4mIFscwRkSriizLGBgYUIWuQCCAkZERAHOVi01NTRAEAXv37oUoihBFET6fD+vWrdO5B0REZ4ZhjIh0IUkSwuGw5nYRU1NTAACr1QqfzwdRFLFt2zZlatHj8cBms+ncAyKi5cEwRkQrKpVKob29XTXKFQwGMTMzAwBwOBxK0LrpppuUysX6+npWLhJRzuNvOSJaFjMzMwgGg6pRrlAohHQ6DQAoKSmBIAjYsmUL7rrrLiV0bdiwgZWLRLRmMYwR0RmZmprS3BS1q6sLsiwDAKqrqyGKIq6++mrcf//9yqhXeXk5KxeJiD6EYYyINI2MjGiGrr6+PuUet9sNURSxa9cuZZRLEAQUFRXp13AioizDMEa0hsmyjL6+Ps3KxbGxMQCAyWSCx+OBIAi4++67lVEun88Hu92ucw+IiLIfwxjRGpDJZNDV1bVgLdf819FoFABgs9ng9/shCAKuvfZaJXQ1NTXBarXq3AMiotzFMEaUQ5LJJEKhkGblYiKRAAA4nU4laN16663K1263GyaTSeceEBGtPQxjRFkoHo+jtbVVFbra29uRyWQAAGVlZRAEAZdddhnuueceZU1XTU0NF9ETEa0iDGNEq9jExITmIvpwOKzcU1NTA1EUsX37dhw4cEBZRF9eXq5fw4mIaMkYxoh0JssyhoeHNUPXwMAAAMBgMKC+vh6iKOK2225TRrn8fj8KCwt17gEREZ0LhjGi80SWZfT09GhWLp48eRIAYDab4fV6IQgC9u3bp4Qun8+H/Px8nXtAREQrgWGMaJml02l0dXWpQldLSwtisRgAID8/X6lc3LFjh7KIvrGxERaLReceEBHR+cQwRnSWEokE2traVIErGAwimUwCAAoLCyEIAi688ELs3r1bCV11dXU8/oeIiAAwjBF9pOnpac3KxY6ODkiSBACoqKiAIAi4/PLLsX//fmV6sbq6mpWLRER0WgxjRH8xPj6uGuUKBAKIRCLKPS6XC4Ig4IYbblBGuQRBQGlpqY4tJyKibMYwRmuKLMsYHBzUrFwcGhoCABiNRjQ0NEAURezZs2dB5aLT6dS5B0RElGsYxignSZKESCSiGuVqaWnBxMQEAMBiscDr9UIURXzyk59UQpfX60VeXp6+HSAiojWDYYyyWjqdRkdHh2qUq7W1FfF4HABgt9uVjVBvvPFGZXqxoaEBZjP/FyAiIn3xbyLKCrOzswgGg6o1XW1tbUilUgCAoqIiiKKIzZs3Y+/evUrocrlcrFwkIqJVi2GMVpVoNKqaVgwEAujq6lIqF6uqqiAIAj71qU/h3nvvVUJXZWUlKxeJiCjrMIyRLkZHRzVDV29vr3JPXV0dBEHATTfdtKBysbi4WMeWExERLS+GMVoxsiyjv79fs3JxZGQEAGAymdDY2AhRFHHnnXcuOP5n3bp1OveAiIho5TGM0TmTJAnhcFgzdE1NTQEArFYrfD4fRFHEVVddpYQuj8cDm82mcw+IiIj0wzBGS5ZKpdDe3q5ZuTg7OwsAcDgcStDatWuXUsVYX1/PykUiIiINa+pvx1gijfBYDMm0BKvZCHepAw7bmvoRLMnMzAyCwaAqdIVCIaTTaQBASUkJRFHEli1b8PnPf15Z07VhwwYuoiciIjoDOZ9EQkNRHDwWQXNwGJHxOORTrhkA1JbYsc1Xgb2X1MJTubZ2V5+cnFQtom9paUFXVxdkee4ntX79egiCgGuuuQb333+/ErrKy8sZuoiIiJaBQZ7/W/c0pqamUFhYiMnJSRQUFJyPdp2znvE4Hn7lOI62j8JkNCAjLd7N+etXNJXhiV2b4Cqxn8eWrryRkRHVKFcgEEB/fz8AwGAwwO12L6hYFEURfr8fRUVF+jaeiIgoSy01P+VkGDv8dgSPHjmBtCSfNoR9mMlogNlowNd3bsTtW2pXsIXLT5Zl9Pb2ah50PTY2BgAwm81oampShS6fzwe7PbcCKBERkd6Wmp9ybpry2eYQnnqt7azem/lLeHvo5eMYnU7gvm2eZW7ductkMujq6lKNcrW2tiIajQIA8vLylMrFa6+9VgleTU1NsFqtOveAiIiITpVTYezw25GzDmIf9tRrbShfZ8NunUbIkskkQqGQapQrGAwikUgAAJxOJwRBwMaNG/GZz3xGCV1utxsmk0mXdhMREdGZWTVhbKmLwZubm+F2u1FfX6+8Zjab4SwoxIy9AlbXBXBetAPmwgrVe9MTQ5j4wyEket5HOjoGY54DluIa5NVdiKIr9qru/9qRE9jaWHbaNWSSJOGpp57C97//fQwMDMDr9eKrX/0q9uzZs6T+xGIxzcrF9vZ2ZDIZAEBZWRlEUcRll12Gffv2KdOL69ev5yJ6IiKiLLdqwtgLL7yw4Pvnn38er7/+uup1QRAwMzMDANizZw+uv/56SJKEZ371HlqP/xnRt48g+s4RlO74GzjEK5X3pU72Y/AnB2Cw2OC48FqYCyuQmR5HcrADk2/9i2YYS0syHn7lOF7Yd8mi7X7kkUfw5JNPYv/+/diyZQt++ctf4nOf+xwMBgNuv/125b6JiQnNTVHD4bByT01NDURRxPbt23HgwAFlj67y8vIz+lkSERFR9li1C/jvu+8+/OM//iO0mhcOh1FfX49vf/vbePDBBxEaiuLap/8vACA9OYyhl/4O6ckhVN/1P2GtbAAAjL32fUy/92vUfPF/qUbNMrEJmBxFi7bltwc+iaYK9bYXfX19qK+vx1/91V/h2WefhSzLGBoawnXXXYeenh7s2bNHGfUaHBwEMDcC2NDQsGAR/fw/2VAcQUREREuzphbwHzwWUbanMBdWoOyGBzD4wlcweeznKN/5FQBA+uQATM4yzenLDwex3u/dA2t5HQou/Qwmfvcc/P8jjNoNNXjsscdw1113QZIk9PT04Nvf/jZSqRQGBwdx+eWXIxAI4OTJk8pzfvWrX+Hiiy/G/v37lfDl9XqRn5+/oj8PIiIiyh45Ecaag8MLtrCw1QgwF1VjNvye8pq5sAKz4fcwE/4P5Ls/9pHPTJ0cwMgvvol1F16His1XIfHeEXz+85/Ht771LXR3dyMWiyn3dnR0YOPGjdixYwdEUYTD4cD27dvxpS99Cffff/+y9pWIiIhyS9aHselEGpHxuOp1S3kdZkJvQUrEYbTZ4bz4RsTeb8bw4UdgqWhAXu0FyKu9EHn1F8FoyVO9Pz3ei8q9TyLPdQFkWUb0j/8Cg8EAo9GIxx57DKIo4qmnnkJ3dzf+/Oc/L3hvPD7XnvlNVYmIiIgWY9S7AeeqeywGrUVvRutcwJKSc8HIWl6H6nu+C8fGbchMDiH6zhGMvPwP6H3mTkTf+7Xq/ZayWuS5LgAwt87rj8fbsWnTJng8Hjz44IO4/vrrYTAYYLPZVO/Ny5v77PlCAyIiIqLFZP3IWDItab4uJWcBAEbrB9tSWEpqUHbjlyFLGaRGezDT8SdMvfVzjP/6WZiLqpDvvki511SwsIIxmZFQXFy8YE1Yfn6+sufXqWZnZ5XrRERERKeT9SNjVrN2F1Ij3TDai2C0qfcIMxhNsFa4UXjZZ1F+yyMAgNiJ3y+8x7DwufOfc2p1Z3V1NQYHB1UVnwMDAwDmDtkmIiIiOp2sD2PuUgc+vO1poq8F6YkB5Ndv/sj3W6vnjjzKTI8veo/hL5/zYRdddBHi8ThaWloWvH7s2DHlOhEREdHpZH0Yc9jMqD1lh/z05DBG//VpwGRGwSW3KK/P9rwPOZNWvX+m4x0Ac1OYi6kttcNhU8/o3nTTTbBYLPje976nvCbLMv7pn/4JNTU12Lp169l0iYiIiNaQrF4z9u677+LFF19ESV8EJwIRzPa3IR58EzAAZZ/+MqwVHxyZNPXWz5EcbEe+byus5W4AQHKoA7H3fwdjnhPOLTdpfobJaMA2r3pvMgDYsGEDHnjgAWW/sS1btuAXv/gFjh49ioMHD/J8SCIiIvpIWR3GDh06hEOHDsFsNkMy58Ncsh7OLTs1z6YsuOyziAd+j9me9xE78XvIqQRM64phFz6Jwk/cDktRleZnZCQZd1y6+GHhTz75JIqLi/GDH/wAP/nJT+DxePDiiy/ic5/73LL2lYiIiHLTqj0O6Uzd+dwxvNk5tmDz13NlMhqwtaH0tGdTEhEREWlZan7K+jVj857YtQlm44eX8p8bs9GAJ3ZtWtZnEhEREZ0qZ8KYq8SOr+/cuKzPfHznRrhK1FtjEBERES2XnAljAHD7llo8eJ13WZ71let82L1l8bViRERERMshqxfwa7lvmwdl62x49MgJpCX5jNaQmYwGmI0GPL5zI4MYERERnRc5NTI27/YttfjtgSuxtaEUwFzIOp3561sbSvHbA1cyiBEREdF5k3MjY/NcJXa8sO8ShIaiOHgsgua2YUTG4gsOFTdgbkPXbd4K3HFpLZoqnHo1l4iIiNaonNnaYiliiTTCYzEk0xKsZiPcpQ7NnfWJiIiIztVS89OaSiIOmxkb1xfq3QwiIiIiRU6uGSMiIiLKFksaGZufyZyamlrRxhARERHlivnc9FErwpYUxqLRKADA5XKdY7OIiIiI1pZoNIrCwsWXSS1pAb8kSejv74fT6YTBsLxHDhERERHlIlmWEY1GsX79ehiNi68MW1IYIyIiIqKVwQX8RERERDpiGCMiIiLSEcMYERERkY4YxoiIiIh0xDBGREREpCOGMSIiIiIdMYwRERER6ej/A5AD0aKABDiyAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "f, ax = plt.subplots(figsize=(6,2), layout=\"constrained\")\n", + "simul_circuit.cg2.visualize_graph(ax=ax)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Speed Benchmark**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To benchamrk the performance of graph based calculation, We consider the three different types of model that are typically used in EIS and 2nd-NLEIS analysis. \n", + "\n", + "This includes: \n", + "- simple RC circuits\n", + "- simple porous electrode model\n", + "- TLM model (requires matrix calculation). \n", + "\n", + "The parameters are chosen to be simple for ths benchmark.\n", + "\n", + "In each cases, two scenarios are considered:\n", + "- Scenario 1: `f = np.geomspace(1e-3,1e4,100)`, representing the typical range and number of frequencies used in impedance measurements \n", + "- Scenario 2: `f = np.geomspace(1e-5,1e5,int(1e4))`, representing impractical high-resolution measurements.\n", + "\n", + "Lastly, We also compare the speed of calculation with respect to directly calling these functions.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### **Simple RC Circuit**" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "from nleis.nleis_elements_pair import RC\n", + "p = [1,1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Scenario 1" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "direct_time: 0.00043966699740849435,\n", + "eval_time: 0.0029851250001229346,\n", + "graph_time: 0.0008621250162832439\n", + "eval_time/direct_time: 6.789513467506083,\n", + "graph_time/direct_time: 1.960859062346779\n", + "eval_time/graph_time: 3.462519870948968\n" + ] + } + ], + "source": [ + "f = np.geomspace(1e-3,1e4,100)\n", + "def time_dirct_cal():\n", + " return RC(p,f)\n", + "direct_time_RC_s1 = timeit.timeit(time_dirct_cal, number=10)\n", + "\n", + "eval_circuit = NLEISCustomCircuit(\n", + " 'RC', initial_guess=p, graph=False)\n", + "graph_circuit = NLEISCustomCircuit(\n", + " 'RC', initial_guess=p, graph=True)\n", + "\n", + "def time_predict_eval():\n", + " return eval_circuit.predict(f, max_f=np.inf)\n", + "\n", + "def time_predict_graph():\n", + " return graph_circuit.predict(f, max_f=np.inf)\n", + "\n", + "eval_time_RC_s1 = timeit.timeit(time_predict_eval, number=10)\n", + "graph_time_RC_s1 = timeit.timeit(time_predict_graph, number=10)\n", + "\n", + "print(f'direct_time: {direct_time_RC_s1},\\neval_time: {eval_time_RC_s1},\\ngraph_time: {graph_time_RC_s1}')\n", + "print(f'eval_time/direct_time: {eval_time_RC_s1/direct_time_RC_s1},\\ngraph_time/direct_time: {graph_time_RC_s1/direct_time_RC_s1}')\n", + "print(f'eval_time/graph_time: {eval_time_RC_s1/graph_time_RC_s1}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Scenario 2" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "direct_time: 0.01712454200605862,\n", + "eval_time: 0.24936924999929033,\n", + "graph_time: 0.016025250020902604\n", + "eval_time/direct_time: 14.562097480391834,\n", + "graph_time/direct_time: 0.9358060504761477\n", + "eval_time/graph_time: 15.56102086856832\n" + ] + } + ], + "source": [ + "f = np.geomspace(1e-5,1e5,int(1e4))\n", + "def time_dirct_cal():\n", + " return RC(p,f)\n", + "direct_time_RC_s2 = timeit.timeit(time_dirct_cal, number=10)\n", + "\n", + "eval_circuit = NLEISCustomCircuit(\n", + " 'RC', initial_guess=p, graph=False)\n", + "graph_circuit = NLEISCustomCircuit(\n", + " 'RC', initial_guess=p, graph=True)\n", + "\n", + "def time_predict_eval():\n", + " return eval_circuit.predict(f, max_f=np.inf)\n", + "\n", + "def time_predict_graph():\n", + " return graph_circuit.predict(f, max_f=np.inf)\n", + "\n", + "eval_time_RC_s2 = timeit.timeit(time_predict_eval, number=10)\n", + "graph_time_RC_s2 = timeit.timeit(time_predict_graph, number=10)\n", + "\n", + "print(f'direct_time: {direct_time_RC_s2},\\neval_time: {eval_time_RC_s2},\\ngraph_time: {graph_time_RC_s2}')\n", + "print(f'eval_time/direct_time: {eval_time_RC_s2/direct_time_RC_s2},\\ngraph_time/direct_time: {graph_time_RC_s2/direct_time_RC_s2}')\n", + "print(f'eval_time/graph_time: {eval_time_RC_s2/graph_time_RC_s2}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### **Porous Electrode**" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "from nleis.nleis_elements_pair import TP\n", + "p = [1,1,1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Scenario 1" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "direct_time: 0.0007940419891383499,\n", + "eval_time: 0.002853000012692064,\n", + "graph_time: 0.0008623329922556877\n", + "eval_time/direct_time: 3.5930089991688985,\n", + "graph_time/direct_time: 1.0860042718791778\n", + "eval_time/graph_time: 3.3084667272548582\n" + ] + } + ], + "source": [ + "f = np.geomspace(1e-3,1e4,100)\n", + "\n", + "def time_dirct_cal():\n", + " return TP(p,f)\n", + "direct_time_P_s1 = timeit.timeit(time_dirct_cal, number=10)\n", + "\n", + "eval_circuit = NLEISCustomCircuit(\n", + " 'TP', initial_guess=p, graph=False)\n", + "graph_circuit = NLEISCustomCircuit(\n", + " 'TP', initial_guess=p, graph=True)\n", + "\n", + "def time_predict_eval():\n", + " return eval_circuit.predict(f, max_f=np.inf)\n", + "\n", + "def time_predict_graph():\n", + " return graph_circuit.predict(f, max_f=np.inf)\n", + "\n", + "eval_time_P_s1 = timeit.timeit(time_predict_eval, number=10)\n", + "graph_time_P_s1 = timeit.timeit(time_predict_graph, number=10)\n", + "\n", + "print(f'direct_time: {direct_time_P_s1},\\neval_time: {eval_time_P_s1},\\ngraph_time: {graph_time_P_s1}')\n", + "print(f'eval_time/direct_time: {eval_time_P_s1/direct_time_P_s1},\\ngraph_time/direct_time: {graph_time_P_s1/direct_time_P_s1}')\n", + "print(f'eval_time/graph_time: {eval_time_P_s1/graph_time_P_s1}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Scenario 2" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "direct_time: 0.018269459018483758,\n", + "eval_time: 0.23237816599430516,\n", + "graph_time: 0.01989183301338926\n", + "eval_time/direct_time: 12.71948806799376,\n", + "graph_time/direct_time: 1.08880251972782\n", + "eval_time/graph_time: 11.682089118578999\n" + ] + } + ], + "source": [ + "f = np.geomspace(1e-5,1e5,int(1e4))\n", + "\n", + "def time_dirct_cal():\n", + " return TP(p,f)\n", + "direct_time_P_s2 = timeit.timeit(time_dirct_cal, number=10)\n", + "\n", + "eval_circuit = NLEISCustomCircuit(\n", + " 'TP', initial_guess=p, graph=False)\n", + "graph_circuit = NLEISCustomCircuit(\n", + " 'TP', initial_guess=p, graph=True)\n", + "\n", + "def time_predict_eval():\n", + " return eval_circuit.predict(f, max_f=np.inf)\n", + "\n", + "def time_predict_graph():\n", + " return graph_circuit.predict(f, max_f=np.inf)\n", + "\n", + "eval_time_P_s2 = timeit.timeit(time_predict_eval, number=10)\n", + "graph_time_P_s2 = timeit.timeit(time_predict_graph, number=10)\n", + "\n", + "print(f'direct_time: {direct_time_P_s2},\\neval_time: {eval_time_P_s2},\\ngraph_time: {graph_time_P_s2}')\n", + "print(f'eval_time/direct_time: {eval_time_P_s2/direct_time_P_s2},\\ngraph_time/direct_time: {graph_time_P_s2/direct_time_P_s2}')\n", + "print(f'eval_time/graph_time: {eval_time_P_s2/graph_time_P_s2}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### **Transmission Line Model**" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "from nleis.nleis_elements_pair import TLMn\n", + "p = [1,1,1,1,1,10,0.1,0.2]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Scenario 1" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "direct_time: 0.12619937499403022,\n", + "eval_time: 0.16212762499344535,\n", + "graph_time: 0.13549716700799763\n", + "eval_time/direct_time: 1.2846943576472918,\n", + "graph_time/direct_time: 1.0736754204559826\n", + "eval_time/graph_time: 1.1965388544535096\n" + ] + } + ], + "source": [ + "f = np.geomspace(1e-3,1e4,100)\n", + "def time_dirct_cal():\n", + " return TLMn(p,f)\n", + "direct_time_TLM_s1 = timeit.timeit(time_dirct_cal, number=10)\n", + "\n", + "eval_circuit = NLEISCustomCircuit(\n", + " 'TLMn', initial_guess=p, graph=False)\n", + "graph_circuit = NLEISCustomCircuit(\n", + " 'TLMn', initial_guess=p, graph=True)\n", + "\n", + "def time_predict_eval():\n", + " return eval_circuit.predict(f, max_f=np.inf)\n", + "\n", + "def time_predict_graph():\n", + " return graph_circuit.predict(f, max_f=np.inf)\n", + "\n", + "eval_time_TLM_s1 = timeit.timeit(time_predict_eval, number=10)\n", + "graph_time_TLM_s1 = timeit.timeit(time_predict_graph, number=10)\n", + "\n", + "print(f'direct_time: {direct_time_TLM_s1},\\neval_time: {eval_time_TLM_s1},\\ngraph_time: {graph_time_TLM_s1}')\n", + "print(f'eval_time/direct_time: {eval_time_TLM_s1/direct_time_TLM_s1},\\ngraph_time/direct_time: {graph_time_TLM_s1/direct_time_TLM_s1}')\n", + "print(f'eval_time/graph_time: {eval_time_TLM_s1/graph_time_TLM_s1}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Scenario 2" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "direct_time: 12.843668707995676,\n", + "eval_time: 17.604674500005785,\n", + "graph_time: 15.474227416008944\n", + "eval_time/direct_time: 1.3706889285493795,\n", + "graph_time/direct_time: 1.2048136531562548\n", + "eval_time/graph_time: 1.1376771212366168\n" + ] + } + ], + "source": [ + "f = np.geomspace(1e-5,1e5,int(1e4))\n", + "def time_dirct_cal():\n", + " return TLMn(p,f)\n", + "direct_time_TLM_s2 = timeit.timeit(time_dirct_cal, number=10)\n", + "\n", + "eval_circuit = NLEISCustomCircuit(\n", + " 'TLMn', initial_guess=p, graph=False)\n", + "graph_circuit = NLEISCustomCircuit(\n", + " 'TLMn', initial_guess=p, graph=True)\n", + "\n", + "def time_predict_eval():\n", + " return eval_circuit.predict(f, max_f=np.inf)\n", + "\n", + "def time_predict_graph():\n", + " return graph_circuit.predict(f, max_f=np.inf)\n", + "\n", + "eval_time_TLM_s2 = timeit.timeit(time_predict_eval, number=10)\n", + "graph_time_TLM_s2 = timeit.timeit(time_predict_graph, number=10)\n", + "\n", + "print(f'direct_time: {direct_time_TLM_s2},\\neval_time: {eval_time_TLM_s2},\\ngraph_time: {graph_time_TLM_s2}')\n", + "print(f'eval_time/direct_time: {eval_time_TLM_s2/direct_time_TLM_s2},\\ngraph_time/direct_time: {graph_time_TLM_s2/direct_time_TLM_s2}')\n", + "print(f'eval_time/graph_time: {eval_time_TLM_s2/graph_time_TLM_s2}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### **Summary**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following table provides a summary of the benchmark results.\n", + "\n", + "It is evident that graph-based evaluation can achieve at least a **3X speed-up** in typical calculations involving NumPy-oriented circuit functions such as `RC` and `TP`. \n", + "\n", + "However, for matrix-oriented circuit functions like `TLMn`, graph-based evaluation shows limited performance gains, likely due to the matrix-solving operations dominating the computation time in such cases" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Condition Direct Cal [s] eval() [s] graph [s] eval/direct graph/direct eval/graph\n", + "----------- ---------------- ------------ ------------ ------------- -------------- ------------\n", + "RC S1 0.000439667 0.00298513 0.000862125 6.78951 1.96086 3.46252\n", + "RC S2 0.0171245 0.249369 0.0160253 14.5621 0.935806 15.561\n", + "P S1 0.000794042 0.002853 0.000862333 3.59301 1.086 3.30847\n", + "P S2 0.0182695 0.232378 0.0198918 12.7195 1.0888 11.6821\n", + "TLM S1 0.126199 0.162128 0.135497 1.28469 1.07368 1.19654\n", + "TLM S2 12.8437 17.6047 15.4742 1.37069 1.20481 1.13768\n" + ] + } + ], + "source": [ + "table = [[\"RC S1\",direct_time_RC_s1,eval_time_RC_s1,graph_time_RC_s1,eval_time_RC_s1/direct_time_RC_s1,graph_time_RC_s1/direct_time_RC_s1,eval_time_RC_s1/graph_time_RC_s1],\n", + " [\"RC S2\",direct_time_RC_s2,eval_time_RC_s2,graph_time_RC_s2,eval_time_RC_s2/direct_time_RC_s2,graph_time_RC_s2/direct_time_RC_s2,eval_time_RC_s2/graph_time_RC_s2],\n", + " [\"P S1\",direct_time_P_s1,eval_time_P_s1,graph_time_P_s1,eval_time_P_s1/direct_time_P_s1,graph_time_P_s1/direct_time_P_s1,eval_time_P_s1/graph_time_P_s1],\n", + " [\"P S2\",direct_time_P_s2,eval_time_P_s2,graph_time_P_s2,eval_time_P_s2/direct_time_P_s2,graph_time_P_s2/direct_time_P_s2,eval_time_P_s2/graph_time_P_s2],\n", + " [\"TLM S1\",direct_time_TLM_s1,eval_time_TLM_s1,graph_time_TLM_s1,eval_time_TLM_s1/direct_time_TLM_s1,graph_time_TLM_s1/direct_time_TLM_s1,eval_time_TLM_s1/graph_time_TLM_s1],\n", + " [\"TLM S2\",direct_time_TLM_s2,eval_time_TLM_s2,graph_time_TLM_s2,eval_time_TLM_s2/direct_time_TLM_s2,graph_time_TLM_s2/direct_time_TLM_s2,eval_time_TLM_s2/graph_time_TLM_s2],]\n", + "print(tabulate(table,headers=[\"Condition\",\"Direct Cal [s]\", \"eval() [s]\", 'graph [s]', 'eval/direct', 'graph/direct', 'eval/graph']))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "nleis", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.10" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/index.rst b/docs/source/index.rst index 5d928e4..7254b7d 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -14,7 +14,8 @@ Welcome to :code:`nleis.py`'s documentation! nleis_fitting visualization data-processing - + release-notes + Indices and tables ================== diff --git a/docs/source/release-notes.rst b/docs/source/release-notes.rst new file mode 100644 index 0000000..cf8f671 --- /dev/null +++ b/docs/source/release-notes.rst @@ -0,0 +1,45 @@ +==================== +Release Notes +==================== + +.. Version 0.2 +.. --------------------------- + + +Version 0.1.1 (2025-01-06) +--------------------------- +This is the official release for the JOSS paper. + +**What's Changed** + +- Documentation updates by @dt-schwartz and @yuefan98 +- Bug fixes by @yuefan98 in https://github.com/yuefan98/nleis.py/pull/25 + +**Full Changelog**: https://github.com/yuefan98/nleis.py/compare/v0.1...v0.1.1 + +Version 0.1 (2024-09-26) +------------------------- +We are excited to announce the first official release of nleis.py! This release marks a significant step forward for nonlinear impedance analysis and will be submitted to JOSS for peer review. + +**Key Features:** + +- Simultaneous fitting and analysis of EIS and 2nd-NLEIS. +- Full support for nonlinear equivalent circuit (nECM) modeling and analysis. +- Various linear and nonlinear circuit element pairs derived from existing literature. +- Seamless integration with impedance.py for expanded impedance analysis capabilities. + +**Improvements:** + +- Comprehensive [documentation](https://nleispy.readthedocs.io/en/latest/), including a Getting Started guide and API reference. +- Improved documentation for supported circuit elements. +- Improved code handling for better performance and readability. + +**Bug Fixes** + +- Initial testing and issue resolution to ensure smooth functionality. + +**New Contributors** + +- Special thanks to @mdmurbach for joining the team and enhancing the package quality. + +**Full Changelog**: https://github.com/yuefan98/nleis.py/commits/v0.1 \ No newline at end of file diff --git a/nleis/fitting.py b/nleis/fitting.py index 43cb0a6..5f2708f 100644 --- a/nleis/fitting.py +++ b/nleis/fitting.py @@ -613,6 +613,9 @@ def extract_circuit_elements(circuit): class CircuitGraph: + ''' + A class to represent a circuit as a directed graph. + ''' # regular expression to find parallel and difference blocks _parallel_difference_block_expression = re.compile(r'(?:p|d)\([^()]*\)') @@ -620,6 +623,8 @@ class CircuitGraph: _whitespce = re.compile(r"\s+") def __init__(self, circuit, constants=None): + ''' + Initialize the CircuitGraph object.''' # remove all whitespace from the circuit string self.circuit = self._whitespce.sub("", circuit) # parse the circuit string and initialize the graph @@ -630,6 +635,9 @@ def __init__(self, circuit, constants=None): self.constants = constants if constants is not None else dict() def parse_circuit(self): + ''' + Parse the circuit string and initialize the graph. + ''' # initialize the node counters for each type of block self.snum = 1 self.pnum = 1 @@ -688,6 +696,9 @@ def parse_circuit(self): # function to add series elements to the graph def add_series_elements(self, elem): + ''' + Add series elements to the graph. + ''' selem = elem.split("-") if len(selem) > 1: node = f"s{self.snum}" @@ -702,11 +713,16 @@ def add_series_elements(self, elem): # function to visualize the graph def visualize_graph(self, **kwargs): + ''' + Visualize the graph.''' pos = nx.multipartite_layout(self.graph, subset_key="layer") nx.draw_networkx(self.graph, pos=pos, **kwargs) # function to compute the impedance of the circuit def compute(self, f, *parameters): + ''' + Compute the impedance of the circuit at the given frequencies. + ''' node_results = {} pindex = 0 for node in self.execution_order: @@ -733,6 +749,9 @@ def compute(self, f, *parameters): # To enable comparision def __eq__(self, other): + ''' + Compare two CircuitGraph objects for equality. + ''' if not isinstance(other, CircuitGraph): return False # Compare the internal graph attributes @@ -742,14 +761,25 @@ def __eq__(self, other): # To enable direct calling def __call__(self, f, *parameters): + ''' + Compute the impedance of the circuit at the given frequencies. + And convert it to a long array for curve_fit. + ''' Z = self.compute(f, *parameters) return np.hstack([Z.real, Z.imag]) def compute_long(self, f, *parameters): + ''' + Compute the impedance of the circuit at the given frequencies. + And convert it to a long array for curve_fit. + ''' Z = self.compute(f, *parameters) return np.hstack([Z.real, Z.imag]) def calculate_circuit_length(self): + ''' + calculate the number of parameters in the circuit + ''' n_params = [ getattr(Zfunc, "num_params", 0) for node, Zfunc in self.graph.nodes(data="Z") @@ -758,4 +788,7 @@ def calculate_circuit_length(self): def format_parameter_name(name, j, n_params): + ''' + Format the parameter name for the given element. + ''' return f"{name}_{j}" if n_params > 1 else f"{name}" diff --git a/nleis/nleis_elements_pair.py b/nleis/nleis_elements_pair.py index fb6503d..d01870b 100644 --- a/nleis/nleis_elements_pair.py +++ b/nleis/nleis_elements_pair.py @@ -1079,6 +1079,7 @@ def TDCn(p, f): def A_matrices_TLMn(N, Rpore, Z12t): """ Construct the matrix `Ax` for the TLMn model + Parameters ---------- N : int @@ -1087,10 +1088,12 @@ def A_matrices_TLMn(N, Rpore, Z12t): Pore electrolyte resistance Z12t : np.complex128 The single element impedance at 2ω + Returns ------- Ax : np.ndarray The matrix `Ax` for the TLMn model + """ Ax = np.zeros((N, N), dtype=np.complex128) From 669c3ae0e17b0599286e78d3c63cea0beb6d2e14 Mon Sep 17 00:00:00 2001 From: Yuefan Ji Date: Mon, 20 Jan 2025 13:18:19 -0800 Subject: [PATCH 6/6] Code improvement: * simplify the code for parse_circuit * Update compute methods to return impedance value * add networkx dependency to environment.yml --- environment.yml | 1 + nleis/fitting.py | 28 +++++++++++++--------------- nleis/nleis_tests/test_graph.py | 2 +- 3 files changed, 15 insertions(+), 16 deletions(-) diff --git a/environment.yml b/environment.yml index 7466603..409e5f3 100644 --- a/environment.yml +++ b/environment.yml @@ -16,5 +16,6 @@ dependencies: - tqdm - tabulate - impedance + - networkx diff --git a/nleis/fitting.py b/nleis/fitting.py index 5f2708f..73e1047 100644 --- a/nleis/fitting.py +++ b/nleis/fitting.py @@ -338,7 +338,7 @@ def opt_function(x): np.hstack([Z.real, Z.imag])) def opt_function_graph(x): - return rmse(cg(f, *x), np.hstack([Z.real, Z.imag])) + return rmse(cg.compute_long(f, *x), np.hstack([Z.real, Z.imag])) class BasinhoppingBounds(object): """ Adapted from the basinhopping documetation @@ -667,22 +667,21 @@ def parse_circuit(self): for pd in pd_blocks: operator = pd[0] pd_elem = pd[2:-1].split(",") + if operator == "p": - pnode = f"p{self.pnum}" + nnum = self.pnum self.pnum += 1 - self.graph.add_node(pnode, Z=circuit_elements["p"]) - for elem in pd_elem: - elem = self.add_series_elements(elem) - self.graph.add_edge(elem, pnode) - parsing_circuit = parsing_circuit.replace(pd, pnode) elif operator == "d": - dnode = f"d{self.dnum}" + nnum = self.dnum self.dnum += 1 - self.graph.add_node(dnode, Z=circuit_elements["d"]) - for elem in pd_elem: - elem = self.add_series_elements(elem) - self.graph.add_edge(elem, dnode) - parsing_circuit = parsing_circuit.replace(pd, dnode) + + node = f"{operator}{nnum}" + self.graph.add_node(node, Z=circuit_elements[operator]) + for elem in pd_elem: + elem = self.add_series_elements(elem) + self.graph.add_edge(elem, node) + parsing_circuit = parsing_circuit.replace(pd, node) + pd_blocks = self._parallel_difference_block_expression.findall( parsing_circuit) @@ -763,10 +762,9 @@ def __eq__(self, other): def __call__(self, f, *parameters): ''' Compute the impedance of the circuit at the given frequencies. - And convert it to a long array for curve_fit. ''' Z = self.compute(f, *parameters) - return np.hstack([Z.real, Z.imag]) + return Z def compute_long(self, f, *parameters): ''' diff --git a/nleis/nleis_tests/test_graph.py b/nleis/nleis_tests/test_graph.py index ce2b5a3..1ecd552 100644 --- a/nleis/nleis_tests/test_graph.py +++ b/nleis/nleis_tests/test_graph.py @@ -129,7 +129,7 @@ def test_CircuitGraph(): assert cg.calculate_circuit_length() == 4 # test for __call__ - assert np.allclose(cg.compute_long( + assert np.allclose(cg.compute( frequencies, *params), cg(frequencies, *params)) # test for __eq__