From c63993b93821be642f453b6b60a520d0d2552f05 Mon Sep 17 00:00:00 2001 From: vyrjana <41105805+vyrjana@users.noreply.github.com> Date: Tue, 1 Oct 2024 22:45:01 +0300 Subject: [PATCH] Merged commits from various branches --- CHANGELOG.md | 33 + .../LICENSE-DRT-from-Loewner-framework.txt | 22 + LICENSES/LICENSE-pyDRTtools.txt | 2 +- LICENSES/README.md | 5 + docs/source/apidocs_drt.rst | 20 +- docs/source/apidocs_fitting.rst | 4 +- docs/source/guide_drt.rst | 215 +++++- docs/source/guide_fitting.rst | 325 +++++++- docs/source/guide_installing.rst | 24 +- docs/source/guide_kramers_kronig.rst | 27 +- docs/source/substitutions.rst | 40 +- src/pyimpspec/__init__.py | 4 + src/pyimpspec/analysis/__init__.py | 7 +- src/pyimpspec/analysis/drt/__init__.py | 20 +- src/pyimpspec/analysis/drt/bht.py | 60 +- src/pyimpspec/analysis/drt/lm.py | 714 ++++++++++++++++++ src/pyimpspec/analysis/drt/mrq_fit.py | 4 +- src/pyimpspec/analysis/drt/peak_analysis.py | 520 +++++++++++++ src/pyimpspec/analysis/drt/result.py | 78 +- src/pyimpspec/analysis/drt/tr_nnls.py | 13 +- src/pyimpspec/analysis/drt/tr_rbf.py | 21 +- src/pyimpspec/analysis/fitting.py | 258 +++++-- .../algorithms/representation.py | 4 +- src/pyimpspec/analysis/kramers_kronig/cnls.py | 14 +- .../analysis/kramers_kronig/exploratory.py | 8 +- .../analysis/kramers_kronig/result.py | 4 +- .../analysis/kramers_kronig/single.py | 4 +- src/pyimpspec/analysis/zhit/__init__.py | 8 +- src/pyimpspec/circuit/__init__.py | 8 +- src/pyimpspec/circuit/circuit_builder.py | 8 +- src/pyimpspec/circuit/diagrams/circuitikz.py | 2 +- src/pyimpspec/circuit/diagrams/schemdraw.py | 2 +- src/pyimpspec/circuit/kramers_kronig.py | 2 + src/pyimpspec/circuit/parser.py | 2 +- src/pyimpspec/circuit/registry.py | 36 +- src/pyimpspec/cli/args.py | 62 +- src/pyimpspec/cli/drt.py | 39 +- src/pyimpspec/cli/utility.py | 2 +- src/pyimpspec/data/__init__.py | 2 +- src/pyimpspec/data/data_set.py | 14 +- src/pyimpspec/mock_data.py | 76 +- src/pyimpspec/plot/mpl/drt.py | 5 + src/pyimpspec/plot/mpl/gamma.py | 314 +++++++- src/pyimpspec/plot/mpl/imaginary.py | 8 + src/pyimpspec/plot/mpl/kramers_kronig.py | 2 +- src/pyimpspec/plot/mpl/magnitude.py | 5 + src/pyimpspec/plot/mpl/nyquist.py | 9 + src/pyimpspec/plot/mpl/phase.py | 9 + src/pyimpspec/plot/mpl/real.py | 8 + tests/test_cli.py | 33 + tests/test_drt.py | 344 ++++++++- tests/test_element.py | 18 +- tests/test_fitting.py | 121 +++ 53 files changed, 3352 insertions(+), 237 deletions(-) create mode 100644 LICENSES/LICENSE-DRT-from-Loewner-framework.txt create mode 100644 src/pyimpspec/analysis/drt/lm.py create mode 100644 src/pyimpspec/analysis/drt/peak_analysis.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 68edf07..3c05b79 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,36 @@ +# 5.1.0 (2024/MM/DD) + +- Added support for analyzing the peaks in DRT results by fitting skew normal distributions: + - Added an `analyze_peaks` method to the `DRTResult` class. + - Added a `DRTPeaks` class, which can be used to, e.g., calculate the area of a peak. +- Added an implementation of the Loewner method for calculating the distribution of relaxation times: + - Added a `calculate_drt_lm` function. + - Added an `LMResult` class. + - Added `--model-order` and `--model-order-method` CLI arguments. + - Updated the `plot.mpl.plot_gamma` function to support plotting `LMResult` instances. +- Added the ability to define constraints when fitting circuits: + - Added a `generate_fit_identifiers` function. + - Added a `FitIdentifiers` class. +- Added support for plotting DRT results as gamma vs f. + - Added CLI argument for plotting DRT results as gamma vs f. +- Added circuits for generating mock data (`CIRCUIT_13`, `CIRCUIT_13_INVALID`, `CIRCUIT_14`, and `CIRCUIT_14_INVALID`). +- Updated the TR-RBF implementation to be based off of a newer version of pyDRTtools. +- Updated plotting functions to support drawing smoother lines if the input result has a `circuit` property. +- Updated an exception message to provide more information if a variable that is being fitted is out of bounds. +- Updated documentation. +- Updated the `generate_mock_data` function so that it attempts to cast arguments to the appropriate type if some other types of values (e.g., integers) are provided instead. +- Updated the `circuit.registry.register_element` function to support a new `private` keyword argument. + - Updated the `circuit.registry.get_elements` function to not include by default `Element` classes that were registered with `private=True`. + - Updated the `KramersKronigRC` and `KramersKronigAdmittanceRC` classes to be registered with `private=True`. +- Fixed a bug that caused methods such as `DRTResult.get_peaks` to miss peaks at the time constant extremes. +- Fixed a bug that caused an element's parameters in `FitResult.to_parameters_dataframe` to not be in a sorted order. +- Fixed the previously unimplemented `FitResult.get_parameters` method. +- Fixed a bug that caused `FitResult.to_parameters_dataframe` to return negative standard errors when the fitted value was negative. +- Fixed a bug that could cause `FitResult.to_parameters_dataframe` to divide by zero. +- Fixed a bug where an exception would be raised when whitespace was included between keyword arguments when passing circuit identifiers or CDCs via the CLI (e.g., `` would previously raise an exception whereas `` would not). +- Refactored some of the code. + + # 5.0.2 (2024/09/10) - Updated the documentation (e.g., various function and method docstrings have been updated). Notably, the functions related to handling progress updates are now included in the API documentation. diff --git a/LICENSES/LICENSE-DRT-from-Loewner-framework.txt b/LICENSES/LICENSE-DRT-from-Loewner-framework.txt new file mode 100644 index 0000000..d339679 --- /dev/null +++ b/LICENSES/LICENSE-DRT-from-Loewner-framework.txt @@ -0,0 +1,22 @@ +MIT License + +Copyright (c) 2023 projectsEECandDRI + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + diff --git a/LICENSES/LICENSE-pyDRTtools.txt b/LICENSES/LICENSE-pyDRTtools.txt index aaf4f60..c01535f 100644 --- a/LICENSES/LICENSE-pyDRTtools.txt +++ b/LICENSES/LICENSE-pyDRTtools.txt @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2020 ciuccislab +Copyright (c) 2023 ciuccislab Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/LICENSES/README.md b/LICENSES/README.md index 9404026..b8cb91a 100644 --- a/LICENSES/README.md +++ b/LICENSES/README.md @@ -8,6 +8,11 @@ - License: GPLv3 or later - Dependency. +# DRT-from-Loewner-framework +- https://github.com/projectsEECandDRI/DRT-from-Loewner-framework +- License: MIT +- Ported to Python and modified. + # DRT-python-code - https://github.com/akulikovsky/DRT-python-code - License: GPLv3 or later diff --git a/docs/source/apidocs_drt.rst b/docs/source/apidocs_drt.rst index fbca888..f121351 100644 --- a/docs/source/apidocs_drt.rst +++ b/docs/source/apidocs_drt.rst @@ -8,7 +8,16 @@ A collection of functions and classes for calculating the distribution of relaxa Wrapper function and base class ------------------------------- .. automodule:: pyimpspec - :members: calculate_drt, DRTResult + :members: calculate_drt + +.. automodule:: pyimpspec + :members: DRTResult + + +Classes related to peak analysis +-------------------------------- +.. automodule:: pyimpspec + :members: DRTPeaks, DRTPeak Method functions and classes @@ -23,6 +32,15 @@ BHT method :members: BHTResult +LM method +~~~~~~~~~~ +.. automodule:: pyimpspec.analysis.drt + :members: calculate_drt_lm + +.. automodule:: pyimpspec.analysis.drt + :members: LMResult + + m(RQ)fit method ~~~~~~~~~~~~~~~ .. automodule:: pyimpspec.analysis.drt diff --git a/docs/source/apidocs_fitting.rst b/docs/source/apidocs_fitting.rst index c322a74..898ee16 100644 --- a/docs/source/apidocs_fitting.rst +++ b/docs/source/apidocs_fitting.rst @@ -9,12 +9,12 @@ A collection of functions and classes for fitting equivalent circuits to data se Function -------- .. automodule:: pyimpspec - :members: fit_circuit + :members: fit_circuit, generate_fit_identifiers Classes ------- .. automodule:: pyimpspec - :members: FitResult, FittedParameter + :members: FitResult, FittedParameter, FitIdentifiers .. raw:: latex diff --git a/docs/source/guide_drt.rst b/docs/source/guide_drt.rst index 8721463..c9c66ba 100644 --- a/docs/source/guide_drt.rst +++ b/docs/source/guide_drt.rst @@ -19,6 +19,7 @@ The supported methods Implementations based on the following approaches are included in pyimpspec: - **BHT**: The Bayesian Hilbert transform method (see `Liu et al. (2020) `_) that was originally implemented in `DRTtools `_ and `pyDRTtools `_. +- **LM**: The Loewner method (see `Mayo and Antoulas (2007) `_, `Sorrentino et al. (2022) `_, `Rüther et al. (2023) `_, and `Sorrentino et al. (2023) `_) that was implemented in `DRT-from-Loewner-framework `_. - **m(RQ)fit**: The multi-(RQ)-fit method (see `Boukamp (2015) `_ and `Boukamp and Rolle (2017) `_). - **TR-NNLS**: Tikhonov regularization and non-negative least squares (see `Kulikovsky (2021) `_) that was originally implemented in `DRT-python-code `_. - **TR-RBF**: Tikhonov regularization and radial basis function (or piecewise linear) discretization (see `Wan et al. (2015) `_, `Ciucci and Chen (2015) `_, `Effat and Ciucci (2017) `_, and `Maradesa et al. (2023) `_) that was originally implemented in DRTtools_ and pyDRTtools_. @@ -51,10 +52,12 @@ Each method has its own function that can be used but there is also a wrapper fu ... ) >>> from pyimpspec.analysis.drt import ( ... BHTResult, # Result of the BHT method + ... LMResult, # Result of the Loewner method ... MRQFitResult, # Result of the m(RQ)fit method ... TRNNLSResult, # Result of the TR-NNLS method ... TRRBFResult, # Result of the TR-RBF method ... calculate_drt_bht, # BHT method + ... calculate_drt_lm, # Loewner method ... calculate_drt_mrq_fit, # m(RQ)fit method ... calculate_drt_tr_nnls, # TR-NNLS method ... calculate_drt_tr_rbf, # TR-RBF method @@ -76,6 +79,7 @@ Below are some figures demonstrating the results of the methods listed above whe ) from pyimpspec.analysis.drt import ( calculate_drt_bht, + calculate_drt_lm, calculate_drt_mrq_fit, calculate_drt_tr_nnls, calculate_drt_tr_rbf, @@ -87,11 +91,18 @@ Below are some figures demonstrating the results of the methods listed above whe ax.set_ylim(-100, 900) data = generate_mock_data("CIRCUIT_1", noise=5e-2, seed=42)[0] + figure, axes = mpl.plot_nyquist(data, colors=dict(impedance="black"), markers=dict(impedance="o")) + figure.tight_layout() + drt = calculate_drt_bht(data) figure, axes = mpl.plot_gamma(drt) adjust_limits(axes[0]) figure.tight_layout() + drt = calculate_drt_lm(data) + figure, axes = mpl.plot_gamma(drt) + figure.tight_layout() + circuit = parse_cdc("R(RQ)(RQ)") fit = fit_circuit(circuit, data) drt = calculate_drt_mrq_fit(data, fit.circuit, fit=fit) @@ -115,16 +126,216 @@ Below are some figures demonstrating the results of the methods listed above whe \clearpage +.. note:: + + The Loewner method can provide an estimate for the high-frequency resistance, but the high-frequency (i.e., low time constant) resistive-capacitive peak may not always provide this value directly. + + +.. plot:: + + from pyimpspec import ( + generate_mock_data, + parse_cdc, + ) + from pyimpspec.analysis.drt import calculate_drt_lm + from pyimpspec import mpl + + cdc = "R{R=140}(R{R=230}C{C=1e-6})(R{R=576}C{C=1e-4})(R{R=150}L{L=4e1})" + circuit = parse_cdc(cdc) + drawing = circuit.to_drawing() + drawing.draw() + + data = generate_mock_data(cdc, noise=0)[0] + drt = calculate_drt_lm(data) + + figure, axes = mpl.plot_nyquist(data, colors=dict(impedance="black"), markers=dict(impedance="o")) + mpl.plot_nyquist(drt, colors=dict(impedance="red"), markers=dict(impedance="+"), figure=figure, axes=axes) + figure.tight_layout() + + figure, axes = mpl.plot_gamma(drt) + figure.tight_layout() + + +.. code:: + + | tau, RC (s) | gamma, RC (ohm) | tau, RL (s) | gamma, RL (ohm) | + |--------------:|------------------:|--------------:|------------------:| + | 8.00717e-12 | 289.999 | 0.266667 | 150.012 | + | 0.00023 | 230 | nan | nan | + | 0.0576 | 576.007 | nan | nan | + + +This example shows that some additional processing of the obtained values may be necessary. +In this case, the high-frequency resistive-capacitive peak does not directly tell us the value of :math:`R_1`. +However, we can still obtain an estimated value: :math:`R_1 \approx 289.999\ \Omega - 150.012\ \Omega \approx 140\ \Omega`. +Other estimated values such as capacitances and inductances can also be obtained: + +- :math:`R_2 \approx 230\ \Omega` and :math:`C_1 \approx 0.0023\ {\rm s} / 230\ \Omega \approx 10^{-6}\ {\rm F}` +- :math:`R_3 \approx 576\ \Omega` and :math:`C_2 \approx 0.0576\ {\rm s} / 576.007\ \Omega \approx 10^{-4}\ {\rm F}` +- :math:`R_4 \approx 150\ \Omega` and :math:`L_2 \approx 0.266667\ {\rm s} \times 150.012\ \Omega \approx 40\ {\rm H}` + +The above is a best-case scenario, but even with a bit of added noise (:math:`\sigma = 0.05\ \% \times |Z|`) one can still obtain decent estimates despite the additional peaks: + +- :math:`R1 \approx 289.925\ \Omega - 149.085\ \Omega \approx 141\ \Omega`. +- :math:`R_2 \approx 230\ \Omega` and :math:`C_1 \approx 0.00230559\ {\rm s} / 230.233\ \Omega \approx 10^{-6}\ {\rm F}` +- :math:`R_3 \approx 574\ \Omega` and :math:`C_2 \approx 0.0574022\ {\rm s} / 574.438\ \Omega \approx 10^{-4}\ {\rm F}` +- :math:`R_4 \approx 149\ \Omega` and :math:`L_2 \approx 0.271608\ {\rm s} \times 149.085\ \Omega \approx 40\ {\rm H}` + + +.. plot:: + + from pyimpspec import ( + generate_mock_data, + parse_cdc, + ) + from pyimpspec.analysis.drt import calculate_drt_lm + from pyimpspec import mpl + + cdc = "R{R=140}(R{R=230}C{C=1e-6})(R{R=576}C{C=1e-4})(R{R=150}L{L=4e1})" + data = generate_mock_data(cdc, noise=5e-2, seed=42)[0] + drt = calculate_drt_lm(data) + + figure, axes = mpl.plot_gamma(drt) + figure.tight_layout() + + +.. code:: + + | tau, RC (s) | gamma, RC (ohm) | tau, RL (s) | gamma, RL (ohm) | + |--------------:|------------------:|--------------:|------------------:| + | 4.60645e-10 | 289.925 | 3.05543e-06 | 0.00891355 | + | 5.28003e-06 | 0.0198689 | 3.05543e-06 | 0.00891355 | + | 5.28003e-06 | 0.0198689 | 1.63833e-05 | 0.0126733 | + | 5.92308e-06 | 0.0270496 | 1.63833e-05 | 0.0126733 | + | 5.92308e-06 | 0.0270496 | 2.75326e-05 | 0.0099311 | + | 1.01729e-05 | 0.00475814 | 2.75326e-05 | 0.0099311 | + | 1.01729e-05 | 0.00475814 | 0.000109749 | 0.0314605 | + | 3.70669e-05 | 0.0147584 | 0.000109749 | 0.0314605 | + | 3.70669e-05 | 0.0147584 | 0.000166522 | 0.0195471 | + | 7.28352e-05 | 0.0196474 | 0.000166522 | 0.0195471 | + | 7.28352e-05 | 0.0196474 | 0.000272624 | 0.0393392 | + | 0.000230559 | 230.233 | 0.000272624 | 0.0393392 | + | 0.000771132 | 0.00893735 | 0.000512795 | 0.0399306 | + | 0.000771132 | 0.00893735 | 0.000512795 | 0.0399306 | + | 0.00105728 | 0.00459818 | 0.00789638 | 0.0116901 | + | 0.00105728 | 0.00459818 | 0.00789638 | 0.0116901 | + | 0.00185517 | 0.0331227 | 0.00956206 | 1.02824 | + | 0.00185517 | 0.0331227 | 0.271608 | 149.085 | + | 0.00370024 | 0.0500698 | 0.658958 | 0.0536769 | + | 0.00370024 | 0.0500698 | 0.658958 | 0.0536769 | + | 0.00546967 | 0.044606 | 1.45096 | 0.0433677 | + | 0.00546967 | 0.044606 | 1.45096 | 0.0433677 | + | 0.0126403 | 0.150534 | nan | nan | + | 0.0126403 | 0.150534 | nan | nan | + | 0.0285097 | 0.00708604 | nan | nan | + | 0.0285097 | 0.00708604 | nan | nan | + | 0.0455447 | 0.116258 | nan | nan | + | 0.0455447 | 0.116258 | nan | nan | + | 0.0574022 | 574.438 | nan | nan | + | 0.0675327 | 0.0681987 | nan | nan | + | 0.0675327 | 0.0681987 | nan | nan | + | 0.131324 | 0.0287684 | nan | nan | + | 0.131324 | 0.0287684 | nan | nan | + | 0.228107 | 0.0382958 | nan | nan | + | 0.228107 | 0.0382958 | nan | nan | + | 0.276414 | 0.0018757 | nan | nan | + | 0.276414 | 0.0018757 | nan | nan | + | 1.00071 | 1.19876 | nan | nan | + + +.. raw:: latex + + \clearpage + + +Peak analysis +------------- + +DRT results (aside from those obtained using the Loewner method) can be analyzed using skew normal distributions (see `Danzer (2019) `_ and `Plank et al. (2024) `_). +This facilitates the estimation of the polarization contribution of each peak based on the area of that peak. + +.. plot:: + + from pyimpspec import ( + generate_mock_data, + parse_cdc, + ) + from pyimpspec.analysis.drt import calculate_drt_tr_rbf + from pyimpspec import mpl + + cdc = "R{R=100}(R{R=300}C{C=5e-6})(R{R=450}C{C=1e-5})" + circuit = parse_cdc(cdc) + drawing = circuit.to_drawing() + drawing.draw() + + data = generate_mock_data(cdc, noise=5e-2, seed=42)[0] + drt = calculate_drt_tr_rbf(data) + + figure, axes = mpl.plot_nyquist(data, colors=dict(impedance="black"), markers=dict(impedance="o")) + mpl.plot_nyquist(drt, colors=dict(impedance="red"), markers=dict(impedance="+"), figure=figure, axes=axes) + figure.tight_layout() + + figure, axes = mpl.plot_gamma(drt) + figure.tight_layout() + + peaks = drt.analyze_peaks() + figure, axes = mpl.plot_gamma(drt) + figure.tight_layout() + + peaks.to_peaks_dataframe().to_markdown(index=False) + + +.. code:: + + | tau (s) | gamma (ohm) | R_peak (ohm) | + |------------:|--------------:|---------------:| + | 1.24844e-05 | 2.68972e-08 | 4.33529e-07 | + | 0.00155063 | 520.733 | 280.396 | + | 0.00419888 | 831.637 | 470.508 | + | 9.76195 | 0.364749 | 0.552689 | + + +.. list-table:: True and estimated values of parallel RC elements. + :header-rows: 1 + + * - Component + - True value + - Estimated value + * - :math:`R_2` + - :math:`300\ \Omega` + - :math:`280\ \Omega` + * - :math:`C_1` + - :math:`5.0 \times 10^{-6}\ {\rm F}` + - :math:`5.5 \times 10^{-6}\ {\rm F}` + * - :math:`R_3` + - :math:`450\ \Omega` + - :math:`471\ \Omega` + * - :math:`C_2` + - :math:`1.0 \times 10^{-5}\ {\rm F}` + - :math:`8.9 \times 10^{-6}\ {\rm F}` + + +.. raw:: latex + + \clearpage + + References: -- `Boukamp, B.A., 2015, Electrochim. Acta, 154, 35-46 `_ - `Boukamp, B.A. and Rolle, A, 2017, Solid State Ionics, 302, 12-18 `_ +- `Boukamp, B.A., 2015, Electrochim. Acta, 154, 35-46 `_ - `Ciucci, F. and Chen, C., 2015, Electrochim. Acta, 167, 439-454 `_ +- `Danzer, M.A., 2019, Batteries, 5, 53 `_ - `Effat, M. B. and Ciucci, F., 2017, Electrochim. Acta, 247, 1117-1129 `_ - `Kulikovsky, A., 2021, J. Electrochem. Soc., 168, 044512 `_ - `Liu, J., Wan, T. H., and Ciucci, F., 2020, Electrochim. Acta, 357, 136864 `_ -- `Wan, T. H., Saccoccio, M., Chen, C., and Ciucci, F., 2015, Electrochim. Acta, 184, 483-499 `_ - `Maradesa, A., Py, B., Wan, T.H., Effat, M.B., and Ciucci F., 2023, J. Electrochem. Soc, 170, 030502 `_ +- `Mayo, A.J. and Antoulas, A.C., 2007, Linear Algebra Its Appl. 425, 634–662 `_, +- `Plank, C., Rüther, T., Jahn, L., Schamel, M., Schmidt, J.P., Ciucci, F., and Danzer, M.A., 2024, Journal of Power Sources, 594, 233845 `_ +- `Rüther, T., Gosea, I.V., Jahn, L., Antoulas, A.C., and Danzer, M.A., 2023, Batteries, 9, 132 `_ +- `Sorrentino, A., Patel, B., Gosea, I.V., Antoulas, A.C., and Vidaković-Koch, T., 2022, SSRN `_ +- `Sorrentino, A., Patel, B., Gosea, I.V., Antoulas, A.C., and Vidaković-Koch, T., 2023, J. Power Sources, 585, 223575 `_ +- `Wan, T. H., Saccoccio, M., Chen, C., and Ciucci, F., 2015, Electrochim. Acta, 184, 483-499 `_ .. raw:: latex diff --git a/docs/source/guide_fitting.rst b/docs/source/guide_fitting.rst index 1f53b98..2cb2c8b 100644 --- a/docs/source/guide_fitting.rst +++ b/docs/source/guide_fitting.rst @@ -14,7 +14,7 @@ Each element in the circuit has a known expression for its impedance and so does Performing a fit ------------------- +---------------- The |fit_circuit| function performs the fitting and returns a |FitResult| object. @@ -190,6 +190,329 @@ The contents of ``parameters`` and ``statistics`` in the example above would be | W_4 | Y | 0.000400664 | 0.0303443 | S*s^(1/2) | No | + +Using mathematical constraints +------------------------------ + +It is possible to make use of `lmfit`_'s support for `mathemical constraints `_. +We will fit a circuit to the following impedance spectrum. + + +.. plot:: + + from pyimpspec import mpl + from pyimpspec import generate_mock_data + + data = generate_mock_data("CIRCUIT_5", noise=5e-2, seed=42)[0] + figure, axes = mpl.plot_nyquist(data, colors={"impedance": "black"}) + figure.tight_layout() + + +The circuit that we will use is given by the following CDC: ``R(RQ)(RQ)(RQ)``. + + +.. plot:: + + from pyimpspec import Circuit, parse_cdc + circuit: Circuit = parse_cdc("R(RQ)(RQ)(RQ)") + drawing = circuit.to_drawing() + drawing.draw() + + +If we simply fit this circuit to the impedance spectrum, then we will get a fit that is not great but manages to somewhat resemble the data. +Also, the fitted ``(RQ)`` units may not align nicely with what we can see in the Nyquist plot from left to right. +For example, ``R4`` and ``Q3`` correspond to the semicircle in the middle based on the values of their parameters when compared to the values of the parameters of the other ``(RQ)`` units. + + +.. plot:: + + from pyimpspec import ( + Circuit, + DataSet, + Element, + FitResult, + fit_circuit, + generate_mock_data, + mpl, + parse_cdc, + ) + from typing import ( + Dict, + List, + ) + + data: DataSet = generate_mock_data("CIRCUIT_5", noise=5e-2, seed=42)[0] + + circuit: Circuit = parse_cdc("R(RQ)(RQ)(RQ)") + + fit: FitResult = fit_circuit( + circuit, + data, + method="least_squares", + weight="boukamp", + ) + print(fit.to_parameters_dataframe().to_markdown(index=False)) + print(fit.to_statistics_dataframe().to_markdown(index=False)) + + figure, axes = mpl.plot_fit( + fit, + data=data, + colored_axes=True, + legend=False, + title="Without constraints", + ) + figure.tight_layout() + +.. code:: + + | Element | Parameter | Value | Std. err. (%) | Unit | Fixed | + |:----------|:------------|--------------:|----------------:|:-------|:--------| + | R_1 | R | 151.126 | 0.387059 | ohm | No | + | R_2 | R | 259.693 | 0.944897 | ohm | No | + | Q_1 | Y | 1.07068e-07 | 0.126065 | S*s^n | No | + | Q_1 | n | 0.950525 | 0.0555558 | | No | + | R_3 | R | 581.496 | 1.97511 | ohm | No | + | Q_2 | Y | 1.8216e-05 | 6.33997 | S*s^n | No | + | Q_2 | n | 0.85963 | 1.65839 | | No | + | R_4 | R | 726.502 | 0.695414 | ohm | No | + | Q_3 | Y | 4.18597e-07 | 0.325163 | S*s^n | No | + | Q_3 | n | 0.940476 | 0.115395 | | No | + + | Label | Value | + |:-------------------------------|:--------------------| + | Log pseudo chi-squared | -2.6523611472485586 | + | Log chi-squared | -6.941315663341954 | + | Log chi-squared (reduced) | -8.798648159773222 | + | Akaike info. criterion | -1651.9545159942043 | + | Bayesian info. criterion | -1627.8873235215617 | + | Degrees of freedom | 72 | + | Number of data points | 82 | + | Number of function evaluations | 209 | + | Method | least_squares | + | Weight | boukamp | + + +However, we can use constraints to enforce the following orders: ``R2 ≤ R4 ≤ R3`` and ``Q1 ≤ Q2 ≤ Q3``. To accomplish this, we need to provide the |fit_circuit| function with two additional arguments: ``constraint_expressions`` and ``constraint_variables``. + +The ``constraint_expressions`` argument maps the names/identifiers, which are provided to `lmfit`_ during the fitting process, of the constrained parameters to strings that contain the expressions that define the constraints. +The |generate_fit_identifiers| function can be used to obtain the names/identifiers of an element's parameters that will be available within the context where the constraint expressions are parsed and evaluated by `lmfit`_. + +The ``constraint_variables`` argument can be used to define any additional variables that would be used by the constraint expressions. +See `lmfit's documentation about Parameters `_ for information about valid keyword arguments. + +.. doctest:: + + >>> from pyimpspec import ( + ... Circuit, + ... DataSet, + ... Element, + ... FitIdentifiers, + ... FitResult, + ... fit_circuit, + ... generate_fit_identifiers, + ... generate_mock_data, + ... parse_cdc, + ... ) + >>> from typing import ( + ... Dict, + ... List, + ... ) + >>> + >>> data: DataSet = generate_mock_data("CIRCUIT_5", noise=5e-2, seed=42)[0] + >>> + >>> circuit: Circuit = parse_cdc("R(RQ)(RQ)(RQ)") + >>> identifiers: Dict[Element, FitIdentifiers] = generate_fit_identifiers(circuit) + >>> elements: List[Element] = circuit.get_elements() + >>> R1, R2, Q1, R3, Q2, R4, Q3 = elements + >>> + >>> fit: FitResult = fit_circuit( + ... circuit, + ... data, + ... method="least_squares", + ... weight="boukamp", + ... constraint_expressions={ + ... identifiers[R3].R: f"{identifiers[R2].R} + alpha", + ... identifiers[R4].R: f"{identifiers[R3].R} - beta", + ... identifiers[Q2].Y: f"{identifiers[Q1].Y} + gamma", + ... identifiers[Q3].Y: f"{identifiers[Q2].Y} + delta", + ... }, + ... constraint_variables=dict( + ... alpha=dict( + ... value=500, + ... min=0, + ... ), + ... beta=dict( + ... value=300, + ... min=0, + ... ), + ... gamma=dict( + ... value=1e-8, + ... min=0, + ... ), + ... delta=dict( + ... value=2e-7, + ... min=0, + ... ), + ... ), + ... ) + >>> + >>> R1, R2, Q1, R3, Q2, R4, Q3 = fit.circuit.get_elements() + >>> ( + ... (R2.get_value("R") < R4.get_value("R") < R3.get_value("R")) + ... and (Q1.get_value("Y") < Q2.get_value("Y") < Q3.get_value("Y")) + ... ) + True + + +Now we can obtain a better fit and also have all of the elements in an order that matches what we can observe in the Nyquist plot from left to right. + + +.. plot:: + + from pyimpspec import ( + Circuit, + DataSet, + Element, + FitIdentifiers, + FitResult, + fit_circuit, + generate_fit_identifiers, + generate_mock_data, + mpl, + parse_cdc, + ) + from typing import ( + Dict, + List, + ) + + data: DataSet = generate_mock_data("CIRCUIT_5", noise=5e-2, seed=42)[0] + + circuit: Circuit = parse_cdc("R(RQ)(RQ)(RQ)") + identifiers: Dict[Element, FitIdentifiers] = generate_fit_identifiers(circuit) + elements: List[Element] = circuit.get_elements() + R1, R2, Q1, R3, Q2, R4, Q3 = elements + + fit: FitResult = fit_circuit( + circuit, + data, + method="least_squares", + weight="boukamp", + constraint_expressions={ + identifiers[R3].R: f"{identifiers[R2].R} + alpha", + identifiers[R4].R: f"{identifiers[R3].R} - beta", + identifiers[Q2].Y: f"{identifiers[Q1].Y} + gamma", + identifiers[Q3].Y: f"{identifiers[Q2].Y} + delta", + }, + constraint_variables=dict( + alpha=dict( + value=500, + min=0, + ), + beta=dict( + value=300, + min=0, + ), + gamma=dict( + value=1e-8, + min=0, + ), + delta=dict( + value=2e-7, + min=0, + ), + ), + ) + print(fit.to_parameters_dataframe().to_markdown(index=False)) + print(fit.to_statistics_dataframe().to_markdown(index=False)) + + R1, R2, Q1, R3, Q2, R4, Q3 = fit.circuit.get_elements() + assert R2.get_value("R") < R4.get_value("R") < R3.get_value("R") + assert Q1.get_value("Y") < Q2.get_value("Y") < Q3.get_value("Y") + + refined_fit: FitResult = fit_circuit( + fit.circuit, + data, + method="least_squares", + weight="boukamp", + ) + + print(refined_fit.to_parameters_dataframe().to_markdown(index=False)) + print(refined_fit.to_statistics_dataframe().to_markdown(index=False)) + + figure, axes = mpl.plot_fit( + fit, + data=data, + colored_axes=True, + legend=False, + title="With constraints", + ) + figure.tight_layout() + + +.. code:: + + | Element | Parameter | Value | Std. err. (%) | Unit | Fixed | + |:----------|:------------|--------------:|----------------:|:-------|:--------| + | R_1 | R | 143.461 | 0.0917482 | ohm | No | + | R_2 | R | 229.453 | 0.288944 | ohm | No | + | Q_1 | Y | 1.78147e-07 | 0.00643312 | S*s^n | No | + | Q_1 | n | 0.912666 | 0.0138234 | | No | + | R_3 | R | 854.409 | nan | ohm | Yes | + | Q_2 | Y | 7.82691e-07 | nan | S*s^n | Yes | + | Q_2 | n | 0.856438 | 0.0285731 | | No | + | R_4 | R | 475.654 | nan | ohm | Yes | + | Q_3 | Y | 1.60541e-05 | nan | S*s^n | Yes | + | Q_3 | n | 0.952466 | 0.19436 | | No | + + | Label | Value | + |:-------------------------------|:--------------------| + | Log pseudo chi-squared | -4.2481804461232695 | + | Log chi-squared | -10.082848279294568 | + | Log chi-squared (reduced) | -11.940180775725837 | + | Akaike info. criterion | -2245.113501987264 | + | Bayesian info. criterion | -2221.0463095146215 | + | Degrees of freedom | 72 | + | Number of data points | 82 | + | Number of function evaluations | 569 | + | Method | least_squares | + | Weight | boukamp | + + +.. note:: + + Keep in mind that an estimate for the error of a constrained parameter's value will not be available. However, this can be remedied by using the obtained fitted circuit to perform yet another fit without any constraints. + + +.. code:: + + | Element | Parameter | Value | Std. err. (%) | Unit | Fixed | + |:----------|:------------|--------------:|----------------:|:-------|:--------| + | R_1 | R | 143.48 | 0.0892796 | ohm | No | + | R_2 | R | 229.69 | 0.283135 | ohm | No | + | Q_1 | Y | 1.78048e-07 | 0.00759739 | S*s^n | No | + | Q_1 | n | 0.912652 | 0.013558 | | No | + | R_3 | R | 853.869 | 0.0916479 | ohm | No | + | Q_2 | Y | 7.80619e-07 | 0.0126453 | S*s^n | No | + | Q_2 | n | 0.856809 | 0.0280145 | | No | + | R_4 | R | 475.967 | 0.230118 | ohm | No | + | Q_3 | Y | 1.60529e-05 | 0.697925 | S*s^n | No | + | Q_3 | n | 0.952274 | 0.190482 | | No | + + | Label | Value | + |:-------------------------------|:--------------------| + | Log pseudo chi-squared | -4.262336239351527 | + | Log chi-squared | -10.111392253122428 | + | Log chi-squared (reduced) | -11.968724749553697 | + | Akaike info. criterion | -2250.5029461349936 | + | Bayesian info. criterion | -2226.435753662351 | + | Degrees of freedom | 72 | + | Number of data points | 82 | + | Number of function evaluations | 79 | + | Method | least_squares | + | Weight | boukamp | + + .. raw:: latex \clearpage diff --git a/docs/source/guide_installing.rst b/docs/source/guide_installing.rst index 1972f45..2831056 100644 --- a/docs/source/guide_installing.rst +++ b/docs/source/guide_installing.rst @@ -60,9 +60,9 @@ For example, open a terminal and run the following command to confirm that pip ( Such a virtual environment needs to be activated before running a script that imports a package installed inside the virtual environment. The system-wide Python environment may also be `externally managed `_ in order to prevent the user from accidentally breaking that environment since the operating system depends upon the packages in that environment. - A third-party tool called `pipx `_ can automatically manage such virtual environments but it is primarily for installing programs that provide, e.g., a command-line interface (CLI) or a graphical user interface (GUI). + A third-party tool called `pipx `_ can automatically manage such virtual environments, but it is primarily for installing programs that provide, e.g., a command-line interface (CLI) or a graphical user interface (GUI). These programs can then be run without having to manually activate the virtual environment since pipx handles that. - The virtual environment would still need to be activated before running a script that imports DearEIS and makes use of DearEIS's application programming interface (API). + The virtual environment would still need to be activated before running a script that imports pyimpspec and makes use of pyimpspec's application programming interface (API). If using pipx, then run the following command to make sure that pipx is available. If pipx is not available, then follow the `instructions to install pipx `_. @@ -75,40 +75,42 @@ If pipx is not available, then follow the `instructions to install pipx Tuple[TimeConstants, Gammas, TimeConstants, Gammas]: """ - Get the time constants (in seconds) of peaks with magnitudes greater than the threshold. + Get the time constants (in seconds) and gammas (in ohms) of peaks with magnitudes greater than the threshold. The threshold and the magnitudes are all relative to the magnitude of the highest peak. Parameters @@ -284,7 +290,7 @@ def to_scores_dataframe( rows: Optional[List[str]] = None, ) -> "DataFrame": # noqa: F821 """ - Get the scores for the data set as a |DataFrame| object that can be used to generate, e.g., a Markdown table. + Get the scores for the data set as a `pandas.DataFrame`_ object that can be used to generate, e.g., a Markdown table. Parameters ---------- @@ -296,7 +302,7 @@ def to_scores_dataframe( Returns ------- - |DataFrame| + `pandas.DataFrame`_ """ from pandas import DataFrame @@ -357,6 +363,51 @@ def to_scores_dataframe( } ) + def analyze_peaks( + self, + num_peaks: int = 0, + peak_positions: Optional[Union[List[float], NDArray[float64]]] = None, + disallow_skew: bool = False, + ) -> Tuple[DRTPeaks, DRTPeaks]: + """ + Analyze the peaks present in a distribution of relaxation times using skew normal distributions. + + Parameters + ---------- + num_peaks: int, optional + If greater than zero, then analyze only that number of peaks (sorted from highest to lowest gamma values). + + peak_positions: Optional[Union[List[float], NDArray[float64]]], optional + Analyze only the peaks at the provided positions. + + disallow_skew: bool, optional + If true, then normal distributions are used instead of skew normal distributions. + + Returns + ------- + Tuple[|DRTPeaks|, |DRTPeaks|] + The first and second |DRTPeaks| instance corresponds to the DRT generated based on the real or imaginary part, respectively, of the impedance spectrum. + """ + peaks_real: DRTPeaks = _analyze_peaks( + self.time_constants, + self.real_gammas, + num_peaks=num_peaks, + peak_positions=peak_positions, + disallow_skew=disallow_skew, + suffix="real", + ) + peaks_imag: DRTPeaks = _analyze_peaks( + self.time_constants, + self.imaginary_gammas, + num_peaks=num_peaks, + peak_positions=peak_positions, + disallow_skew=disallow_skew, + suffix="imag.", + ) + object.__setattr__(self, "_peak_analysis", (peaks_real, peaks_imag)) + + return (peaks_real, peaks_imag) + def _compute_res_scores( res: NDArray[float64], @@ -1095,7 +1146,6 @@ def _calculate_model_impedance( ) -# TODO: Add support for admittance? def calculate_drt_bht( data: DataSet, rbf_type: str = "gaussian", @@ -1173,7 +1223,7 @@ def calculate_drt_bht( Returns ------- - BHTResult + |BHTResult| """ if not isinstance(rbf_type, str): raise TypeError(f"Expected a string instead of {rbf_type=}") diff --git a/src/pyimpspec/analysis/drt/lm.py b/src/pyimpspec/analysis/drt/lm.py new file mode 100644 index 0000000..053406c --- /dev/null +++ b/src/pyimpspec/analysis/drt/lm.py @@ -0,0 +1,714 @@ +# pyimpspec is licensed under the GPLv3 or later (https://www.gnu.org/licenses/gpl-3.0.html). +# Copyright 2024 pyimpspec developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# The licenses of pyimpspec's dependencies and/or sources of portions of code are included in +# the LICENSES folder. + +# This module implements the Loewner method +# 10.1016/j.laa/2007.03.008 +# 10.2139/ssrn.4217752 +# 10.3390/batteries9020132 +# 10.1016/j.jpowsour.2023.233575 + +# Based on source code from +# https://github.com/projectsEECandDRI/DRT-from-Loewner-framework +# (commit 19a40ea7d53d6d62ea5d7eeb63584818f2e36a7b) +# with some modifications. + +from dataclasses import dataclass +from numpy import ( + argwhere, + array, + complex128, + concatenate, + flip, + float64, + full, + int64, + log10 as log, + nan, + pi, + sqrt, + zeros, +) +from numpy.linalg import ( + solve, + svd, + matrix_rank, +) +from numpy.typing import NDArray +from pyimpspec.data import DataSet +from pyimpspec.analysis.kramers_kronig.single import ( + KramersKronigResult, + perform_kramers_kronig_test, +) +from pyimpspec.analysis.utility import ( + _calculate_pseudo_chisqr, + _calculate_residuals, + _get_default_num_procs, +) +from .result import DRTResult +from .peak_analysis import DRTPeaks +from pyimpspec.progress import Progress +from pyimpspec.typing import ( + ComplexImpedance, + ComplexImpedances, + Frequencies, + Gamma, + Gammas, + Indices, + TimeConstant, + TimeConstants, +) +from pyimpspec.typing.helpers import ( + Dict, + List, + NDArray, + Optional, + Tuple, + Union, + float64, + _is_integer, +) + + +_DEBUG: bool = bool(0) + +_MODEL_ORDER_METHODS: List[str] = [ + "matrix_rank", + "pseudo_chisqr", +] + + +@dataclass(frozen=True) +class LMResult(DRTResult): + """ + An object representing the results of calculating the distribution of relaxation times in a data set using the Loewner method. + + Parameters + ---------- + time_constants: |TimeConstants| + The time constants. + + gammas: Gammas + All gamma values. Positive and negative values correspond to resistive-capacitive and resistive-inductive peaks, respectively. + + frequencies: |Frequencies| + The frequencies of the impedance spectrum. + + impedances: |ComplexImpedances| + The impedance produced by the model. + + residuals: |ComplexResiduals| + The residuals of the real parts of the model and the data set. + + pseudo_chisqr: float + The pseudo chi-squared value, |pseudo chi-squared|, of the modeled impedance (eq. 14 in Boukamp, 1995). + + singular_values: NDArray[float64] + The singular values as a function of the model order :math:`k \\in [1..n]` where :math:`n` is the closest value divisible by 2 that is less than the total number of frequencies in the |DataSet| that was analyzed. + """ + + gammas: Gammas + singular_values: NDArray[float64] + + @property + def _resistive_capacitive_time_constants(self) -> NDArray[float64]: + indices: NDArray[int64] = argwhere(self.gammas >= 0.0).flatten() + + return self.time_constants[indices] + + @property + def _resistive_capacitive_gammas(self) -> NDArray[float64]: + indices: NDArray[int64] = argwhere(self.gammas >= 0.0).flatten() + + return self.gammas[indices] + + @property + def _resistive_inductive_time_constants(self) -> NDArray[float64]: + indices: NDArray[int64] = argwhere(self.gammas < 0.0).flatten() + + return self.time_constants[indices] + + @property + def _resistive_inductive_gammas(self) -> NDArray[float64]: + indices: NDArray[int64] = argwhere(self.gammas < 0.0).flatten() + + return abs(self.gammas[indices]) + + def get_label(self) -> str: + return f"LM (k={len(self.time_constants)})" + + def get_gammas(self) -> Tuple[Gammas, Gammas]: + """ + Get the gamma values for the resistive-capacitive and resistive-inductive peaks. + + Returns + ------- + Tuple[|Gammas|, |Gammas|] + """ + return ( + self._resistive_capacitive_gammas, + self._resistive_inductive_gammas, + ) + + def get_drt_data(self) -> Tuple[TimeConstants, Gammas, TimeConstants, Gammas]: + """ + Get the data necessary to plot this DRT result as a DRT plot: the time constants and corresponding gamma values for the resistive-capacitive and the resistive-inductive peaks. + + Returns + ------- + Tuple[|TimeConstants|, |Gammas|, |TimeConstants|, |Gammas|] + """ + return ( + self._resistive_capacitive_time_constants, + self._resistive_capacitive_gammas, + self._resistive_inductive_time_constants, + self._resistive_inductive_gammas, + ) + + def to_peaks_dataframe( + self, + threshold: float = 0.0, + columns: Optional[List[str]] = None, + ) -> "DataFrame": # noqa: F821 + from pandas import DataFrame + + if columns is None: + columns = [ + "tau, RC (s)", + "gamma, RC (ohm)", + "tau, RL (s)", + "gamma, RL (ohm)", + ] + elif not isinstance(columns, list): + raise TypeError(f"Expected a list of strings instead of {columns=}") + elif len(columns) != 4: + raise ValueError(f"Expected a list with 4 items instead of {len(columns)=}") + elif not all(map(lambda s: isinstance(s, str), columns)): + raise TypeError(f"Expected a list of strings instead of {columns=}") + elif len(set(columns)) != 4: + raise ValueError( + f"Expected a list of 4 unique strings instead of {columns=}" + ) + + def get_peak_indices(threshold, gammas): + if gammas.size > 0: + max_gamma = max(gammas) + else: + max_gamma = 1.0 + + return argwhere((gammas / max_gamma) > threshold).flatten() + + def pad( + tau: TimeConstants, + gamma: Gammas, + width: int, + ) -> Tuple[TimeConstants, Gammas]: + tmp_tau: TimeConstants = full(width, nan, dtype=TimeConstant) + tmp_tau[: tau.size] = tau + tmp_gamma: Gammas = full(width, nan, dtype=Gamma) + tmp_gamma[: gamma.size] = gamma + return ( + tmp_tau, + tmp_gamma, + ) + + indices_RC: Indices = get_peak_indices( + threshold, + self._resistive_capacitive_gammas, + ) + indices_RL: Indices = get_peak_indices( + threshold, + self._resistive_inductive_gammas, + ) + + width: int = max( + ( + indices_RC.size, + indices_RL.size, + ) + ) + + + tau_RC: TimeConstants + gamma_RC: Gammas + tau_RC, gamma_RC = pad( + tau=self._resistive_capacitive_time_constants[indices_RC], + gamma=self._resistive_capacitive_gammas[indices_RC], + width=width, + ) + + tau_RL: TimeConstants + gamma_RL: Gammas + tau_RL, gamma_RL = pad( + tau=self._resistive_inductive_time_constants[indices_RL], + gamma=self._resistive_inductive_gammas[indices_RL], + width=width, + ) + + return DataFrame.from_dict( + { + columns[0]: tau_RC, + columns[1]: gamma_RC, + columns[2]: tau_RL, + columns[3]: gamma_RL, + } + ) + + def to_statistics_dataframe( + self, + ) -> "DataFrame": # noqa: F821 + from pandas import DataFrame + + statistics: Dict[str, Union[int, float, str]] = { + "Log pseudo chi-squared": log(self.pseudo_chisqr), + "Model order (k)": len(self.time_constants), + } + + return DataFrame.from_dict( + { + "Label": list(statistics.keys()), + "Value": list(statistics.values()), + } + ) + + def get_peaks( + self, + threshold: float = 0.0, + ) -> Tuple[TimeConstants, Gammas, TimeConstants, Gammas]: + """ + Get the time constants (in seconds) of peaks with magnitudes greater than the threshold. + The threshold and the magnitudes are all relative to the magnitude of the highest peak. + + Parameters + ---------- + threshold: float, optional + The minimum peak height threshold (relative to the height of the tallest peak) for a peak to be included. + + Returns + ------- + Tuple[|TimeConstants|, |Gammas|, |TimeConstants|, |Gammas|] + """ + def get_peak_indices(threshold: float, gamma: Gammas) -> Indices: + if len(gamma) == 0: + return array([], dtype=int64) + + max_g: float = max(gamma) + if max_g == 0.0: + return array([], dtype=int64) + + indices: Indices = array(list(range(len(gamma))), dtype=int64) + + return array( + list( + filter( + lambda i: gamma[i] / max_g > threshold and gamma[i] > 0.0, indices + ) + ), + dtype=int64, + ) + + indices_RC: Indices = get_peak_indices(threshold, self._resistive_capacitive_gammas) + indices_RL: Indices = get_peak_indices(threshold, self._resistive_inductive_gammas) + + return ( + self._resistive_capacitive_time_constants[indices_RC], + self._resistive_capacitive_gammas[indices_RC], + self._resistive_inductive_time_constants[indices_RL], + self._resistive_inductive_gammas[indices_RL], + ) + + def analyze_peaks( + self, + num_peaks: int = 0, + peak_positions: Optional[Union[List[float], NDArray[float64]]] = None, + disallow_skew: bool = False, + ) -> Tuple[DRTPeaks, DRTPeaks]: + """ + Analyze the peaks present in a distribution of relaxation times using skew normal distributions. + + **Peak analysis of the DRT obtained using the Loewner method is not supported.** + + Parameters + ---------- + num_peaks: int, optional + If greater than zero, then analyze only that number of peaks (sorted from highest to lowest gamma values). + + peak_positions: Optional[Union[List[float], NDArray[float64]]], optional + Analyze only the peaks at the provided positions. + + disallow_skew: bool, optional + If true, then normal distributions are used instead of skew normal distributions. + + Returns + ------- + Tuple[|DRTPeaks|, |DRTPeaks|] + Neither |DRTPeaks| instance actually contains any results corresponding to peaks. + """ + return ( + DRTPeaks( + self._resistive_capacitive_time_constants, + peaks=[], + suffix="RC", + ), + DRTPeaks( + self._resistive_inductive_time_constants, + peaks=[], + suffix="RL", + ), + ) + + +def _interleave_complex_conjugate(values: NDArray[complex128]) -> NDArray[complex128]: + dst: NDArray[complex128] = zeros(len(values) * 2, dtype=complex128) + dst[0::2] += values.T + dst[1::2] += values.conj().T + + return dst + + +def _generate_loewner_matrices( + s: NDArray[complex128], + Z: ComplexImpedances, +) -> Tuple[ + NDArray[float64], + NDArray[float64], + NDArray[float64], + NDArray[float64], +]: + from scipy.linalg import block_diag + + # Establish the number of points to use + n: int = len(s) + while n % 2 != 0: + n -= 1 + + # Generate the vectors of values for the two disjointed sets + la: NDArray[complex128] = _interleave_complex_conjugate(s[0::2][:n // 2]) + W: NDArray[complex128] = _interleave_complex_conjugate(Z[0::2][:n // 2]) + mu: NDArray[complex128] = _interleave_complex_conjugate(s[1::2][:n // 2]) + V: NDArray[complex128] = _interleave_complex_conjugate(Z[1::2][:n // 2]) + + # Generate the Loewner and shifted Loewner matrices + L: NDArray[complex128] = zeros((n, n), dtype=complex128) + Ls: NDArray[complex128] = zeros((n, n), dtype=complex128) + + i: int + for i in range(0, len(mu)): + j: int + for j in range(0, len(la)): + L[i][j] = (V[i] - W[j]) / (mu[i] - la[j]) + Ls[i][j] = (mu[i] * V[i] - la[j] * W[j]) / (mu[i] - la[j]) + + # Transform the complex-valued matrices into real-valued matrices + block: NDArray[complex128] = array( + [ + [1, 1j], + [1, -1j] + ] + ) / sqrt(2) + P: NDArray[complex128] = block_diag(*[block for _ in range(0, n//2)]) + Pstar: NDArray[complex128] = P.conj().T + + L = Pstar @ L @ P + Ls = Pstar @ Ls @ P + W = (W.T @ P).T + V = Pstar @ V + + return ( + L.real, + Ls.real, + V.real, + W.real, + ) + +def _generate_reduced_order_model( + L: NDArray[float64], + Ls: NDArray[float64], + V: NDArray[float64], + W: NDArray[float64], + k: int, +) -> Tuple[ + NDArray[float64], + NDArray[float64], + NDArray[float64], + NDArray[float64], + NDArray[float64], +]: + singular_values: NDArray[float64] + _, singular_values, _ = svd(concatenate((L, Ls), 1)) + + Y_L: NDArray[float64] + X_L: NDArray[float64] + Y_L, _, X_L = svd(L) + + Yk: NDArray[float64] = Y_L[:, :k] + Xk: NDArray[float64] = X_L[:, :k] + + Ek = -Yk.T @ L @ Xk + Ak = -Yk.T @ Ls @ Xk + Bk = Yk.T @ V + Ck = W.T @ Xk + + return ( + Ek, + Ak, + Bk, + Ck, + singular_values, + ) + + +def _calculate_model_impedance( + s: NDArray[complex128], + Ek: NDArray[float64], + Ak: NDArray[float64], + Bk: NDArray[float64], + Ck: NDArray[float64], +) -> ComplexImpedances: + Z_fit: ComplexImpedances = zeros(s.shape, dtype=ComplexImpedance) + + i: int + for i in range(0, len(s)): + A: NDArray[complex128] = s[i] * Ek - Ak + B: NDArray[complex128] = solve(A, Bk) + Z_fit[i] = Ck @ B + + return Z_fit + + + +def _extract_peaks( + Ek: NDArray[float64], + Ak: NDArray[float64], + Bk: NDArray[float64], + Ck: NDArray[float64], +) -> Tuple[TimeConstants, Gammas]: + from scipy.linalg import eig + + eigenvalues: NDArray[complex128] + eigenvectors: NDArray[complex128] + eigenvalues, eigenvectors = eig(Ak, Ek) + + Bt: NDArray[complex128] = solve(eigenvectors, solve(Ek, Bk)) + Ct: NDArray[complex128] = Ck @ eigenvectors + residues: NDArray[complex128] = Bt * Ct.T + + time_constants: TimeConstants = abs(-1 / eigenvalues) + gammas: Gammas = (-residues / eigenvalues).real + + return ( + time_constants, + gammas, + ) + + +def calculate_drt_lm( + data: DataSet, + model_order: int = 0, + model_order_method: str = "matrix_rank", + num_procs: int = -1, + **kwargs, +) -> LMResult: + """ + Calculates the distribution of relaxation times (DRT) using the Loewner method (LM). + + References: + + - `Mayo, A.J. and Antoulas, A.C., 2007, Linear Algebra Its Appl. 425, 634–662 `_, + - `Sorrentino, A., Patel, B., Gosea, I.V., Antoulas, A.C., and Vidaković-Koch, T., 2022, SSRN `_ + - `Rüther, T., Gosea, I.V., Jahn, L., Antoulas, A.C., and Danzer, M.A., 2023, Batteries, 9, 132 `_ + - `Sorrentino, A., Patel, B., Gosea, I.V., Antoulas, A.C., and Vidaković-Koch, T., 2023, J. Power Sources, 585, 223575 `_ + + Parameters + ---------- + data: DataSet + The data set to use in the calculations. + + model_order: int, optional + The order of the model (k). + + model_order_method: str, optional + How to automatically pick the order of the model if the model order is not specified: + + - "matrix_rank" + - "pseudo_chisqr" + + num_procs: int, optional + The maximum number of parallel processes to use. + A value less than 1 results in an attempt to figure out a suitable value based on, e.g., the number of cores detected. + Additionally, a negative value can be used to reduce the number of processes by that much (e.g., to leave one core for a GUI thread). + + Returns + ------- + LMResult + """ + f: Frequencies = data.get_frequencies() + n: int = len(f) + while n % 2 != 0: + n -= 1 + + if n <= 2: + raise ValueError(f"Expected more frequencies (i.e., {len(f)=} > 2)") + + if not _is_integer(model_order): + raise TypeError(f"Expected an integer instead of {model_order=}") + elif not (model_order <= n): + raise ValueError(f"Expected {model_order=} <= {n=}") + + if not isinstance(model_order_method, str): + raise TypeError(f"Expected a string instead of {model_order_method=}") + elif model_order_method not in _MODEL_ORDER_METHODS: + raise ValueError(f"Unsupported {model_order_method=} encountered instead of one of the following:\n- {'\n- '.join(_MODEL_ORDER_METHODS)}") + + if not _is_integer(num_procs): + raise TypeError(f"Expected an integer instead of {num_procs=}") + elif num_procs < 1: + num_procs = max((_get_default_num_procs() - abs(num_procs), 1)) + + prog: Progress + with Progress("Preparing matrices", total=3) as prog: + s: NDArray[complex128] = flip(1j * 2 * pi * f) + Z_exp: ComplexImpedances = flip(data.get_impedances()) + + L: NDArray[float64] + Ls: NDArray[float64] + V: NDArray[float64] + W: NDArray[float64] + L, Ls, V, W = _generate_loewner_matrices(s, Z_exp) + prog.increment() + prog.set_message("Generating reduced model") + + Z_exp = flip(Z_exp) + + reduced_model_has_been_generated: bool = False + Ek: NDArray[float64] + Ak: NDArray[float64] + Bk: NDArray[float64] + Ck: NDArray[float64] + Z_fit: ComplexImpedances + pseudo_chisqr: float + + # If the model order has not been specified, then pick one + k: int = model_order + if k < 1 and model_order_method == "matrix_rank": + k = matrix_rank(concatenate((L, Ls), 1)) + elif k < 1 and model_order_method == "pseudo_chisqr": + test: KramersKronigResult = perform_kramers_kronig_test( + data, + test=kwargs.get("test", "complex"), + num_RC=kwargs.get("num_RC", 0), + add_capacitance=kwargs.get("add_capacitance", True), + add_inductance=kwargs.get("add_inductance", True), + admittance=kwargs.get("admittance", None), + min_log_F_ext=kwargs.get("min_log_F_ext", -1.0), + max_log_F_ext=kwargs.get("max_log_F_ext", 1.0), + log_F_ext=kwargs.get("log_F_ext", 0.0), + num_F_ext_evaluations=kwargs.get("num_F_ext_evaluations", 20), + rapid_F_ext_evaluations=kwargs.get("rapid_F_ext_evaluations", True), + cnls_method=kwargs.get("cnls_method", "leastsq"), + max_nfev=kwargs.get("max_nfev", 0), + timeout=kwargs.get("timeout", 60), + num_procs=kwargs.get("num_procs", -1), + ) + + k_values: List[int] = [] + pseudo_chisqr_values: List[float] = [] + + for k in range(1, n+1): + k_values.append(k) + Ek, Ak, Bk, Ck, singular_values = _generate_reduced_order_model( + L, + Ls, + V, + W, + k=k, + ) + reduced_model_has_been_generated = True + + Z_fit = _calculate_model_impedance( + s, + Ek, + Ak, + Bk, + Ck, + ) + + Z_fit = flip(Z_fit) + pseudo_chisqr = _calculate_pseudo_chisqr(Z_exp=Z_exp, Z_fit=Z_fit) + pseudo_chisqr_values.append(pseudo_chisqr) + if pseudo_chisqr <= test.pseudo_chisqr: + break + + if _DEBUG and reduced_model_has_been_generated: + import matplotlib.pyplot as plt + from pandas import Series + + fig, ax = plt.subplots() + + ax.set_xlim(0, n + 1) + ax.axhline(log(test.pseudo_chisqr), color="black", linestyle=":") + ax.scatter(k_values, log(pseudo_chisqr_values), color="blue", marker="o") + + ax2 = ax.twinx() + std = Series(log(pseudo_chisqr_values)).rolling(3, center=True, min_periods=1).std() + ax2.scatter(k_values, std, color="red", marker="+") + + plt.show() + elif k < 1: + raise NotImplementedError(f"Unsupported {model_order_method=}") + + if not reduced_model_has_been_generated: + Ek, Ak, Bk, Ck, singular_values = _generate_reduced_order_model( + L, + Ls, + V, + W, + k=k, + ) + + Z_fit = _calculate_model_impedance( + s, + Ek, + Ak, + Bk, + Ck, + ) + + Z_fit = flip(Z_fit) + pseudo_chisqr = _calculate_pseudo_chisqr(Z_exp=Z_exp, Z_fit=Z_fit) + + reduced_model_has_been_generated = True + + assert reduced_model_has_been_generated + + prog.increment() + prog.set_message("Calculating time constants") + + time_constants: TimeConstants + gammas: Gammas + time_constants, gammas = _extract_peaks(Ek, Ak, Bk, Ck) + + return LMResult( + time_constants=time_constants, + gammas=gammas, + frequencies=f, + impedances=Z_fit, + residuals=_calculate_residuals(Z_exp, Z_fit), + pseudo_chisqr=pseudo_chisqr, + singular_values=singular_values, + ) diff --git a/src/pyimpspec/analysis/drt/mrq_fit.py b/src/pyimpspec/analysis/drt/mrq_fit.py index 7315c25..e3aa7e5 100644 --- a/src/pyimpspec/analysis/drt/mrq_fit.py +++ b/src/pyimpspec/analysis/drt/mrq_fit.py @@ -109,7 +109,7 @@ class MRQFitResult(DRTResult): pseudo_chisqr: float The pseudo chi-squared value, |pseudo chi-squared|, of the modeled impedance (eq. 14 in Boukamp, 1995). - circuit: Circuit + circuit: |Circuit| The fitted circuit. """ @@ -522,7 +522,7 @@ def calculate_drt_mrq_fit( Returns ------- - MRQFitResult + |MRQFitResult| """ if not isinstance(circuit, Circuit): raise TypeError(f"Expected a Circuit instead of {circuit=}") diff --git a/src/pyimpspec/analysis/drt/peak_analysis.py b/src/pyimpspec/analysis/drt/peak_analysis.py new file mode 100644 index 0000000..eab8909 --- /dev/null +++ b/src/pyimpspec/analysis/drt/peak_analysis.py @@ -0,0 +1,520 @@ +# pyimpspec is licensed under the GPLv3 or later (https://www.gnu.org/licenses/gpl-3.0.html). +# Copyright 2024 pyimpspec developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# The licenses of pyimpspec's dependencies and/or sources of portions of code are included in +# the LICENSES folder. + +# This module implements analysis of peaks in DRT spectra by fitting skew +# normal distributions.See 10.3390/batteries5030053 for more information. + +# TODO: Implement tests for this feature + +from pyimpspec.progress import Progress +from dataclasses import dataclass +from scipy.signal import find_peaks +from numpy import ( + argmin, + argwhere, + copy, + exp, + float64, + isclose, + log10 as log, + logspace, + sign, + zeros, +) +from numpy.typing import NDArray +from pyimpspec.typing import ( + Gamma, + Gammas, + Indices, + TimeConstants, +) +from pyimpspec.typing.helpers import ( + Dict, + List, + Optional, + Tuple, + Union, + _is_boolean, + _is_floating_array, + _is_floating_list, + _is_integer, + _is_integer_array, + _is_integer_list, +) + + +def _skew_normal(x, h, p, a, s): + dp = x - p + den = 2 * s**2 + num = (dp * (1 + a * sign(dp)))**2 + + return h * exp(-num/den) + + +def _function(x, parameters: "Parameters") -> NDArray[float64]: + params = parameters.valuesdict() + assert len(params) % 2 == 0, (len(params), params.keys()) + + gammas: NDArray[float64] = zeros(len(x), dtype=float64) + + i = 0 + while True: + try: + h = params[f"h_{i}"] + except KeyError: + break + p = params[f"p_{i}"] + alpha = params[f"alpha_{i}"] + sigma = params[f"sigma_{i}"] + + gammas += _skew_normal(x, h, p, alpha, sigma) + i += 1 + + return gammas + + +def _residual( + parameters: "Parameters", + x: NDArray[float64], + y: NDArray[float64], +): + return _function(x, parameters) - y + + +@dataclass(frozen=True) +class DRTPeak: + """ + An object that represents a peak in a distribution of relaxation times. + The peak is modeled using a (skew) normal distribution. + + Parameters + ---------- + position: float64 + The relative position of the (skew) normal distribution. + + height: float64 + The relative height of the (skew) normal distribution. + + alpha: float64 + The skew of the (skew) normal distribution. + + sigma: float64 + The standard deviation of the (skew) normal distribution. + + x_offset: float64 + The offset used to translate time constants to the relative x-axis. + + x_scale: float64 + The scaling factor used to translate time constants to the relative x-axis. + + y_offset: float64 + The offset used to translate the relative y-axis values to gamma (ohm). + + y_scale: float64 + The scaling factor used to translate the relative y-axis values to gamma (ohm). + """ + position: float64 # Relative x-axis + height: float64 # Relative y-axis + alpha: float64 # Skew + sigma: float64 # Standard deviation + x_offset: float64 # For translating from time constant to relative x-axis + x_scale: float64 # For translating from time constant to relative x-axis + y_offset: float64 # For translating relative y-axis to gamma (ohm) + y_scale: float64 # For translating relative y-axis to gamma (ohm) + + def get_gammas( + self, + time_constants: TimeConstants, + ) -> Gammas: + """ + Calculate the gamma values (ohm) corresponding to the given time constants. + + Parameters + ---------- + time_constants: |TimeConstants| + + Returns + ------- + |Gammas| + """ + x: NDArray[float64] = (log(time_constants) - self.x_offset) / self.x_scale + + return (_skew_normal( + x=x, + h=self.height, + p=self.position, + a=self.alpha, + s=self.sigma, + ) * self.y_scale) + self.y_offset + + +@dataclass(frozen=True) +class DRTPeaks: + """ + An object that represents a collection of peaks in a distribution of relaxation times. + The peaks are modeled using a (skew) normal distribution. + + Parameters + ---------- + time_constants: |TimeConstants| + The time constants of the DRT. + + peaks: List[|DRTPeak| + A list of the (skew) normal distributions fitted to the DRT peaks. + + suffix: str + The suffix used in the labels when plotting the results. + """ + time_constants: TimeConstants + peaks: List[DRTPeak] + suffix: str + + def __iter__(self): + return iter(self.peaks) + + def get_num_peaks(self) -> int: + """ + Get the number of peaks. + + Returns + ------- + int + """ + return len(self.peaks) + + def get_time_constants( + self, + num_per_decade: int = 100, + ) -> TimeConstants: + """ + Get either the original time constants (``num_per_decade < 1``) or an interpolated set of time constants. + + Parameters + ---------- + num_per_decade: int, optional + The number of points per decade to use when generating the time constants. + + Returns + ------- + |TimeConstants| + """ + if not _is_integer(num_per_decade): + raise TypeError(f"Expected an integer instead of {num_per_decade=}") + + if num_per_decade < 1: + return self.time_constants + + log_tau_min = log(min(self.time_constants)) + log_tau_max = log(max(self.time_constants)) + + return logspace( + log_tau_min, + log_tau_max, + num=int(round(log_tau_max - log_tau_min) * num_per_decade) + 1, + ) + + def get_gammas( + self, + peak_indices: Optional[List[int]] = None, + num_per_decade: int = 100, + ) -> Gammas: + """ + Calculate the gamma values (ohm) for one or more peaks. + + Parameters + ---------- + peak_indices: Optional[List[int]], optional + If indices are specified, then only those peaks will be used. + + num_per_decade: int, optional + The number of points per decade to use when generating the gamma values. + + Returns + ------- + |Gammas| + """ + num_peaks: int = self.get_num_peaks() + if peak_indices is None: + pass + elif not (_is_integer_list(peak_indices) or _is_integer_array(peak_indices)): + raise TypeError(f"Expected an array of integers instead of {peak_indices=}") + elif not all(map(lambda i: 0 <= i < num_peaks, peak_indices)): + raise ValueError(f"Expected 0 <= {peak_indices=} < {num_peaks}") + + time_constants: TimeConstants = self.get_time_constants( + num_per_decade=num_per_decade, + ) + + gammas: Gammas = zeros(time_constants.shape, dtype=Gamma) + + i: int + peak: DRTPeak + for i, peak in enumerate(self.peaks): + if peak_indices and i not in peak_indices: + continue + + gammas += peak.get_gammas(time_constants=time_constants) + + return gammas + + def get_peak_area( + self, + index: int, + ) -> float64: + """ + Calculate the area of a peak. + + Parameters + ---------- + index: int + The index (zero-based) of the peak for which to calculate the area. + + Returns + ------- + float64 + """ + from scipy.integrate import quad + + num_peaks: int = self.get_num_peaks() + if not _is_integer(index): + raise TypeError(f"Expected an integer instead of {index=}") + elif not (0 <= index < num_peaks): + raise ValueError(f"Expected 0 <= {index=} < {num_peaks}") + + peak: DRTPeak = self.peaks[index] + + return quad( + func=lambda x: peak.get_gammas(10**x), + a=log(min(self.time_constants)), + b=log(max(self.time_constants)), + )[0] / log(exp(1)) + + def to_peaks_dataframe( + self, + peak_indices: Optional[List[int]] = None, + ) -> "DataFrame": # noqa: F821 + """ + Generate a `pandas.DataFrame`_ object containing information about the peaks. + + Parameters + ---------- + peak_indices: Optional[List[int]], optional + If indices are specified, then only those peaks will be used. + + Returns + ------- + `pandas.DataFrame`_ + """ + from pandas import DataFrame + + num_peaks: int = self.get_num_peaks() + if peak_indices is None: + pass + elif not (_is_integer_list(peak_indices) or _is_integer_array(peak_indices)): + raise TypeError(f"Expected an array of integers instead of {peak_indices=}") + elif not all(map(lambda i: 0 <= i < num_peaks, peak_indices)): + raise ValueError(f"Expected 0 <= {peak_indices=} < {num_peaks}") + + positions: List[float64] = [] + heights: List[float64] = [] + areas: List[float64] = [] + + i: int + peak: DRTPeak + for i, peak in enumerate(self.peaks): + if peak_indices and i not in peak_indices: + continue + + positions.append(10**((peak.position * peak.x_scale) + peak.x_offset)) + heights.append((peak.height * peak.y_scale) + peak.y_offset) + areas.append(self.get_peak_area(i)) + + suffix: str = "" + if self.suffix != "": + suffix = f", {self.suffix}" + + return DataFrame.from_dict({ + "tau" + suffix + " (s)": positions, + "gamma" + suffix + " (ohm)": heights, + "R_peak" + suffix + " (ohm)": areas, + }) + + +def _generate_parameters( + peaks: List[Tuple[float64, float64]], + disallow_skew: bool, +) -> "Parameters": + from lmfit.parameter import Parameters + + parameters: Parameters = Parameters() + + i: int + x: float64 + y: float64 + for i, (x, y) in enumerate(peaks): + assert 0.0 <= x <= 1.0, x + parameters.add( + name=f"h_{i}", + value=y, + min=0.0, + max=2.0, + ) + + assert 0.0 <= x <= 1.0 + kw = dict( + name=f"p_{i}", + value=x, + min=x * 0.99, + max=x * 1.01, + ) + + if isclose(kw["min"], 0.0): + kw["min"] = -1e-10 + + if isclose(kw["max"], 1.0): + kw["max"] = 1.0 + 1e-10 + + parameters.add(**kw) + + parameters.add( + name=f"alpha_{i}", + value=0.0, + min=-1.0, + max=1.0, + vary=not disallow_skew, + ) + + parameters.add( + name=f"sigma_{i}", + value=0.05, + min=1e-10, + max=1e2, + ) + + return parameters + + +def _analyze_peaks( + time_constants: TimeConstants, + gammas: Gammas, + num_peaks: int, + peak_positions: Optional[Union[List[float], NDArray[float64]]], + disallow_skew: bool = False, + suffix: str = "", +) -> DRTPeaks: + from lmfit.minimizer import minimize, MinimizerResult + from lmfit.parameter import Parameters + + if not _is_floating_array(time_constants): + raise TypeError(f"Expected an array of floats instead of {time_constants=}") + elif not _is_floating_array(gammas): + raise TypeError(f"Expected an array of floats instead of {gammas=}") + elif not _is_integer(num_peaks): + raise TypeError(f"Expected an integer instead of {num_peaks=}") + elif not ( + peak_positions is None + or _is_floating_array(peak_positions) + or _is_floating_list(peak_positions) + ): + raise TypeError(f"Expected either None or a list/array of floats instead of {peak_positions=}") + elif not _is_boolean(disallow_skew): + raise TypeError(f"Expected a boolean instead of {disallow_skew=}") + elif not isinstance(suffix, str): + raise TypeError(f"Expected a string instead of {suffix=}") + + prog: Progress + with Progress( + "Scaling data", + total=4, + ) as prog: + x: NDArray[float64] = log(time_constants) + x_offset: float64 = min(x) + x -= x_offset + x_scale: float64 = max(x) + x /= x_scale + assert isclose(min(x), 0.0) + assert isclose(max(x), 1.0) + prog.increment() + + y: NDArray[float64] = copy(gammas) + indices: Indices = argwhere(y < 0.0).flatten() + y[indices] = 0.0 + y_offset: float64 = min(y) + y = y - y_offset + y_scale: float64 = max(y) + y /= y_scale + assert isclose(min(y), 0.0), min(y) + assert isclose(max(y), 1.0), max(y) + prog.increment() + + prog.set_message("Detecting peaks") + peaks: List[Tuple[float64, float64]] = [] + i: int + + if peak_positions is not None and len(peak_positions) > 0: + pos: float + for pos in peak_positions: + i = argmin(abs(time_constants - pos)) + peaks.append((x[i], y[i])) + else: + y_ext: NDArray[float64] = zeros(len(y) + 2, dtype=float64) + y_ext[1:-1] += y + + for i in find_peaks(y_ext)[0]: + i -= 1 + peaks.append((x[i], y[i])) + + if num_peaks > 0: + peaks.sort(key=lambda t: t[1], reverse=True) + peaks = peaks[:num_peaks] + + peaks.sort(key=lambda t: t[0]) + + prog.increment() + + prog.set_message(f"Fitting {len(peaks)} peaks") + parameters: Parameters = _generate_parameters(peaks, disallow_skew=disallow_skew) + + fit: MinimizerResult = minimize( + _residual, + parameters, + method="leastsq", + args=(x, y), + ) + params: Dict[str, float64] = fit.params.valuesdict() + + drt_peaks: List[DRTPeak] = [] + + for i in range(0, len(peaks)): + drt_peaks.append(DRTPeak( + position=params[f"p_{i}"], + height=params[f"h_{i}"], + alpha=params[f"alpha_{i}"], + sigma=params[f"sigma_{i}"], + x_offset=x_offset, + x_scale=x_scale, + y_offset=y_offset, + y_scale=y_scale, + )) + + return DRTPeaks( + time_constants=time_constants, + peaks=drt_peaks, + suffix=suffix, + ) diff --git a/src/pyimpspec/analysis/drt/result.py b/src/pyimpspec/analysis/drt/result.py index b3898d0..a6156d3 100644 --- a/src/pyimpspec/analysis/drt/result.py +++ b/src/pyimpspec/analysis/drt/result.py @@ -26,6 +26,12 @@ angle, array, int64, + zeros, +) +from .peak_analysis import ( + DRTPeaks, + DRTPeak, + _analyze_peaks, ) from pyimpspec.typing import ( ComplexImpedances, @@ -40,9 +46,14 @@ ) from pyimpspec.typing.helpers import ( List, + NDArray, Optional, Tuple, + Union, + float64, + _is_boolean, _is_floating, + _is_floating_array, ) @@ -119,7 +130,7 @@ def to_peaks_dataframe( columns: Optional[List[str]] = None, ) -> "DataFrame": # noqa: F821 """ - Get the peaks as a |DataFrame| object that can be used to generate, e.g., a Markdown table. + Get the peaks as a `pandas.DataFrame`_ object that can be used to generate, e.g., a Markdown table. Parameters ---------- @@ -131,11 +142,11 @@ def to_peaks_dataframe( Returns ------- - |DataFrame| + `pandas.DataFrame`_ """ pass - def _get_peak_indices(self, threshold: float, gamma: Gammas) -> Indices: + def _get_peak_indices(self, threshold: float, gammas: Gammas) -> Indices: from scipy.signal import find_peaks if not _is_floating(threshold): @@ -143,18 +154,23 @@ def _get_peak_indices(self, threshold: float, gamma: Gammas) -> Indices: elif not (0.0 <= threshold <= 1.0): raise ValueError(f"Expected a value in the range [0.0, 1.0] instead of {threshold=}") - indices: Indices = find_peaks(gamma)[0] + padded_gammas: Gammas = zeros(gammas.size + 2, dtype=gammas.dtype) + padded_gammas[1:-1] = gammas + + indices: Indices = find_peaks(padded_gammas)[0] if not indices.any(): return array([], dtype=int64) - max_g: float = max(gamma) + max_g: float = max(gammas) if max_g == 0.0: return array([], dtype=int64) + indices -= 1 # Because of the padding + return array( list( filter( - lambda _: gamma[_] / max_g > threshold and gamma[_] > 0.0, indices + lambda _: gammas[_] / max_g > threshold and gammas[_] > 0.0, indices ) ), dtype=int64, @@ -206,10 +222,56 @@ def to_statistics_dataframe( self, ) -> "DataFrame": # noqa: F821 """ - Get the statistics related to the DRT as a |DataFrame| object. + Get the statistics related to the DRT as a `pandas.DataFrame`_ object. Returns ------- - |DataFrame| + `pandas.DataFrame`_ """ pass + + def analyze_peaks( + self, + num_peaks: int = 0, + peak_positions: Optional[Union[List[float], NDArray[float64]]] = None, + disallow_skew: bool = False, + ) -> DRTPeaks: + """ + Analyze the peaks present in a distribution of relaxation times using skew normal distributions. + + Parameters + ---------- + num_peaks: int, optional + If greater than zero, then analyze only that number of peaks (sorted from highest to lowest gamma values). + + peak_positions: Optional[Union[List[float], NDArray[float64]]], optional + Analyze only the peaks at the provided positions. + + disallow_skew: bool, optional + If true, then normal distributions are used instead of skew normal distributions. + + Returns + ------- + |DRTPeaks| + """ + args = self.get_drt_data() + if not isinstance(args, tuple): + raise TypeError(f"Expected a tuple instead of {args=}") + elif len(args) != 2: + raise ValueError(f"Expected a tuple with two values instead of {args=}") + elif not _is_floating_array(args[0]): + raise TypeError(f"Expected an array of floats instead of {args[0]=}") + elif not _is_floating_array(args[1]): + raise TypeError(f"Expected an array of floats instead of {args[1]=}") + elif not _is_boolean(disallow_skew): + raise TypeError(f"Expected a boolean instead of {disallow_skew=}") + + peaks: DRTPeaks = _analyze_peaks( + *args, + num_peaks=num_peaks, + peak_positions=peak_positions, + disallow_skew=disallow_skew, + ) + object.__setattr__(self, "_peak_analysis", peaks) + + return peaks diff --git a/src/pyimpspec/analysis/drt/tr_nnls.py b/src/pyimpspec/analysis/drt/tr_nnls.py index acef00c..186d6fa 100644 --- a/src/pyimpspec/analysis/drt/tr_nnls.py +++ b/src/pyimpspec/analysis/drt/tr_nnls.py @@ -134,15 +134,14 @@ def to_peaks_dataframe( f"Expected a list of 2 unique strings instead of {columns=}" ) - indices: Indices = self._get_peak_indices( - threshold, - self.gammas, # type: ignore - ) + time_constants: TimeConstants + gammas: Gammas + time_constants, gammas = self.get_peaks(threshold=threshold) return DataFrame.from_dict( { - columns[0]: self.time_constants[indices], # type: ignore - columns[1]: self.gammas[indices], # type: ignore + columns[0]: time_constants, # type: ignore + columns[1]: gammas, # type: ignore } ) @@ -445,7 +444,7 @@ def calculate_drt_tr_nnls( Returns ------- - TRNNLSResult + |TRNNLSResult| """ if not isinstance(mode, str): raise TypeError(f"Expected a string instead of {mode=}") diff --git a/src/pyimpspec/analysis/drt/tr_rbf.py b/src/pyimpspec/analysis/drt/tr_rbf.py index f8855cb..47a7f61 100644 --- a/src/pyimpspec/analysis/drt/tr_rbf.py +++ b/src/pyimpspec/analysis/drt/tr_rbf.py @@ -22,7 +22,7 @@ # - 10.1016/j.electacta.2015.03.123 # - 10.1016/j.electacta.2017.07.050 # Based on code from https://github.com/ciuccislab/pyDRTtools. -# pyDRTtools commit: 3694b9b4cef9b29d623bef7300280810ec351d46 +# pyDRTtools commit: 1653298d52183c36ec941197ae59399b9dc85579 from dataclasses import dataclass from time import time @@ -346,6 +346,11 @@ def _generate_truncated_multivariate_gaussians( R: NDArray[float64] = cholesky(M) R = R.T # change the lower matrix to upper matrix + # Symmetrize the matrix M + M = 0.5 * (M + M.T) + if not _is_positive_definite(M): + M = _nearest_positive_definite(M) + mu: NDArray[float64] if cov: # Using M as a covariance matrix mu = mu_r @@ -1206,8 +1211,6 @@ def _x_to_gamma( def _prepare_complex_matrices( A_re: NDArray[float64], A_im: NDArray[float64], - b_re: NDArray[float64], - b_im: NDArray[float64], M: NDArray[float64], f: Frequencies, num_freqs: int, @@ -1256,8 +1259,6 @@ def _prepare_complex_matrices( def _prepare_real_matrices( A_re: NDArray[float64], A_im: NDArray[float64], - b_re: NDArray[float64], - b_im: NDArray[float64], M: NDArray[float64], num_freqs: int, num_taus: int, @@ -1307,8 +1308,6 @@ def _prepare_real_matrices( def _prepare_imaginary_matrices( A_re: NDArray[float64], A_im: NDArray[float64], - b_re: NDArray[float64], - b_im: NDArray[float64], M: NDArray[float64], f: Frequencies, num_freqs: int, @@ -1750,7 +1749,7 @@ def calculate_drt_tr_rbf( Returns ------- - TRRBFResult + |TRRBFResult| """ global _SOLVER_IMPORTED @@ -1901,8 +1900,6 @@ def calculate_drt_tr_rbf( A_re, A_im, M, num_RL = _prepare_complex_matrices( A_re, A_im, - b_re, - b_im, M, f, num_freqs, @@ -1913,8 +1910,6 @@ def calculate_drt_tr_rbf( A_re, A_im, M, num_RL = _prepare_real_matrices( A_re, A_im, - b_re, - b_im, M, num_freqs, num_taus, @@ -1923,8 +1918,6 @@ def calculate_drt_tr_rbf( A_re, A_im, M, num_RL = _prepare_imaginary_matrices( A_re, A_im, - b_re, - b_im, M, f, num_freqs, diff --git a/src/pyimpspec/analysis/fitting.py b/src/pyimpspec/analysis/fitting.py index f7b0d1b..cc1b775 100644 --- a/src/pyimpspec/analysis/fitting.py +++ b/src/pyimpspec/analysis/fitting.py @@ -21,7 +21,9 @@ from dataclasses import dataclass from multiprocessing import Pool from multiprocessing.context import TimeoutError as MPTimeoutError +from types import SimpleNamespace from typing import ( + Any, Callable, Dict, List, @@ -68,11 +70,43 @@ from pyimpspec.typing.helpers import ( _is_boolean, _is_complex_array, + _is_floating, _is_floating_array, _is_integer, ) +class FitIdentifiers(SimpleNamespace): + """ + An object that maps an |Element|'s parameters to the name/identifier used by the corresponding `lmfit Parameter `_. The values can be accessed in two ways: + + - ``name: str = identifiers.R`` + - ``name: str = identifiers["R"]`` + + In this example, the |FitIdentifiers| instance ``identifiers`` was generated based on a |Resistor| instance. Thus, ``identifiers`` has a property called ``R``. + + It is also possible to do the following as if ``identifiers`` was a dictionary: + + - ``{symbol: identifiers[symbol] for symbol in identifiers}`` + - ``{symbol: name for symbol, name in identifiers.items()}`` + """ + + def __iter__(self): + return iter(self.__dict__) + + def items(self): + return self.__dict__.items() + + def keys(self): + return self.__dict__.keys() + + def values(self): + return self.__dict__.values() + + def __getitem__(self, key: str): + return getattr(self, key) + + @dataclass(frozen=True) class FittedParameter: """ @@ -172,13 +206,13 @@ class FitResult: Parameters ---------- - circuit: Circuit + circuit: |Circuit| The fitted circuit. - parameters: Dict[str, Dict[str, FittedParameter]] + parameters: Dict[str, Dict[str, |FittedParameter|]] Fitted parameters and their estimated standard errors (if possible to estimate). - minimizer_result: |MinimizerResult| + minimizer_result: `lmfit.MinimizerResult`_ The results of the fit as provided by the `lmfit.minimize`_ function. frequencies: |Frequencies| @@ -355,14 +389,14 @@ def get_parameters(self) -> Dict[str, Dict[str, FittedParameter]]: ------- Dict[str, Dict[str, FittedParameter]] """ - # TODO: Deprecated or unimplemented? + return self.parameters def to_parameters_dataframe( self, running: bool = False, ) -> "DataFrame": # noqa: F821 """ - Get the fitted parameters and the corresponding estimated errors as a |DataFrame| object. + Get the fitted parameters and the corresponding estimated errors as a `pandas.DataFrame`_ object. Parameters ---------- running: bool, optional @@ -370,7 +404,7 @@ def to_parameters_dataframe( Returns ------- - |DataFrame| + `pandas.DataFrame`_ """ from pandas import DataFrame @@ -402,8 +436,8 @@ def to_parameters_dataframe( parameters = self.parameters[element_name] parameter_label: str - parameter: FittedParameter - for parameter_label, parameter in parameters.items(): + for parameter_label in sorted(parameters.keys()): + parameter: FittedParameter = parameters[parameter_label] element_names.append( self.circuit.get_element_name( element, @@ -416,8 +450,12 @@ def to_parameters_dataframe( fitted_values.append(parameter.value) stderr_values.append( - parameter.stderr / parameter.value * 100 - if parameter.stderr is not None and not parameter.fixed + abs(parameter.stderr / parameter.value * 100) + if ( + parameter.stderr is not None + and not parameter.fixed + and parameter.value != 0.0 + ) else nan ) fixed.append("Yes" if parameter.fixed else "No") @@ -436,11 +474,11 @@ def to_parameters_dataframe( def to_statistics_dataframe(self) -> "DataFrame": # noqa: F821 """ - Get the statistics related to the fit as a |DataFrame| object. + Get the statistics related to the fit as a `pandas.DataFrame`_ object. Returns ------- - |DataFrame| + `pandas.DataFrame`_ """ from lmfit.minimizer import MinimizerResult from pandas import DataFrame @@ -468,67 +506,117 @@ def to_statistics_dataframe(self) -> "DataFrame": # noqa: F821 def _to_lmfit( - identifiers: Dict[int, Element], + identifiers: Dict[Element, FitIdentifiers], + constraint_expressions: Dict[str, str], + constraint_variables: Dict[str, Dict[str, Any]], ) -> "Parameters": # noqa: F821 from lmfit import Parameters - if not isinstance(identifiers, dict): - raise TypeError(f"Expected a dictionary instead of {identifiers=}") - result: Parameters = Parameters() - ident: int + name: str + kwargs: dict + for name, kwargs in constraint_variables.items(): + if "value" in kwargs: + value: float = kwargs["value"] + minimum: Union[int, float] = kwargs.get("min", -inf) + maximum: Union[int, float] = kwargs.get("max", inf) + + if not (minimum <= value <= maximum): + raise ValueError( + f"Expected {minimum} <= {value} <= {maximum} for {name=}" + ) + + result.add(name=name, **kwargs) + + queued_kwargs: List[Dict[str, Any]] = [] + element: Element - for ident, element in identifiers.items(): + mapping: FitIdentifiers + for element, mapping in identifiers.items(): lower_limits: Dict[str, float] = element.get_lower_limits() upper_limits: Dict[str, float] = element.get_upper_limits() fixed: Dict[str, bool] = element.are_fixed() symbol: str - value: float for symbol, value in element.get_values().items(): if not (lower_limits[symbol] <= value <= upper_limits[symbol]): raise ValueError( - f"Expected {lower_limits[symbol]} <= {value} <= {upper_limits[symbol]}" + f"Expected {lower_limits[symbol]=} <= {value} <= {upper_limits[symbol]=} for {symbol=}" ) - result.add( - f"{symbol}_{ident}", - value, + name: str = mapping[symbol] + kwargs = dict( + name=name, + value=value, min=lower_limits[symbol], max=upper_limits[symbol], vary=not fixed[symbol], ) + # Parameters with constraint expressions may need to be added in a + # specific order since a NameError may be raised if the expression + # makes use of a parameter that has not yet been defined. + if name in constraint_expressions: + kwargs["expr"] = constraint_expressions[name] + queued_kwargs.append(kwargs) + else: + result.add(**kwargs) + + # Brute force any queued parameters with constraint expressions until + # either all of them slot nicely into place or we run out of attempts. + num_attempts_left: int = len(queued_kwargs) + + while num_attempts_left > 0: + kwargs = queued_kwargs.pop(0) + try: + result.add(**kwargs) + num_attempts_left = len(queued_kwargs) + except NameError: + queued_kwargs.append(kwargs) + num_attempts_left -= 1 + + if len(queued_kwargs) > 0: + raise FittingError( + "Failed to generate parameters for lmfit. Make sure that all" + + " constraint expressions are valid and all the required" + + " variables have been defined!" + ) + return result def _from_lmfit( parameters: "Parameters", # noqa: F821 - identifiers: Dict[int, Element], + identifiers: Dict[Element, Dict[str, str]], ): from lmfit import Parameters if not isinstance(parameters, Parameters): raise TypeError(f"Expected a Parameters instead of {parameters=}") - if not isinstance(identifiers, dict): - raise TypeError(f"Expected a dictionary instead of {identifiers=}") + lookup: Dict[str, Tuple[Element, str]] = {} + fitted_values: Dict[Element, Dict[str, float]] = {} - result: Dict[int, Dict[str, float]] = {_: {} for _ in identifiers} + element: Element + mapping: Dict[str, str] + for element, mapping in identifiers.items(): + symbol: str + key: str + for symbol, key in mapping.items(): + lookup[key] = (element, symbol) + fitted_values[element] = {} - key: str value: float for key, value in parameters.valuesdict().items(): - ident: int - symbol: str - symbol, ident = key.rsplit("_", 1) # type: ignore - ident = int(ident) - result[ident][symbol] = float(value) + if key not in lookup: + continue - element: Element - for ident, element in identifiers.items(): - element.set_values(**result[ident]) + element, symbol = lookup[key] + fitted_values[element][symbol] = value + + for element, values in fitted_values.items(): + element.set_values(**values) def _residual( @@ -627,6 +715,7 @@ def _extract_parameters( from lmfit.minimizer import MinimizerResult parameters: Dict[str, Dict[str, FittedParameter]] = {} + internal_identifiers: Dict[int, Element] = { v: k for k, v in circuit.generate_element_identifiers(running=True).items() } @@ -691,15 +780,28 @@ def _fit_process( from lmfit import minimize from lmfit.minimizer import MinimizerResult - circuit: Circuit + original_circuit: Circuit f: Frequencies Z_exp: ComplexImpedances method: str weight: str max_nfev: int auto: bool - circuit, f, Z_exp, method, weight, max_nfev, auto = args + constraint_expressions: Dict[str, str] + constraint_variables: Dict[str, Dict[str, Any]] + ( + original_circuit, + f, + Z_exp, + method, + weight, + max_nfev, + auto, + constraint_expressions, + constraint_variables, + ) = args + circuit = deepcopy(original_circuit) if not isinstance(circuit, Circuit): raise TypeError(f"Expected a Circuit instead of {circuit=}") @@ -721,10 +823,10 @@ def _fit_process( if not _is_boolean(auto): raise TypeError(f"Expected a boolean instead of {auto=}") + identifiers: Dict[Element, FitIdentifiers] + identifiers = generate_fit_identifiers(circuit) + weight_func: Callable = _WEIGHT_FUNCTIONS[weight] - identifiers: Dict[int, Element] = { - v: k for k, v in circuit.generate_element_identifiers(running=True).items() - } with catch_warnings(): if auto: @@ -734,7 +836,7 @@ def _fit_process( try: fit: MinimizerResult = minimize( _residual, - _to_lmfit(identifiers), + _to_lmfit(identifiers, constraint_expressions, constraint_variables), method, args=( circuit, @@ -746,7 +848,7 @@ def _fit_process( max_nfev=None if max_nfev < 1 else max_nfev, ) - except (Exception, Warning): # TODO + except (Exception, Warning): # TODO: Be more specific return ( circuit, inf, @@ -799,6 +901,34 @@ def validate_circuit(circuit: Circuit): element_names.add(name) +def generate_fit_identifiers(circuit: Circuit) -> Dict[Element, FitIdentifiers]: + """ + Generate the identifiers that are provided to `lmfit`_ when fitting the provided circuit. These identifiers can be used to define optional `constraints `_ when calling |fit_circuit|. + + Parameters + ---------- + circuit: |Circuit| + The circuit to process. + + Returns + ------- + Dict[|Element|, |FitIdentifiers|] + A dictionary that maps each element to a |FitIdentifiers| instance that maps each parameter of that element to the identifier that is used by `lmfit`_ during fitting. + """ + identifiers: Dict[Element, Dict[str, str]] = {} + + element: Element + ident: int + for element, ident in circuit.generate_element_identifiers(running=True).items(): + identifiers[element] = {} + + symbol: str + for symbol in element.get_values().keys(): + identifiers[element][symbol] = f"{symbol}_{ident}" + + return {element: FitIdentifiers(**mapping) for element, mapping in identifiers.items()} + + def fit_circuit( circuit: Circuit, data: DataSet, @@ -807,16 +937,19 @@ def fit_circuit( max_nfev: int = -1, num_procs: int = -1, timeout: int = 0, + constraint_expressions: Optional[Dict[str, str]] = None, + constraint_variables: Optional[Dict[str, dict]] = None, + **kwargs, ) -> FitResult: """ Fit a circuit to a data set. Parameters ---------- - circuit: Circuit + circuit: |Circuit| The circuit to fit to a data set. - data: DataSet + data: |DataSet| The data set that the circuit will be fitted to. method: str, optional @@ -843,9 +976,17 @@ def fit_circuit( The amount of time in seconds that a single fit is allowed to take before being timed out. If this values is less than one, then no time limit is imposed. + constraint_expressions: Optional[Dict[str, str]], optional + A mapping of the optional `constraints `_ to apply. + + constraint_variables: Optional[Dict[str, dict]], optional + A mapping of variable names to their keyword arguments (see `lmfit.Parameters `_ for more information). These variables may be referenced in the optional `constraints `_. + + **kwargs + Returns ------- - FitResult + |FitResult| """ from lmfit.minimizer import MinimizerResult @@ -883,6 +1024,24 @@ def fit_circuit( f"Expected an integer equal to or greater than zero instead of {timeout=}" ) + if constraint_expressions is None: + pass + elif not isinstance(constraint_expressions, dict): + raise TypeError(f"Expected either None or a dictionary instead of {constraint_expressions=}") + elif not all(map(lambda k: isinstance(k, str), constraint_expressions.keys())): + raise TypeError(f"Expected all keys to be strings instead of {constraint_expressions=}") + elif not all(map(lambda v: isinstance(v, str), constraint_expressions.values())): + raise TypeError(f"Expected all values to be strings instead of {constraint_expressions=}") + + if constraint_variables is None: + pass + elif not isinstance(constraint_variables, dict): + raise TypeError(f"Expected either None or a dictionary instead of {constraint_variables=}") + elif not all(map(lambda k: isinstance(k, str), constraint_variables.keys())): + raise TypeError(f"Expected all keys to be strings instead of {constraint_variables=}") + elif not all(map(lambda v: isinstance(v, dict), constraint_variables.values())): + raise TypeError(f"Expected all values to be dictionaries instead of {constraint_variables=}") + num_steps: int = (len(_METHODS) if method == "auto" else 1) * ( len(_WEIGHT_FUNCTIONS) if weight == "auto" else 1 ) @@ -915,15 +1074,20 @@ def fit_circuit( ) Z_exp: ComplexImpedances = data.get_impedances() + if len(Z_exp) != len(f): + raise ValueError(f"Expected {len(Z_exp)=} == {len(f)=}") + args = ( ( - deepcopy(circuit), + circuit, f, Z_exp, method, weight, max_nfev, True, + constraint_expressions or {}, + constraint_variables or {}, ) for (method, weight) in method_weight_combos ) diff --git a/src/pyimpspec/analysis/kramers_kronig/algorithms/representation.py b/src/pyimpspec/analysis/kramers_kronig/algorithms/representation.py index dd2201b..a03c212 100644 --- a/src/pyimpspec/analysis/kramers_kronig/algorithms/representation.py +++ b/src/pyimpspec/analysis/kramers_kronig/algorithms/representation.py @@ -442,12 +442,12 @@ def suggest( Parameters ---------- - suggestions: List[Tuple[KramersKronigResult, Dict[int, float], int, int]] + suggestions: List[Tuple[|KramersKronigResult|, Dict[int, float], int, int]] A list obtained by processing List[|KramersKronigResult|] for different representations with |suggest_num_RC| and collecting the return values. Returns ------- - Tuple[KramersKronigResult, Dict[int, float], int, int] + Tuple[|KramersKronigResult|, Dict[int, float], int, int] A tuple containing: diff --git a/src/pyimpspec/analysis/kramers_kronig/cnls.py b/src/pyimpspec/analysis/kramers_kronig/cnls.py index b09f199..df36ab8 100644 --- a/src/pyimpspec/analysis/kramers_kronig/cnls.py +++ b/src/pyimpspec/analysis/kramers_kronig/cnls.py @@ -28,8 +28,10 @@ pi, ) from pyimpspec.analysis.fitting import ( + FitIdentifiers, _from_lmfit, _to_lmfit, + generate_fit_identifiers, ) from numpy.typing import NDArray from pyimpspec.circuit.base import Element @@ -137,13 +139,17 @@ def _test_wrapper(args: tuple) -> Tuple[int, Circuit]: admittance, ) - identifiers: Dict[int, Element] = { - v: k for k, v in circuit.generate_element_identifiers(running=True).items() - } + identifiers: Dict[Element, FitIdentifiers] + identifiers = generate_fit_identifiers(circuit) + fit: MinimizerResult fit = minimize( _complex_residual, - _to_lmfit(identifiers), + _to_lmfit( + identifiers, + constraint_expressions={}, + constraint_variables={}, + ), method, args=( circuit, diff --git a/src/pyimpspec/analysis/kramers_kronig/exploratory.py b/src/pyimpspec/analysis/kramers_kronig/exploratory.py index 7cfc2d2..75aee76 100644 --- a/src/pyimpspec/analysis/kramers_kronig/exploratory.py +++ b/src/pyimpspec/analysis/kramers_kronig/exploratory.py @@ -1025,7 +1025,7 @@ def evaluate_log_F_ext( Parameters ---------- - data: DataSet + data: |DataSet| The data set to be tested. test: str, optional @@ -1074,7 +1074,7 @@ def evaluate_log_F_ext( Returns ------- - List[Tuple[float, List[KramersKronigResult], float]] + List[Tuple[float, List[|KramersKronigResult|], float]] A list of tuples containing: @@ -1365,7 +1365,7 @@ def perform_exploratory_kramers_kronig_tests( Parameters ---------- - data: DataSet + data: |DataSet| The data set to be tested. test: str, optional @@ -1416,7 +1416,7 @@ def perform_exploratory_kramers_kronig_tests( Returns ------- - Tuple[List[KramersKronigResult], Tuple[KramersKronigResult, Dict[int, float], int, int]] + Tuple[List[|KramersKronigResult|], Tuple[|KramersKronigResult|, Dict[int, float], int, int]] A tuple containing a list of |KramersKronigResult| and the corresponding result of |suggest_num_RC| for the suggested extension of the range of time constants and the suggested representation of the immittance spectrum to test. """ diff --git a/src/pyimpspec/analysis/kramers_kronig/result.py b/src/pyimpspec/analysis/kramers_kronig/result.py index d68057b..62aa072 100644 --- a/src/pyimpspec/analysis/kramers_kronig/result.py +++ b/src/pyimpspec/analysis/kramers_kronig/result.py @@ -603,7 +603,7 @@ def to_statistics_dataframe( extended_statistics: int = 3, ) -> "DataFrame": # noqa: F821 r""" - Get the statistics related to the test as a |DataFrame| object. + Get the statistics related to the test as a `pandas.DataFrame`_ object. Parameters ---------- @@ -626,7 +626,7 @@ def to_statistics_dataframe( Returns ------- - |DataFrame| + `pandas.DataFrame`_ """ from pandas import DataFrame diff --git a/src/pyimpspec/analysis/kramers_kronig/single.py b/src/pyimpspec/analysis/kramers_kronig/single.py index b10fbbe..35e4e54 100644 --- a/src/pyimpspec/analysis/kramers_kronig/single.py +++ b/src/pyimpspec/analysis/kramers_kronig/single.py @@ -67,7 +67,7 @@ def perform_kramers_kronig_test( Parameters ---------- - data: DataSet + data: |DataSet| The data set to be tested. test: str, optional @@ -134,7 +134,7 @@ def perform_kramers_kronig_test( Returns ------- - KramersKronigResult + |KramersKronigResult| A single linear Kramers-Kronig test result representing the suggested extension of the range of time constants, the suggested number of RC elements (i.e., time constants), and the suggested representation of the immittance spectrum to test. """ diff --git a/src/pyimpspec/analysis/zhit/__init__.py b/src/pyimpspec/analysis/zhit/__init__.py index f44f4f5..20957d8 100644 --- a/src/pyimpspec/analysis/zhit/__init__.py +++ b/src/pyimpspec/analysis/zhit/__init__.py @@ -182,11 +182,11 @@ def get_residuals_data( def to_statistics_dataframe(self) -> "DataFrame": # noqa: F821 """ - Get the statistics related to the modulus reconstruction as a |DataFrame| object. + Get the statistics related to the modulus reconstruction as a `pandas.DataFrame`_ object. Returns ------- - |DataFrame| + `pandas.DataFrame`_ """ from pandas import DataFrame @@ -234,7 +234,7 @@ def perform_zhit( Parameters ---------- - data: DataSet + data: |DataSet| The data set for which the modulus of the impedance should be reconstructed. smoothing: str, optional @@ -274,7 +274,7 @@ def perform_zhit( Returns ------- - ZHITResult + |ZHITResult| """ if not isinstance(smoothing, str): raise TypeError(f"Expected a string instead of {smoothing=}") diff --git a/src/pyimpspec/circuit/__init__.py b/src/pyimpspec/circuit/__init__.py index 4bc09e7..ca93251 100644 --- a/src/pyimpspec/circuit/__init__.py +++ b/src/pyimpspec/circuit/__init__.py @@ -56,7 +56,7 @@ def parse_cdc(cdc: str) -> Circuit: Returns ------- - Circuit + |Circuit| """ if not isinstance(cdc, str): raise TypeError(f"Expected a string instead of {cdc=}") @@ -74,10 +74,10 @@ def simulate_spectrum( Parameters ---------- - circuit: Circuit + circuit: |Circuit| The circuit to use when calculating impedances at various frequencies. - frequencies: Optional[Frequencies], optional + frequencies: Optional[|Frequencies|], optional Excitation frequencies in hertz. If no frequencies are provided, then a frequency range of 10 mHz to 100 kHz with 10 points per decade will be used. @@ -86,7 +86,7 @@ def simulate_spectrum( Returns ------- - DataSet + |DataSet| """ if not isinstance(circuit, Circuit): raise TypeError(f"Expected a Circuit instead of {circuit=}") diff --git a/src/pyimpspec/circuit/circuit_builder.py b/src/pyimpspec/circuit/circuit_builder.py index d210997..26c7014 100644 --- a/src/pyimpspec/circuit/circuit_builder.py +++ b/src/pyimpspec/circuit/circuit_builder.py @@ -60,7 +60,7 @@ def series(self) -> "CircuitBuilder": Returns ------- - CircuitBuilder + |CircuitBuilder| """ series: "CircuitBuilder" = CircuitBuilder(parallel=False) self._elements.append(series) @@ -73,7 +73,7 @@ def parallel(self) -> "CircuitBuilder": Returns ------- - CircuitBuilder + |CircuitBuilder| """ parallel: "CircuitBuilder" = CircuitBuilder(parallel=True) self._elements.append(parallel) @@ -91,7 +91,7 @@ def add(self, element: Element): Parameters ---------- - element: Element + element: |Element| The element to add to the current series or parallel connection. """ if not isinstance(element, Element): @@ -135,6 +135,6 @@ def to_circuit(self) -> Circuit: Returns ------- - Circuit + |Circuit| """ return Parser().process(self._to_string()) diff --git a/src/pyimpspec/circuit/diagrams/circuitikz.py b/src/pyimpspec/circuit/diagrams/circuitikz.py index 1d2c970..7748799 100644 --- a/src/pyimpspec/circuit/diagrams/circuitikz.py +++ b/src/pyimpspec/circuit/diagrams/circuitikz.py @@ -74,7 +74,7 @@ def to_circuitikz( running: bool, optional Whether or not to use running counts as the lower indices of elements. - custom_labels: Optional[Dict[Element, str]], optional + custom_labels: Optional[Dict[|Element|, str]], optional A mapping of elements to their custom labels that are used instead of the automatically generated labels. The labels can make use of LaTeX's math mode. diff --git a/src/pyimpspec/circuit/diagrams/schemdraw.py b/src/pyimpspec/circuit/diagrams/schemdraw.py index 43fa2e6..036c0db 100644 --- a/src/pyimpspec/circuit/diagrams/schemdraw.py +++ b/src/pyimpspec/circuit/diagrams/schemdraw.py @@ -70,7 +70,7 @@ def to_drawing( running: bool, optional Whether or not to use running counts as the lower indices of elements. - custom_labels: Optional[Dict[Element, str]], optional + custom_labels: Optional[Dict[|Element|, str]], optional A mapping of elements to their custom labels that are used instead of the automatically generated labels. The labels can make use of LaTeX's math mode. diff --git a/src/pyimpspec/circuit/kramers_kronig.py b/src/pyimpspec/circuit/kramers_kronig.py index 5162815..0a401ab 100644 --- a/src/pyimpspec/circuit/kramers_kronig.py +++ b/src/pyimpspec/circuit/kramers_kronig.py @@ -68,6 +68,7 @@ def _impedance(self, f: Frequencies, R: float, tau: float) -> ComplexImpedances: ), ], ), + private=True, ) @@ -105,4 +106,5 @@ def _impedance(self, f: Frequencies, C: float, tau: float) -> ComplexImpedances: ), ], ), + private=True, ) diff --git a/src/pyimpspec/circuit/parser.py b/src/pyimpspec/circuit/parser.py index 6410010..9ba95c6 100644 --- a/src/pyimpspec/circuit/parser.py +++ b/src/pyimpspec/circuit/parser.py @@ -88,7 +88,7 @@ class Parser: def __init__(self): self._tokens: List[Token] = [] self._stack: List[Stackable] = [] - self._valid_elements: Dict[str, Type[Element]] = get_elements() + self._valid_elements: Dict[str, Type[Element]] = get_elements(private=True) if not isinstance(self._valid_elements, dict): raise TypeError(f"Expected a dictionary instead of {self._valid_elements=}") diff --git a/src/pyimpspec/circuit/registry.py b/src/pyimpspec/circuit/registry.py index 5357491..5d96eb6 100644 --- a/src/pyimpspec/circuit/registry.py +++ b/src/pyimpspec/circuit/registry.py @@ -59,6 +59,7 @@ _ELEMENTS: Dict[str, Type[Element]] = {} _DEFAULT_ELEMENTS: Dict[str, Type[Element]] = {} _DEFAULT_ELEMENT_PARAMETERS: Dict[str, Dict[str, float]] = {} +_PRIVATE_ELEMENTS: Dict[str, Type[Element]] = {} def _initialized(): @@ -138,7 +139,7 @@ class SubcircuitDefinition: description: str A description of the subcircuit. Can be an empty string. - value: Optional[Connection] + value: Optional[|Connection|] The default value for the parameter. Can be a connection such as Series or Parallel, or None. """ @@ -156,7 +157,7 @@ class ElementDefinition: Parameters ---------- - Class: Type[Element] + Class: Type[|Element|] The class representing the circuit element. symbol: str @@ -177,7 +178,7 @@ class ElementDefinition: For example, "(2*pi*f*C*I)^-1" for a capacitor is valid and will show up in the generated docstring as :math:`Z = \frac{1}{j 2 \pi f C}`. This is also used to verify that the circuit element's ``_impedance`` method outputs the same answer as the circuit element's impedance equation. - parameters: List[ParameterDefinition] + parameters: List[|ParameterDefinition|] A list of objects that define the circuit element's parameter(s). """ @@ -507,7 +508,7 @@ def _initialize_element( return (symbol, Class) -def get_elements(default_only: bool = False) -> Dict[str, Type[Element]]: +def get_elements(default_only: bool = False, private: bool = False) -> Dict[str, Type[Element]]: """ Returns a dictionary that maps element symbols to their corresponding classes. @@ -516,16 +517,21 @@ def get_elements(default_only: bool = False) -> Dict[str, Type[Element]]: default_only: bool, optional Return only the elements that are included by default in pyimpspec (i.e., exclude user-defined elements). + private: bool, optional + Include elements that have been marked as private. + Returns ------- - Dict[str, Type[Element]] + Dict[str, Type[|Element|]] """ if not _is_boolean(default_only): raise TypeError(f"Expected a boolean instead of {default_only=}") keys: List[str] = sorted(list((_DEFAULT_ELEMENTS if default_only else _ELEMENTS).keys())) + if not private: + keys = [k for k in keys if k not in _PRIVATE_ELEMENTS] - elements: Dict[str, Type[Element]] = {_: _ELEMENTS[_] for _ in keys} + elements: Dict[str, Type[Element]] = {k: _ELEMENTS[k] for k in keys} return elements @@ -536,7 +542,7 @@ def register_element(definition: ElementDefinition, **kwargs): Parameters ---------- - definition: ElementDefinition + definition: |ElementDefinition| An ElementDefinition instance that defines a class, which inherits from the Element class and implements a method called `_impedance` that takes a float value (i.e., the excitation frequency in hertz) and returns a complex value. **kwargs @@ -554,6 +560,8 @@ def register_element(definition: ElementDefinition, **kwargs): ) _ELEMENTS[symbol] = Class + if kwargs.get("private", False) is True: + _PRIVATE_ELEMENTS[symbol] = Class def reset_default_parameter_values(elements: Optional[Union[Type[Element], List[Type[Element]]]] = None): @@ -562,7 +570,7 @@ def reset_default_parameter_values(elements: Optional[Union[Type[Element], List[ Parameters ---------- - elements: Optional[Union[Type[Element], List[Type[Element]]]], optional + elements: Optional[Union[Type[|Element|], List[Type[|Element|]]]], optional Specific element class(es) to reset. If none are provided, then all the elements that are included by default in pyimpspec are reset. """ if elements is None: @@ -589,6 +597,14 @@ def reset_default_parameter_values(elements: Optional[Union[Type[Element], List[ def reset(elements: bool = True, default_parameters: bool = True): """ Remove all user-defined elements from the registry. + + Parameters + ---------- + elements: bool, optional + If true, then remove any user-defined elements. + + default_parameters: bool, optional + If true, reset the default parameters of all elements. """ if elements: _ELEMENTS.clear() @@ -604,7 +620,7 @@ def remove_elements(elements: Union[Type[Element], List[Type[Element]]]): Parameters ---------- - elements: Union[Type[Element], List[Type[Element]]] + elements: Union[Type[|Element|], List[Type[|Element|]]] The element(s) to remove. """ if isinstance(elements, list): @@ -631,4 +647,6 @@ def remove_elements(elements: Union[Type[Element], List[Type[Element]]]): for key in list(_ELEMENTS.keys()): if _ELEMENTS[key] is element: _ELEMENTS.pop(key) + if key in _PRIVATE_ELEMENTS and _PRIVATE_ELEMENTS[key] is element: + _PRIVATE_ELEMENTS.pop(key) break diff --git a/src/pyimpspec/cli/args.py b/src/pyimpspec/cli/args.py index 42ddcdd..0d03649 100644 --- a/src/pyimpspec/cli/args.py +++ b/src/pyimpspec/cli/args.py @@ -27,6 +27,7 @@ _RBF_SHAPES, _CROSS_VALIDATION_METHODS, ) +from pyimpspec.analysis.drt.lm import _MODEL_ORDER_METHODS from pyimpspec.analysis.fitting import ( _METHODS as FIT_METHODS, _WEIGHT_FUNCTIONS as FIT_WEIGHTS, @@ -513,7 +514,7 @@ def drt_args(parser: ArgumentParser): metavar="FLOAT", type=float, default=0.5, - help="A maximum limit (between 0.0 and 1.0) for the relative vertical symmetry of the DRT. A high degree of symmetry is common for results where the gamma value oscillates wildly (e.g., due to a small regularization parameter). A low value for the limit should improve the results but may cause the BHT method to take longer to finish. The TR-RBF method only uses this limit when the regularization parameter (lambda) is not provided.", + help="A maximum limit (between 0.0 and 1.0) for the relative vertical symmetry of the DRT. A high degree of symmetry is common for results where the gamma value oscillates wildly (e.g., due to a small regularization parameter). A low value for the limit should improve the results but may cause the BHT method to take longer to finish.", ) parser.add_argument( "--circuit", @@ -558,6 +559,26 @@ def drt_args(parser: ArgumentParser): type=int, help="The maximum number of function evaluations to use when fitting a circuit (m(RQ)fit method only). A value below 1 means no limit. Defaults to -1.", ) + parser.add_argument( + "--model-order", + "-k", + dest="model_order", + metavar="INTEGER", + type=int, + default=0, + help="The model order (k) to use (Loewner method only). Defaults to 0.", + ) + parser.add_argument( + "--model-order-method", + "-km", + dest="model_order_method", + metavar="STRING", + type=str, + default="matrix_rank", + help="The approach to use to automatically pick the model order (k) if a model order is not specified (Loewner method only). " + + ", ".join(sorted(map(lambda _: f"'{_}'", _MODEL_ORDER_METHODS))) + + ". Defaults to 'matrix_rank'.", + ) parser.add_argument( "--num-procs", dest="num_procs", @@ -571,6 +592,38 @@ def drt_args(parser: ArgumentParser): Defaults to -1. """.strip(), ) + parser.add_argument( + "--analyze-peaks", + "-ap", + dest="analyze_peaks", + action="store_true", + help="Perform peak analyses by fitting skew normal distributions.", + ) + parser.add_argument( + "--num-peaks", + "-np", + dest="num_peaks", + type=int, + default=0, + help="The number of peaks to include in the peak analysis. The tallest peaks are prioritized. Only applicable when peak positions are not provided manually.", + ) + parser.add_argument( + "--peak-positions", + "-pp", + dest="peak_positions", + metavar="FLOAT", + type=float, + nargs="*", + default=[], + help="The positions of the peaks to analyze. If not provided, then peaks and their positions are detected automatically.", + ) + parser.add_argument( + "--disallow-skew", + "-ds", + dest="disallow_skew", + action="store_true", + help="If true, then normal distributions are used instead of skew normal distributions when analyzing peaks.", + ) add_output_args(parser) add_plot_args(parser) parser.add_argument( @@ -582,6 +635,13 @@ def drt_args(parser: ArgumentParser): default="drt", help="The type of plot to generate: 'bode', 'drt', 'imaginary', 'magnitude', 'nyquist', 'phase', 'real', or 'real-imaginary'. Defaults to 'drt'.", ) + parser.add_argument( + "--plot-frequency", + "-pF", + dest="plot_frequency", + action="store_true", + help="Plot gamma vs frequency instead of time constant.", + ) def fit_args(parser: ArgumentParser): diff --git a/src/pyimpspec/cli/drt.py b/src/pyimpspec/cli/drt.py index 449a098..3f1ca63 100644 --- a/src/pyimpspec/cli/drt.py +++ b/src/pyimpspec/cli/drt.py @@ -29,6 +29,7 @@ Dict, IO, List, + Optional, Tuple, ) from .utility import ( @@ -53,6 +54,7 @@ def overlay_plot( mpl, parse_cdc, ) + from pyimpspec.analysis.drt import LMResult if not (args.output_name[0] != "" or not args.output): raise ValueError("Expected an output name") @@ -92,6 +94,8 @@ def overlay_plot( gaussian_width=args.gaussian_width, num_per_decade=args.num_per_decade, max_nfev=args.max_nfev, + model_order=args.model_order, + model_order_method=args.model_order_method, num_procs=args.num_procs, ) drts.append( @@ -142,6 +146,7 @@ def overlay_plot( colors={"gamma": color}, figure=figure, axes=axes, + frequency=args.plot_frequency, ) if not args.plot_no_legend: @@ -180,13 +185,16 @@ def individual_plots( args: Namespace, print_func: Callable, ): + from pandas import DataFrame from pyimpspec import ( + DRTPeaks, DRTResult, DataSet, calculate_drt, mpl, parse_cdc, ) + from pyimpspec.analysis.drt import LMResult from pyimpspec.plot.colors import COLOR_BLACK from pyimpspec.plot.mpl.helpers import _color_axis @@ -247,8 +255,23 @@ def individual_plots( gaussian_width=args.gaussian_width, num_per_decade=args.num_per_decade, max_nfev=args.max_nfev, + model_order=args.model_order, + model_order_method=args.model_order_method, num_procs=args.num_procs, ) + + peaks: Optional[DRTPeaks] = None + if ( + args.analyze_peaks + and hasattr(drt, "analyze_peaks") + and not isinstance(drt, LMResult) + ): + peaks = drt.analyze_peaks( + num_peaks=args.num_peaks, + peak_positions=args.peak_positions or None, + disallow_skew=args.disallow_skew, + ) + clear_default_handler_output() label: str = f"{data.get_label()}\n{drt.get_label()}" @@ -257,6 +280,7 @@ def individual_plots( drt, title=label if args.plot_title else "", peak_threshold=args.peak_threshold, + frequency=args.plot_frequency, **kwargs, ) else: @@ -265,6 +289,7 @@ def individual_plots( data=data, title=label if args.plot_title else "", peak_threshold=args.peak_threshold, + frequency=args.plot_frequency, **kwargs, ) @@ -339,13 +364,25 @@ def individual_plots( if args.peak_threshold >= 0.0: fragments.append( format_text( - drt.to_peaks_dataframe(threshold=args.peak_threshold), args + drt.to_peaks_dataframe(threshold=args.peak_threshold), + args, ) ) if args.method == "bht": fragments.append(format_text(drt.to_scores_dataframe(), args)) + if ( + args.analyze_peaks + and peaks is not None + and not isinstance(drt, LMResult) + ): + if isinstance(peaks, tuple): + for df in (p.to_peaks_dataframe() for p in peaks): + fragments.append(format_text(df, args)) + else: + fragments.append(format_text(peaks.to_peaks_dataframe(), args)) + report: str = "\n\n".join(fragments) if args.output: output_path: str = get_output_path( diff --git a/src/pyimpspec/cli/utility.py b/src/pyimpspec/cli/utility.py index 2779161..0213fa6 100644 --- a/src/pyimpspec/cli/utility.py +++ b/src/pyimpspec/cli/utility.py @@ -263,7 +263,7 @@ def _parse_identity(identity: str) -> Tuple[str, dict]: i: int = identity.rfind(":") if ":" in identity and i > max(map(lambda s: identity.rfind(s), ("}", "]", ")"))): - for arg in identity[i + 1:].split(","): + for arg in map(str.strip, identity[i + 1:].split(",")): key, value = arg.split("=") if key in kwarg_types: kwargs[key] = kwarg_types[key](value) diff --git a/src/pyimpspec/data/__init__.py b/src/pyimpspec/data/__init__.py index 550e091..c538c17 100644 --- a/src/pyimpspec/data/__init__.py +++ b/src/pyimpspec/data/__init__.py @@ -136,7 +136,7 @@ def parse_data( Returns ------- - List[DataSet] + List[|DataSet|] """ _validate_path(path) if not (isinstance(file_format, str) or file_format is None): diff --git a/src/pyimpspec/data/data_set.py b/src/pyimpspec/data/data_set.py index 328d7ea..23d43db 100644 --- a/src/pyimpspec/data/data_set.py +++ b/src/pyimpspec/data/data_set.py @@ -689,7 +689,7 @@ def to_dataframe( negative_phase: bool = False, ) -> "DataFrame": # noqa: F821 """ - Create a |DataFrame| instance from this DataSet. + Create a `pandas.DataFrame`_ instance from this DataSet. Parameters ---------- @@ -711,7 +711,7 @@ def to_dataframe( Returns ------- - |DataFrame| + `pandas.DataFrame`_ """ from pandas import DataFrame @@ -963,12 +963,12 @@ def dataframe_to_data_sets( degrees: bool = True, ) -> List[DataSet]: """ - Takes a |DataFrame| object, checks if there are multiple frequency sweeps, and, if necessary, splits the data into multiple DataSet objects. + Takes a `pandas.DataFrame`_ object, checks if there are multiple frequency sweeps, and, if necessary, splits the data into multiple DataSet objects. Parameters ---------- - df: DataFrame - The DataFrame to be converted. + df: `pandas.DataFrame`_ + The `pandas.DataFrame`_ to be converted. The object should contain columns for frequencies and either real and imaginary impedances, or the magnitudes and phases of impedances. Multiple labels are supported for the column headers and they are detected based on whether or not the header starts with one of the supported labels. The comparisons are case insensitive and inverted signs are detected based on if the column header starts with either a hyphen (-) or a minus sign (:math:`-`). @@ -980,7 +980,7 @@ def dataframe_to_data_sets( - Phases of the impedances: 'phase', 'phz', and 'phi'. path: Union[str, pathlib.Path] - The path to the file that was used to create the DataFrame. + The path to the file that was used to create the `pandas.DataFrame`_. label: str, optional The label assigned to the new DataSet. @@ -990,7 +990,7 @@ def dataframe_to_data_sets( Returns ------- - List[DataSet] + List[|DataSet|] """ from pandas import DataFrame diff --git a/src/pyimpspec/mock_data.py b/src/pyimpspec/mock_data.py index 36a30e0..0881184 100644 --- a/src/pyimpspec/mock_data.py +++ b/src/pyimpspec/mock_data.py @@ -114,9 +114,12 @@ def generate_circuit(self, **kwargs) -> Circuit: drift: float = kwargs.get("drift", 1.0) if not _is_floating(drift): - raise TypeError(f"Expected a float instead of {drift=}") - else: - drift *= self.drift + try: + drift = float(drift) + except ValueError: + raise TypeError(f"Expected a float instead of {drift=}") + + drift *= self.drift connection: Connection for connection in circuit.get_connections(recursive=True): @@ -322,6 +325,36 @@ def get_identifier(self) -> str: num_per_decade=10, drift=0.5, ), + MockDefinition( + label="Circuit 13", + cdc="R{R=70}(R{R=200}C{C=2.5e-3})(R{R=100}C{C=1e-4})(R{R=50}L{L=3e2})", + log_max_f=2.0, + log_min_f=-2.0, + num_per_decade=10, + ), + MockDefinition( + label="Circuit 13 invalid", + cdc="R{R=70}(R{R=200}C{C=2.5e-3})(R{R=100}C{C=1e-4})(R{R=50:drift}L{L=3e2})", + log_max_f=2.0, + log_min_f=-2.0, + num_per_decade=10, + drift=0.5, + ), + MockDefinition( + label="Circuit 14", + cdc="R{R=70}(R{R=200}Q{Y=2.5e-3,n=0.9})(R{R=100}Q{Y=1e-4,n=0.85})(R{R=50}La{L=3e2,n=0.95})", + log_max_f=2.0, + log_min_f=-2.0, + num_per_decade=10, + ), + MockDefinition( + label="Circuit 14 invalid", + cdc="R{R=70}(R{R=200}Q{Y=2.5e-3,n=0.9})(R{R=100}Q{Y=1e-4,n=0.85})(R{R=50:drift}La{L=3e2,n=0.95})", + log_max_f=2.0, + log_min_f=-2.0, + num_per_decade=10, + drift=0.5, + ), ] if len(set(d.get_identifier() for d in _definitions)) != len(_definitions): @@ -341,7 +374,8 @@ def _add_noise( seed = seed & (2**32 - 1) # Truncate to 32 bits just to be safe rs: RandomState = RandomState(seed=seed) - Z_noisy: ComplexImpedances = zeros(data.get_num_points(), dtype=ComplexImpedance) + f: Frequencies = data.get_frequencies() + Z_noisy: ComplexImpedances = zeros(len(f), dtype=ComplexImpedance) Z_noisy.real = rs.normal(0, sd) Z_noisy.imag = rs.normal(0, sd) Z_noisy += Z_ideal @@ -377,15 +411,23 @@ def _simulate_spectrum( "log_max_f", 5.0 if definition is None else definition.log_max_f, ) + if not _is_floating(log_max_f): + try: + log_max_f = float(log_max_f) + except ValueError: + raise TypeError(f"Expected a float instead of {log_max_f=}") + log_min_f: float = kwargs.get( "log_min_f", -1.0 if definition is None else definition.log_min_f, ) - if not _is_floating(log_max_f): - raise TypeError(f"Expected a float instead of {log_max_f=}") - elif not _is_floating(log_min_f): - raise TypeError(f"Expected a float instead of {log_min_f=}") - elif log_max_f <= log_min_f: + if not _is_floating(log_min_f): + try: + log_min_f = float(log_min_f) + except ValueError: + raise TypeError(f"Expected a float instead of {log_min_f=}") + + if log_max_f <= log_min_f: raise ValueError(f"Expected {log_min_f=} < {log_max_f=}") num_per_decade: int = kwargs.get( @@ -393,16 +435,24 @@ def _simulate_spectrum( 10 if definition is None else definition.num_per_decade, ) if not _is_integer(num_per_decade): - raise TypeError(f"Expected an integer instead of {num_per_decade=}") - elif num_per_decade < 1: + try: + num_per_decade = int(num_per_decade) + except ValueError: + raise TypeError(f"Expected an integer instead of {num_per_decade=}") + + if num_per_decade < 1: raise ValueError( f"Expected a value greater than zero instead of {num_per_decade=}" ) noise: float = kwargs.get("noise", 0.0) if not _is_floating(noise): - raise TypeError(f"Expected a float instead of {noise=}") - elif noise < 0.0: + try: + noise = float(noise) + except ValueError: + raise TypeError(f"Expected a float instead of {noise=}") + + if noise < 0.0: raise ValueError( f"Expected a value greater than or equal to zero instead of {noise=}" ) diff --git a/src/pyimpspec/plot/mpl/drt.py b/src/pyimpspec/plot/mpl/drt.py index f7ea25c..ce8de2e 100644 --- a/src/pyimpspec/plot/mpl/drt.py +++ b/src/pyimpspec/plot/mpl/drt.py @@ -55,6 +55,7 @@ def plot_drt( axes: Optional[List["Axes"]] = None, # noqa: F821 adjust_axes: bool = True, colored_axes: bool = False, + frequency: bool = False, **kwargs, ) -> Tuple["Figure", List["Axes"]]: # noqa: F821 """ @@ -103,6 +104,9 @@ def plot_drt( colored_axes: bool, optional Color the y-axes. + frequency: bool, optional + Plot gamma as a function of frequency. + **kwargs Returns @@ -224,6 +228,7 @@ def plot_drt( figure=figure, axes=[axes[2]], adjust_axes=adjust_axes, + frequency=frequency, ) plot_residuals( diff --git a/src/pyimpspec/plot/mpl/gamma.py b/src/pyimpspec/plot/mpl/gamma.py index db6584f..3858459 100644 --- a/src/pyimpspec/plot/mpl/gamma.py +++ b/src/pyimpspec/plot/mpl/gamma.py @@ -17,10 +17,20 @@ # The licenses of pyimpspec's dependencies and/or sources of portions of code are included in # the LICENSES folder. -from pyimpspec.analysis.drt import DRTResult +from pyimpspec.analysis.drt import ( + DRTPeaks, + DRTResult, + BHTResult, + LMResult, + MRQFitResult, + TRNNLSResult, + TRRBFResult, +) from numpy import ( complex128, isnan, + pi, + zeros, ) from pyimpspec.typing import ( Gamma, @@ -35,7 +45,11 @@ Tuple, _is_boolean, ) -from pyimpspec.plot.colors import COLOR_BLACK +from pyimpspec.plot.colors import ( + COLOR_BLACK, + COLOR_BLUE, + COLOR_RED, +) from .helpers import ( _color_axis, _configure_log_limits, @@ -51,12 +65,16 @@ def _plot_credible_intervals( color: str, bounds_alpha: float, axis: "Axes", # noqa: F821 + versus_frequency: bool, ): x: TimeConstants y1: Gammas y2: Gammas y3: Gammas x, y1, y2, y3 = drt.get_drt_credible_intervals_data() + if versus_frequency: + x = 1/(2*pi*x) + if y1.size == y2.size == y3.size > 0: mean_label: Optional[str] = None ci_label: Optional[str] = None @@ -67,6 +85,7 @@ def _plot_credible_intervals( else: mean_label = "mean" ci_label = r"$3\sigma$ CI" + axis.fill_between( x, y2, @@ -75,6 +94,7 @@ def _plot_credible_intervals( alpha=bounds_alpha, label=ci_label, ) + axis.plot( x, y1, @@ -84,12 +104,219 @@ def _plot_credible_intervals( ) +def _plot_bht_gammas( + drt: BHTResult, + peak_threshold: float, + label: str, + color: str, + axis: "Axes", # noqa: F821 + versus_frequency: bool, +): + x: TimeConstants + real_y: Gammas + imag_y: Gammas + x, real_y, imag_y = drt.get_drt_data() + + real_label: Optional[str] = None + imaginary_label: Optional[str] = None + + if label.strip() != "": + if label != "": + real_label = f"{label}, real" + imaginary_label = f"{label}, imag." + else: + real_label = "real" + imaginary_label = "imag." + + axis.plot( + x, + real_y, + color=color, + linestyle="-", + label=real_label, + ) + axis.plot( + x, + imag_y, + color=color, + linestyle=":", + label=imaginary_label, + ) + + if peak_threshold >= 0.0: + time_constants_real: TimeConstants + gammas_real: Gammas + time_constants_imag: TimeConstants + gammas_imag: Gammas + ( + time_constants_real, + gammas_real, + time_constants_imag, + gammas_imag, + ) = drt.get_peaks(threshold=peak_threshold) + + for time_constants, gammas in ( + (time_constants_real, gammas_real), + (time_constants_imag, gammas_imag), + ): + for x, y in zip(time_constants, gammas): + if isnan(x) or isnan(y): + continue + + if versus_frequency: + x = 1/(2*pi*x) + + axis.plot( + [x, x], + [0, y], + linestyle=":", + alpha=0.5, + color=color, + ) + + if hasattr(drt, "_peak_analysis"): + peaks_real: DRTPeaks + peaks_imag: DRTPeaks + peaks_real, peaks_imag = drt._peak_analysis + + for peaks in (peaks_real, peaks_imag): + peaks_color = COLOR_BLUE if peaks == peaks_real else COLOR_RED + + baseline = None + total_y = None + for i in range(0, peaks.get_num_peaks()): + x = peaks.get_time_constants(num_per_decade=100) + if versus_frequency: + x = 1/(2*pi*x) + + y = peaks.get_gammas(peak_indices=[i], num_per_decade=100) + if total_y is None: + total_y = zeros(x.shape, dtype=y.dtype) + + total_y += y + + if baseline is None: + baseline = zeros(x.shape, dtype=y.dtype) + + axis.fill_between( + x, + baseline, + y, + alpha=0.25, + ) + + if total_y is not None: + axis.plot( + x, + total_y, + linestyle="--", + color=peaks_color, + label="Sum, " + peaks.suffix + ) + + +def _plot_lm_gammas( + drt: LMResult, + peak_threshold: float, + label: str, + color: str, + axis: "Axes", # noqa: F821 + versus_frequency: bool, +): + x_RC: TimeConstants + y_RC: Gammas + x_RL: TimeConstants + y_RL: Gammas + x_RC, y_RC, x_RL, y_RL = drt.get_drt_data() + if versus_frequency: + x_RC = 1/(2*pi*x_RC) + x_RL = 1/(2*pi*x_RL) + + y_RL *= -1 + + label_RC: Optional[str] = None + label_RL: Optional[str] = None + + if label.strip() != "": + if label != "": + label_RC = f"{label}, RC" + label_RL = f"{label}, RL" + else: + label_RC = "RC" + label_RL = "RL" + + if x_RC.any(): + axis.scatter( + x_RC, + y_RC, + edgecolor=color, + facecolor="none", + marker="o", + label=label_RC, + ) + for x, y in zip(x_RC, y_RC): + axis.plot( + [x, x], + [0.0, y], + color=color, + linestyle="--", + ) + + if x_RL.any(): + axis.scatter( + x_RL, + y_RL, + color=color, + marker="x", + label=label_RL, + ) + for x, y in zip(x_RL, y_RL): + axis.plot( + [x, x], + [0.0, y], + color=color, + linestyle=":", + ) + + x: List[float] = x_RC.tolist() + x_RL.tolist() + if x: + axis.plot( + [min(x), max(x)], + [0.0, 0.0], + color=color, + linestyle="-", + ) + + def _plot_gammas( drt: DRTResult, + peak_threshold: float, label: str, color: str, axis: "Axes", # noqa: F821 + versus_frequency: bool, ): + if isinstance(drt, BHTResult): + _plot_bht_gammas( + drt, + peak_threshold=peak_threshold, + label=label, + color=color, + axis=axis, + versus_frequency=versus_frequency, + ) + return + elif isinstance(drt, LMResult): + _plot_lm_gammas( + drt, + peak_threshold=peak_threshold, + label=label, + color=color, + axis=axis, + versus_frequency=versus_frequency, + ) + return + x: TimeConstants real_y: Gammas = None imag_y: Gammas = None @@ -103,6 +330,10 @@ def _plot_gammas( imag_y = y.imag else: real_y = y + + if versus_frequency: + x = 1/(2*pi*x) + real_label: Optional[str] = None imaginary_label: Optional[str] = None if label.strip() != "": @@ -115,12 +346,14 @@ def _plot_gammas( else: real_label = "real" imaginary_label = "imag." + axis.plot( x, real_y, color=color, label=real_label, ) + if imag_y is not None and imag_y.size > 0: axis.plot( x, @@ -130,10 +363,6 @@ def _plot_gammas( label=imaginary_label, ) - -def _plot_peaks(drt, peak_threshold: float, color, axis): - x: TimeConstant - y: Gamma if peak_threshold >= 0.0 and hasattr(drt, "get_peaks") and callable(drt.get_peaks): # Most DRT results return a tuple with a length of two. However, BHTResult # returns a tuple with a length of four (real tau, real gamma, imaginary tau, @@ -142,6 +371,10 @@ def _plot_peaks(drt, peak_threshold: float, color, axis): for x, y in zip(peaks[0], peaks[1]): if isnan(x) or isnan(y): continue + + if versus_frequency: + x = 1/(2*pi*x) + axis.plot( [x, x], [0, y], @@ -149,10 +382,15 @@ def _plot_peaks(drt, peak_threshold: float, color, axis): alpha=0.5, color=color, ) + if len(peaks) == 4: for x, y in zip(peaks[2], peaks[3]): if isnan(x) or isnan(y): continue + + if versus_frequency: + x = 1/(2*pi*x) + axis.plot( [x, x], [0, y], @@ -161,6 +399,42 @@ def _plot_peaks(drt, peak_threshold: float, color, axis): color=color, ) + if hasattr(drt, "_peak_analysis"): + peaks: DRTPeaks = drt._peak_analysis + peaks_color = COLOR_BLUE + + baseline = None + total_y = None + for i in range(0, peaks.get_num_peaks()): + x = peaks.get_time_constants(num_per_decade=100) + if versus_frequency: + x = 1/(2*pi*x) + + y = peaks.get_gammas(peak_indices=[i], num_per_decade=100) + if total_y is None: + total_y = zeros(x.shape, dtype=y.dtype) + + total_y += y + + if baseline is None: + baseline = zeros(x.shape, dtype=y.dtype) + + axis.fill_between( + x, + baseline, + y, + alpha=0.25, + ) + + if total_y is not None: + axis.plot( + x, + total_y, + linestyle="--", + color=peaks_color, + label="Sum", + ) + def plot_gamma( drt: DRTResult, @@ -173,6 +447,7 @@ def plot_gamma( axes: Optional[List["Axes"]] = None, # noqa: F821 adjust_axes: bool = True, colored_axes: bool = False, + frequency: bool = False, **kwargs, ) -> Tuple["Figure", List["Axes"]]: # noqa: F821 """ @@ -211,6 +486,9 @@ def plot_gamma( colored_axes: bool, optional Color the y-axes. + frequency: bool, optional + Plot gamma as a function of frequency. + **kwargs Returns @@ -240,18 +518,34 @@ def plot_gamma( elif not isinstance(label, str): raise TypeError(f"Expected a string or None instead of {label=}") + if not _is_boolean(frequency): + raise TypeError(f"Expected a boolean instead of {frequency=}") + if hasattr(drt, "get_drt_credible_intervals_data") and callable( drt.get_drt_credible_intervals_data ): - _plot_credible_intervals(drt, label, color, bounds_alpha, axis) + _plot_credible_intervals( + drt=drt, + label=label, + color=color, + bounds_alpha=bounds_alpha, + axis=axis, + versus_frequency=frequency, + ) - _plot_peaks(drt, peak_threshold, color, axis) - _plot_gammas(drt, label, color, axis) + _plot_gammas( + drt=drt, + peak_threshold=peak_threshold, + label=label, + color=color, + axis=axis, + versus_frequency=frequency, + ) if not _is_boolean(adjust_axes): raise TypeError(f"Expected a boolean instead of {adjust_axes=}") elif adjust_axes: - axis.set_xlabel(r"$\tau\ (\rm s)$") + axis.set_xlabel(r"$f\ ({\rm Hz})$" if frequency else r"$\tau\ (\rm s)$") axis.set_ylabel(r"$\gamma\ (\Omega)$") _configure_log_scale(axis, x=True) _configure_log_limits(axis, x=True) diff --git a/src/pyimpspec/plot/mpl/imaginary.py b/src/pyimpspec/plot/mpl/imaginary.py index dc5211c..e26ff15 100644 --- a/src/pyimpspec/plot/mpl/imaginary.py +++ b/src/pyimpspec/plot/mpl/imaginary.py @@ -18,11 +18,13 @@ # the LICENSES folder. from inspect import signature +from pyimpspec.circuit.circuit import Circuit from pyimpspec.data import DataSet from pyimpspec.analysis import ( KramersKronigResult, FitResult, ) +from pyimpspec.analysis.utility import _interpolate from pyimpspec.analysis.drt import DRTResult from pyimpspec.typing import ( Frequencies, @@ -156,6 +158,12 @@ def plot_imaginary( data.get_impedances(num_per_decade=num_per_decade) ** (-1 if admittance else 1) ).imag + elif line and hasattr(data, "circuit") and isinstance(data.circuit, Circuit): + x = _interpolate(data.get_frequencies(), num_per_decade=num_per_decade) + y = ( + data.circuit.get_impedances(x) + ** (-1 if admittance else 1) + ).imag else: x = data.get_frequencies() y = (data.get_impedances() ** (-1 if admittance else 1)).imag diff --git a/src/pyimpspec/plot/mpl/kramers_kronig.py b/src/pyimpspec/plot/mpl/kramers_kronig.py index cb9e487..13e0296 100644 --- a/src/pyimpspec/plot/mpl/kramers_kronig.py +++ b/src/pyimpspec/plot/mpl/kramers_kronig.py @@ -107,7 +107,7 @@ def plot_log_F_ext( adjust_axes: bool, optional Whether or not to adjust the axes (label, scale, limits, etc.). - **kwargs, + **kwargs Returns ------- diff --git a/src/pyimpspec/plot/mpl/magnitude.py b/src/pyimpspec/plot/mpl/magnitude.py index 13a3134..e2d384d 100644 --- a/src/pyimpspec/plot/mpl/magnitude.py +++ b/src/pyimpspec/plot/mpl/magnitude.py @@ -18,11 +18,13 @@ # the LICENSES folder. from inspect import signature +from pyimpspec.circuit.circuit import Circuit from pyimpspec.data import DataSet from pyimpspec.analysis import ( KramersKronigResult, FitResult, ) +from pyimpspec.analysis.utility import _interpolate from pyimpspec.analysis.drt import DRTResult from pyimpspec.typing.helpers import ( Dict, @@ -153,6 +155,9 @@ def plot_magnitude( data.get_impedances(num_per_decade=num_per_decade) ** (-1 if admittance else 1) ) + elif line and hasattr(data, "circuit") and isinstance(data.circuit, Circuit): + x = _interpolate(data.get_frequencies(), num_per_decade=num_per_decade) + y = abs(data.circuit.get_impedances(x) ** (-1 if admittance else 1)) else: x = data.get_frequencies() y = abs(data.get_impedances() ** (-1 if admittance else 1)) diff --git a/src/pyimpspec/plot/mpl/nyquist.py b/src/pyimpspec/plot/mpl/nyquist.py index a14978f..64ec01e 100644 --- a/src/pyimpspec/plot/mpl/nyquist.py +++ b/src/pyimpspec/plot/mpl/nyquist.py @@ -18,11 +18,13 @@ # the LICENSES folder. from inspect import signature +from pyimpspec.circuit.circuit import Circuit from pyimpspec.data import DataSet from pyimpspec.analysis import ( KramersKronigResult, FitResult, ) +from pyimpspec.analysis.utility import _interpolate from pyimpspec.analysis.drt import ( DRTResult, ) @@ -145,6 +147,13 @@ def plot_nyquist( X = data.get_impedances(num_per_decade=num_per_decade) ** ( -1 if admittance else 1 ) + elif line and hasattr(data, "circuit") and isinstance(data.circuit, Circuit): + X = data.circuit.get_impedances( + _interpolate( + data.get_frequencies(), + num_per_decade=num_per_decade, + ), + ) ** (-1 if admittance else 1) else: X = data.get_impedances() ** (-1 if admittance else 1) diff --git a/src/pyimpspec/plot/mpl/phase.py b/src/pyimpspec/plot/mpl/phase.py index bc5d778..e9c4aa2 100644 --- a/src/pyimpspec/plot/mpl/phase.py +++ b/src/pyimpspec/plot/mpl/phase.py @@ -18,11 +18,13 @@ # the LICENSES folder. from inspect import signature +from pyimpspec.circuit.circuit import Circuit from pyimpspec.data import DataSet from pyimpspec.analysis import ( KramersKronigResult, FitResult, ) +from pyimpspec.analysis.utility import _interpolate from pyimpspec.analysis.drt import DRTResult from numpy import angle from pyimpspec.typing import ( @@ -161,6 +163,13 @@ def plot_phase( ** (-1 if admittance else 1), deg=True, ) + elif line and hasattr(data, "circuit") and isinstance(data.circuit, Circuit): + x = _interpolate(data.get_frequencies(), num_per_decade=num_per_decade) + y = angle( + data.circuit.get_impedances(x) + ** (-1 if admittance else 1), + deg=True, + ) else: x = data.get_frequencies() y = angle(data.get_impedances() ** (-1 if admittance else 1), deg=True) diff --git a/src/pyimpspec/plot/mpl/real.py b/src/pyimpspec/plot/mpl/real.py index 6471ef8..777e098 100644 --- a/src/pyimpspec/plot/mpl/real.py +++ b/src/pyimpspec/plot/mpl/real.py @@ -18,11 +18,13 @@ # the LICENSES folder. from inspect import signature +from pyimpspec.circuit.circuit import Circuit from pyimpspec.data import DataSet from pyimpspec.analysis import ( KramersKronigResult, FitResult, ) +from pyimpspec.analysis.utility import _interpolate from pyimpspec.analysis.drt import DRTResult from pyimpspec.typing import ( Frequencies, @@ -160,6 +162,12 @@ def plot_real( data.get_impedances(num_per_decade=num_per_decade) ** (-1 if admittance else 1) ).real + elif line and hasattr(data, "circuit") and isinstance(data.circuit, Circuit): + x = _interpolate(data.get_frequencies(), num_per_decade=num_per_decade) + y = ( + data.circuit.get_impedances(x) + ** (-1 if admittance else 1) + ).real else: x = data.get_frequencies() y = (data.get_impedances() ** (-1 if admittance else 1)).real diff --git a/tests/test_cli.py b/tests/test_cli.py index 61a0cce..2e3d70a 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -167,6 +167,7 @@ def test_output(self): set(), msg=set(self.parsers.keys()) - set(skip), ) + command: str parser: Callable for command, parser in self.parsers.items(): @@ -237,6 +238,7 @@ def test_plot(self): set(), msg=set(self.parsers.keys()) - set(skip), ) + command: str parser: Callable for command, parser in self.parsers.items(): @@ -370,6 +372,7 @@ def test_args(self): parser("-pt", "imaginary"), ] ) + for args in COMMAND_LINE_ARGUMENTS["plot"]: lines: List[str] = [] plot_command( @@ -392,6 +395,7 @@ def test_args(self): parser: Callable = lambda *a: self.parser.parse_args( ["test", *self.data_paths, *a] ) + long: str short: str dest: str @@ -479,6 +483,7 @@ def test_args(self): parser: Callable = lambda *a: self.parser.parse_args( ["zhit", *self.data_paths, *a] ) + long: str short: str dest: str @@ -501,6 +506,7 @@ def test_args(self): else: args_long = parser(long, str(value)) self.assertEqual(getattr(args_long, dest), value) + if short != "": args_short: Namespace if isinstance(value, bool): @@ -508,7 +514,9 @@ def test_args(self): else: args_short = parser(short, str(value)) self.assertEqual(getattr(args_long, dest), getattr(args_short, dest)) + COMMAND_LINE_ARGUMENTS["zhit"].append(args_long) + args: Namespace for args in COMMAND_LINE_ARGUMENTS["zhit"]: lines: List[str] = [] @@ -533,6 +541,7 @@ def test_args(self): parser: Callable = lambda *a: self.parser.parse_args( ["drt", *self.data_paths, *a] ) + long: str short: str dest: str @@ -557,21 +566,33 @@ def test_args(self): ("--threshold", "-t", "peak_threshold", 0.2), ("--max-nfev", "", "max_nfev", 4), ("--num-procs", "", "num_procs", 2), + ("--plot-frequency", "-pF", "plot_frequency", True), + ("--analyze-peaks", "-ap", "analyze_peaks", True), + ("--num-peaks", "-np", "num_peaks", 2), + ("--peak-positions", "-pp", "peak_positions", [1e-1, 0.2]), + ("--disallow-skew", "-ds", "disallow_skew", True), ]: args_long: Namespace if isinstance(value, bool): args_long = parser(long) + elif isinstance(value, list): + args_long = parser(long, *map(str, value)) else: args_long = parser(long, str(value)) self.assertEqual(getattr(args_long, dest), value) + if short != "": args_short: Namespace if isinstance(value, bool): args_short = parser(short) + elif isinstance(value, list): + args_short = parser(short, *map(str, value)) else: args_short = parser(short, str(value)) self.assertEqual(getattr(args_long, dest), getattr(args_short, dest)) + COMMAND_LINE_ARGUMENTS["drt"].append(args_long) + args: Namespace for args in COMMAND_LINE_ARGUMENTS["drt"]: lines: List[str] = [] @@ -596,6 +617,7 @@ def test_args(self): parser: Callable = lambda *a: self.parser.parse_args( ["fit", "R(RC)(RW)", *self.data_paths, *a] ) + long: str short: str dest: str @@ -612,6 +634,7 @@ def test_args(self): else: args_long = parser(long, str(value)) self.assertEqual(getattr(args_long, dest), value) + if short != "": args_short: Namespace if isinstance(value, bool): @@ -619,7 +642,9 @@ def test_args(self): else: args_short = parser(short, str(value)) self.assertEqual(getattr(args_long, dest), getattr(args_short, dest)) + COMMAND_LINE_ARGUMENTS["fit"].append(args_long) + args: Namespace for args in COMMAND_LINE_ARGUMENTS["fit"]: lines: List[str] = [] @@ -638,6 +663,7 @@ def setUpClass(cls): def test_args(self): parser: Callable = lambda *a: self.parser.parse_args(["circuit", *a]) + long: str short: str dest: str @@ -666,6 +692,7 @@ def test_args(self): else: args_long = parser(long, str(value)) self.assertEqual(getattr(args_long, dest), value) + if short != "": args_short: Namespace if isinstance(value, bool): @@ -673,7 +700,9 @@ def test_args(self): else: args_short = parser(short, str(value)) self.assertEqual(getattr(args_long, dest), getattr(args_short, dest)) + COMMAND_LINE_ARGUMENTS["circuit"].append(args_long) + args: Namespace for args in COMMAND_LINE_ARGUMENTS["circuit"]: lines: List[str] = [] @@ -692,6 +721,7 @@ def setUpClass(cls): def test_args(self): parser: Callable = lambda *a: self.parser.parse_args(["config", *a]) + long: str short: str dest: str @@ -706,6 +736,7 @@ def test_args(self): else: args_long = parser(long, str(value)) self.assertEqual(getattr(args_long, dest), value) + if short != "": args_short: Namespace if isinstance(value, bool): @@ -713,7 +744,9 @@ def test_args(self): else: args_short = parser(short, str(value)) self.assertEqual(getattr(args_long, dest), getattr(args_short, dest)) + COMMAND_LINE_ARGUMENTS["config"].append(args_long) + args: Namespace for args in COMMAND_LINE_ARGUMENTS["config"]: lines: List[str] = [] diff --git a/tests/test_drt.py b/tests/test_drt.py index 0369cbc..81690ba 100644 --- a/tests/test_drt.py +++ b/tests/test_drt.py @@ -20,29 +20,42 @@ from typing import ( Callable, Dict, + List, Optional, + Tuple, + Union, ) from unittest import TestCase from numpy import ( allclose, angle, + concatenate, + exp, + float64, + log10 as log, ndarray, + sort, ) from numpy.random import ( seed, ) +from scipy.integrate import quad from pandas import DataFrame from pyimpspec import ( Circuit, + DRTPeak, + DRTPeaks, DRTResult, DataSet, calculate_drt, fit_circuit, + generate_mock_data, parse_cdc, parse_data, ) from pyimpspec.analysis.drt import ( BHTResult, + LMResult, MRQFitResult, TRNNLSResult, TRRBFResult, @@ -58,7 +71,9 @@ from pyimpspec.typing import ( ComplexImpedances, Frequencies, + Gamma, Gammas, + TimeConstant, TimeConstants, ) from test_matplotlib import ( @@ -86,6 +101,141 @@ def progress_callback(*args, **kwargs): PROGRESS.register(progress_callback) +class TestLoewnerMethod(TestCase): + @classmethod + def setUpClass(cls): + cls.data: DataSet = generate_mock_data( + "R{R=140}L{L=1e-4}(R{R=230}C{C=1e-6})(R{R=576}C{C=1e-4})(R{R=150}L{L=4e1})", + )[0] + cls.expected_time_constants = list(map(TimeConstant, ( + 3.6e-8, + 4.2e-8, + 230.0 * 1e-6, + 576.0 * 1e-4, + 4e1 / 150.0, + ))) + cls.expected_gammas_RC = list(map(Gamma, ( + 1445.9, + 230.0, + 576.0, + ))) + cls.expected_gammas_RL = list(map(Gamma, ( + 1155.9, + 150.0, + ))) + cls.result: LMResult = calculate_drt(cls.data, method="lm") + + def calculate_drt(self, *args, **kwargs): + return calculate_drt( + self.data, + method="lm", + *args, + **kwargs, + ) + + def test_calculate_drt_lm(self): + with self.assertRaises(TypeError): + self.calculate_drt(model_order=1.0) + + with self.assertRaises(TypeError): + self.calculate_drt(model_order_method=True) + + with self.assertRaises(TypeError): + self.calculate_drt(model_order_method=1) + + with self.assertRaises(TypeError): + self.calculate_drt(model_order_method=0.5) + + def test_result_get_peaks(self): + peaks: Tuple[TimeConstants, Gammas, TimeConstants, Gammas] + peaks = self.result.get_peaks() + self.assertIsInstance(peaks, tuple) + + time_constants_RC, gammas_RC, time_constants_RL, gammas_RL = peaks + self.assertEqual(time_constants_RC.shape, gammas_RC.shape) + self.assertEqual(time_constants_RL.shape, gammas_RL.shape) + self.assertEqual(time_constants_RC.dtype, TimeConstant) + self.assertEqual(time_constants_RL.dtype, TimeConstant) + self.assertEqual(gammas_RC.dtype, Gamma) + self.assertEqual(gammas_RL.dtype, Gamma) + self.assertTrue((gammas_RC >= 0.0).all()) + self.assertTrue((gammas_RL >= 0.0).all()) + + def test_result_to_statistics_dataframe(self): + df: DataFrame = self.result.to_statistics_dataframe() + self.assertIsInstance(df, DataFrame) + self.assertEqual(len(df.keys()), 2) + self.assertTrue("Label" in df.keys()) + self.assertTrue("Value" in df.keys()) + + md: str = df.to_markdown(index=False) + self.assertTrue("Log pseudo chi-squared" in md) + self.assertTrue("Model order (k)" in md) + + def test_result_to_peaks_dataframe(self): + df: DataFrame = self.result.to_peaks_dataframe() + self.assertIsInstance(df, DataFrame) + self.assertEqual(len(df.keys()), 4, msg=df.keys()) + self.assertTrue("tau, RC (s)" in df.keys()) + self.assertTrue("gamma, RC (ohm)" in df.keys()) + self.assertTrue("tau, RL (s)" in df.keys()) + self.assertTrue("gamma, RL (ohm)" in df.keys()) + + md: str = df.to_markdown(index=False) + lines: List[str] = [line for line in md.split("\n") if line.strip() != ""] + self.assertEqual(len(lines), 5) + + def test_result_get_drt_data(self): + time_constants_RC: TimeConstants + gamma_RC: Gammas + time_constants_RL: TimeConstants + gamma_RL: Gammas + time_constants_RC, gammas_RC, time_constants_RL, gammas_RL = self.result.get_drt_data() + + for time_constants, gammas in ( + (time_constants_RC, gammas_RC), + (time_constants_RL, gammas_RL), + ): + self.assertIsInstance(time_constants, ndarray) + self.assertEqual(time_constants.dtype, TimeConstant) + self.assertIsInstance(gammas, ndarray) + self.assertEqual(gammas.dtype, Gamma) + self.assertEqual(time_constants.shape, gammas.shape) + + def test_result_get_time_constants(self): + time_constants: TimeConstants = self.result.get_time_constants() + self.assertEqual(len(time_constants), 5) + self.assertEqual(time_constants.dtype, TimeConstant) + + i: int + t: TimeConstant + for i, t in enumerate(time_constants): + self.assertAlmostEqual( + t, + self.expected_time_constants[i], + msg=f"{i=}: {t=} != {self.expected_time_constants[i]}", + ) + + def test_result_get_gammas(self): + gammas: Tuple[Gammas, Gammas] = self.result.get_gammas() + self.assertIsInstance(gammas, tuple) + self.assertEqual(len(gammas), 2) + + gammas_RC, gammas_RL = gammas + self.assertEqual(gammas_RC.shape, (3,)) + self.assertEqual(gammas_RL.shape, (2,)) + self.assertEqual(gammas_RC.dtype, Gamma) + self.assertEqual(gammas_RL.dtype, Gamma) + + i: int + g: float64 + for i, g in enumerate(gammas_RC): + self.assertAlmostEqual(g, self.expected_gammas_RC[i], places=1) + + for i, g in enumerate(gammas_RL): + self.assertAlmostEqual(g, self.expected_gammas_RL[i], places=1) + + seed(28041948) # GNU STP DATA: DataSet = parse_data("data.dta")[0] @@ -139,34 +289,34 @@ def test_result_label(self): if method == "mrq-fit": self.assertEqual(drt.get_label(), "R-2(RQ)") else: - self.assertTrue(method.upper() in drt.get_label()) + self.assertTrue(method.upper() in drt.get_label(), msg=method) def test_result_get_frequencies(self): method: str drt: DRTResult for method, drt in self._method_results.items(): - self.assertEqual(len(drt.get_frequencies()), 29) + self.assertEqual(len(drt.get_frequencies()), 29, msg=method) f_mod: Frequencies = drt.get_frequencies() - self.assertIsInstance(f_mod, ndarray) - self.assertEqual(self._f_exp.shape, f_mod.shape) - self.assertTrue(allclose(self._f_exp, f_mod)) + self.assertIsInstance(f_mod, ndarray, msg=method) + self.assertEqual(self._f_exp.shape, f_mod.shape, msg=method) + self.assertTrue(allclose(self._f_exp, f_mod), msg=method) def test_result_get_impedances(self): method: str drt: DRTResult for method, drt in self._method_results.items(): - self.assertEqual(len(drt.get_impedances()), 29) + self.assertEqual(len(drt.get_impedances()), 29, msg=method) Z_mod: ComplexImpedances = drt.get_impedances() - self.assertIsInstance(Z_mod, ndarray) - self.assertEqual(self._Z_exp.shape, Z_mod.shape) + self.assertIsInstance(Z_mod, ndarray, msg=method) + self.assertEqual(self._Z_exp.shape, Z_mod.shape, msg=method) def test_result_get_nyquist_data(self): method: str drt: DRTResult for method, drt in self._method_results.items(): Z_mod: ComplexImpedances = drt.get_impedances() - self.assertTrue(allclose(Z_mod.real, drt.get_nyquist_data()[0])) - self.assertTrue(allclose(-Z_mod.imag, drt.get_nyquist_data()[1])) + self.assertTrue(allclose(Z_mod.real, drt.get_nyquist_data()[0]), msg=method) + self.assertTrue(allclose(-Z_mod.imag, drt.get_nyquist_data()[1]), msg=method) def test_result_get_bode_data(self): method: str @@ -174,9 +324,9 @@ def test_result_get_bode_data(self): for method, drt in self._method_results.items(): f_mod: Frequencies = drt.get_frequencies() Z_mod: ComplexImpedances = drt.get_impedances() - self.assertTrue(allclose(f_mod, drt.get_bode_data()[0])) - self.assertTrue(allclose(abs(Z_mod), drt.get_bode_data()[1])) - self.assertTrue(allclose(-angle(Z_mod, deg=True), drt.get_bode_data()[2])) + self.assertTrue(allclose(f_mod, drt.get_bode_data()[0]), msg=method) + self.assertTrue(allclose(abs(Z_mod), drt.get_bode_data()[1]), msg=method) + self.assertTrue(allclose(-angle(Z_mod, deg=True), drt.get_bode_data()[2]), msg=method) def test_result_get_time_constants(self): method: str @@ -195,9 +345,15 @@ def test_result_get_gammas(self): real_gammas, imag_gammas = drt.get_gammas() self.assertIsInstance(real_gammas, ndarray) self.assertIsInstance(imag_gammas, ndarray) + elif isinstance(drt, LMResult): + gammas_RC: Gammas + gammas_RL: Gammas + gammas_RC, gammas_RL = drt.get_gammas() + self.assertIsInstance(gammas_RC, ndarray) + self.assertIsInstance(gammas_RL, ndarray) else: gamma: Gammas = drt.get_gammas() - self.assertIsInstance(gamma, ndarray) + self.assertIsInstance(gamma, ndarray, msg=method) def test_result_get_drt_data(self): method: str @@ -222,18 +378,36 @@ def test_result_get_drt_data(self): self.assertEqual(imag_gammas.size, drt.get_gammas()[1].size) self.assertTrue(allclose(real_gammas, drt.get_gammas()[0])) self.assertTrue(allclose(imag_gammas, drt.get_gammas()[1])) + elif isinstance(drt, LMResult): + self.assertEqual(len(drt.get_drt_data()), 4) + time_constants_RC: TimeConstants + gammas_RC: Gammas + time_constants_RL: TimeConstants + gammas_RL: Gammas + time_constants_RC, gammas_RC, time_constants_RL, gammas_RL = drt.get_drt_data() + self.assertIsInstance(time_constants_RC, ndarray) + self.assertIsInstance(gammas_RC, ndarray) + self.assertIsInstance(time_constants_RL, ndarray) + self.assertIsInstance(gammas_RL, ndarray) + self.assertEqual(time_constants_RC.shape, gammas_RC.shape) + self.assertEqual(time_constants_RL.shape, gammas_RL.shape) + self.assertTrue(0 <= time_constants_RC.size <= self._f_exp.size) + self.assertTrue(0 <= time_constants_RL.size <= self._f_exp.size) + self.assertTrue(allclose(sort(concatenate((time_constants_RC, time_constants_RL))), drt.get_time_constants())) + self.assertTrue(allclose(gammas_RC, drt.get_gammas()[0])) + self.assertTrue(allclose(gammas_RL, drt.get_gammas()[1])) else: - self.assertEqual(len(drt.get_drt_data()), 2) + self.assertEqual(len(drt.get_drt_data()), 2, msg=method) gamma: Gammas time_constants, gamma = drt.get_drt_data() - self.assertIsInstance(time_constants, ndarray) - self.assertIsInstance(gamma, ndarray) - self.assertEqual(time_constants.shape, gamma.shape) - self.assertTrue(time_constants.size >= self._f_exp.size) - self.assertTrue(gamma.size >= self._f_exp.size) - self.assertTrue(allclose(time_constants, drt.get_time_constants())) - self.assertTrue(allclose(gamma, drt.get_gammas())) - self.assertEqual(gamma.size, drt.get_gammas().size) + self.assertIsInstance(time_constants, ndarray, msg=method) + self.assertIsInstance(gamma, ndarray, msg=method) + self.assertEqual(time_constants.shape, gamma.shape, msg=method) + self.assertTrue(time_constants.size >= self._f_exp.size, msg=method) + self.assertTrue(gamma.size >= self._f_exp.size, msg=method) + self.assertTrue(allclose(time_constants, drt.get_time_constants()), msg=method) + self.assertTrue(allclose(gamma, drt.get_gammas()), msg=method) + self.assertEqual(gamma.size, drt.get_gammas().size, msg=method) def test_result_scores(self): method: str @@ -300,31 +474,43 @@ def test_get_peaks(self): gamma: Gammas time_constants_re: TimeConstants time_constants_im: TimeConstants - gamma_re: Gammas - gamma_im: Gammas + gammas_re: Gammas + gammas_im: Gammas + drt = self._method_results["tr-nnls"] time_constants, gamma = drt.get_peaks() self.assertGreater(time_constants.size, 0) self.assertEqual(time_constants.size, gamma.size) self.assertIsInstance(drt.to_peaks_dataframe(), DataFrame) + drt = self._method_results["tr-rbf"] time_constants, gamma = drt.get_peaks() self.assertGreater(time_constants.size, 0) self.assertEqual(time_constants.size, gamma.size) self.assertIsInstance(drt.to_peaks_dataframe(), DataFrame) + drt = self._method_results["bht"] - time_constants_re, gamma_re, time_constants_im, gamma_im = drt.get_peaks() + time_constants_re, gammas_re, time_constants_im, gammas_im = drt.get_peaks() self.assertGreater(time_constants_re.size, 0) - self.assertEqual(time_constants_re.size, gamma_re.size) + self.assertEqual(time_constants_re.size, gammas_re.size) self.assertGreater(time_constants_im.size, 0) - self.assertEqual(time_constants_im.size, gamma_im.size) + self.assertEqual(time_constants_im.size, gammas_im.size) self.assertIsInstance(drt.to_peaks_dataframe(), DataFrame) + drt = self._method_results["mrq-fit"] time_constants, gamma = drt.get_peaks() self.assertGreater(time_constants.size, 0) self.assertEqual(time_constants.size, gamma.size) self.assertIsInstance(drt.to_peaks_dataframe(), DataFrame) + drt = self._method_results["lm"] + time_constants_RC, gammas_RC, time_constants_RL, gammas_RL = drt.get_peaks() + self.assertGreater(time_constants_RC.size, 0) + self.assertEqual(time_constants_RC.size, gammas_RC.size) + self.assertEqual(time_constants_RL.size, 0) + self.assertEqual(time_constants_RL.size, gammas_RL.size) + self.assertIsInstance(drt.to_peaks_dataframe(), DataFrame) + def test_method_tr_nnls(self): method: str = "tr-nnls" self.assertRaises( @@ -339,9 +525,10 @@ def test_method_tr_nnls(self): method=method, mode="crivens", ) - drt: TRNNLSResult = self.wrapper(method="tr-nnls", mode="real") - self.wrapper(method="tr-nnls", mode="imaginary") - self.wrapper(method="tr-nnls", mode="real", lambda_value=1e-3) + + drt: TRNNLSResult = self.wrapper(method=method, mode="real") + self.wrapper(method=method, mode="imaginary") + self.wrapper(method=method, mode="real", lambda_value=1e-3) def test_method_tr_rbf(self): method: str = "tr-rbf" @@ -421,11 +608,12 @@ def test_method_tr_rbf(self): self.wrapper(method=method, rbf_shape=rbf_shape) def test_credible_intervals(self): - # TODO: Switch to an impedance spectrum with some noise? + # TODO: Switch to an impedance spectrum with some noise? # Calculating credible intervals seems to be faster and thus # the timeout would not need to be raised. Other test cases might # need to be updated in terms of expected values. drt: DRTResult = self.wrapper(method="tr-rbf", credible_intervals=True, timeout=300) + time_constants: TimeConstants mean_gamma: Gammas lower_bound: Gammas @@ -436,6 +624,7 @@ def test_credible_intervals(self): lower_bound, upper_bound, ) = drt.get_drt_credible_intervals_data() + self.assertIsInstance(time_constants, ndarray) self.assertIsInstance(mean_gamma, ndarray) self.assertIsInstance(lower_bound, ndarray) @@ -693,6 +882,97 @@ def test_method_mrq_fit(self): circuit=parse_cdc(self._circuit.to_string()), ) + def test_method_lm(self): + method: str = "lm" + self.assertRaises( + TypeError, + self.wrapper, + method=method, + model_order=0.2, + ) + self.assertRaises( + TypeError, + self.wrapper, + method=method, + model_order="crivens", + ) + self.assertRaises( + TypeError, + self.wrapper, + method=method, + model_order_method=2, + ) + self.assertRaises( + ValueError, + self.wrapper, + method=method, + model_order_method="crivens", + ) + + drt: LMResult = self.wrapper(method=method) + self.wrapper(method=method, model_order=2) + self.wrapper(method=method, model_order_method="pseudo_chisqr") + + def test_analyze_peaks(self): + result: DRTResult + for result in self._method_results.values(): + peaks: Union[DRTPeaks, Tuple[DRTPeaks, ...]] = result.analyze_peaks() + if isinstance(result, BHTResult): + self.assertIsInstance(peaks, tuple) + self.assertEqual(len(peaks), 2) + self.assertIsInstance(peaks[0], DRTPeaks) + self.assertIsInstance(peaks[1], DRTPeaks) + elif isinstance(result, LMResult): + self.assertIsInstance(peaks, tuple) + self.assertEqual(len(peaks), 2) + self.assertIsInstance(peaks[0], DRTPeaks) + self.assertIsInstance(peaks[1], DRTPeaks) + else: + self.assertIsInstance(peaks, DRTPeaks) + + data: DataSet = generate_mock_data( + "R{R=140}(R{R=230}C{C=1e-6})(R{R=576}C{C=1e-4})", + )[0] + expected_values: List[float] = [230.0, 576.0] + result = calculate_drt(data, method="tr-rbf") + + with self.assertRaises(TypeError): + result.analyze_peaks(num_peaks=2.3) + + with self.assertRaises(TypeError): + result.analyze_peaks(peak_positions=(4, 2)) + + with self.assertRaises(TypeError): + result.analyze_peaks(disallow_skew=2) + + method: str + for method in _METHODS: + if method in ("bht", "lm", "mrq-fit"): + continue + + result = calculate_drt(data, method=method) + peaks = result.analyze_peaks(num_peaks=2) + self.assertIsInstance(peaks, DRTPeaks) + self.assertEqual(peaks.get_num_peaks(), 2) + + time_constants: TimeConstants = peaks.get_time_constants() + + i: int + peak: DRTPeak + for i, peak in enumerate(peaks): + area: float = quad( + func=lambda x: peak.get_gammas(10**x), + a=log(min(time_constants)), + b=log(max(time_constants)), + )[0] / log(exp(1)) + + self.assertAlmostEqual(peaks.get_peak_area(i), area) + self.assertAlmostEqual( + area, + expected_values[i], + delta=0.1*expected_values[i], + ) + def test_matplotlib(self): drt: DRTResult for drt in self._method_results.values(): diff --git a/tests/test_element.py b/tests/test_element.py index 74fb055..7e5c559 100644 --- a/tests/test_element.py +++ b/tests/test_element.py @@ -40,21 +40,25 @@ inf, ) from pyimpspec import ( - Capacitor, Connection, - ConstantPhaseElement, Container, - DeLevieFiniteLength, Element, ElementDefinition, + ParameterDefinition, + get_elements, + register_element, +) +from pyimpspec.circuit.elements import ( + Capacitor, + ConstantPhaseElement, + DeLevieFiniteLength, Gerischer, HavriliakNegami, HavriliakNegamiAlternative, Inductor, - KramersKronigRC, KramersKronigAdmittanceRC, + KramersKronigRC, ModifiedInductor, - ParameterDefinition, Resistor, TransmissionLineModel, TransmissionLineModelBlockingCPE, @@ -67,8 +71,6 @@ WarburgOpen, WarburgShort, ZARC, - get_elements, - register_element, ) from pyimpspec.circuit.registry import ( _validate_impedances, @@ -473,7 +475,7 @@ def check( with self.assertRaises(TypeError): get_elements(default_only="test") - elements: Dict[str, Type[Element]] = get_elements() + elements: Dict[str, Type[Element]] = get_elements(private=True) symbols: List[str] = list(elements.keys()) check(symbols, elements, "C", Capacitor) check(symbols, elements, "G", Gerischer) diff --git a/tests/test_fitting.py b/tests/test_fitting.py index 0bfde9f..2b8e46c 100644 --- a/tests/test_fitting.py +++ b/tests/test_fitting.py @@ -34,9 +34,13 @@ from pyimpspec import ( Circuit, DataSet, + Element, + FitIdentifiers, FitResult, FittedParameter, fit_circuit, + generate_mock_data, + generate_fit_identifiers, parse_cdc, parse_data, ) @@ -51,6 +55,7 @@ ) from pyimpspec.typing.helpers import ( Callable, + Dict, List, Optional, Tuple, @@ -431,3 +436,119 @@ def test_matplotlib(self): check_mpl_return_values(self, *mpl.plot_fit(self.result, data=DATA)) with self.assertRaises(AttributeError): mpl.plot_fit(DATA, data=DATA) + + def test_fit_identifiers(self): + data: DataSet = generate_mock_data("CIRCUIT_5", noise=5e-2, seed=42)[0] + + circuit: Circuit = parse_cdc("R(RQ)(RQ)(RQ)") + elements: List[Element] = circuit.get_elements() + R1, R2, Q1, R3, Q2, R4, Q3 = elements + identifiers: Dict[Element, FitIdentifiers] = generate_fit_identifiers(circuit) + + # Make sure the types are as expected + self.assertIsInstance(identifiers, dict) + self.assertTrue(all(map(lambda k: isinstance(k, Element), identifiers.keys()))) + self.assertTrue(all(map(lambda v: isinstance(v, FitIdentifiers), identifiers.values()))) + self.assertTrue(all(map(lambda v: all(map(lambda kv: isinstance(kv[0], str) and isinstance(kv[1], str), v.items())), identifiers.values()))) + + # Make sure that each element was processed + self.assertTrue(all(map(lambda e: e in identifiers, elements))) + + # Make sure that all parameter symbols are valid + for element in elements: + for symbol in identifiers[element]: + element.get_value(symbol) + + # Make sure all identifiers are unique + self.assertEqual(len(set([identifiers[R].R for R in (R1, R2, R3, R4)])), 4) + self.assertEqual(len(set([identifiers[Q].Y for Q in (Q1, Q2, Q3)])), 3) + self.assertEqual(len(set([identifiers[Q].n for Q in (Q1, Q2, Q3)])), 3) + + # Make sure the identifiers have the expected values + self.assertEqual(identifiers[R1].R, "R_0") + + self.assertEqual(identifiers[R2].R, "R_1") + self.assertEqual(identifiers[Q1].Y, "Y_2") + self.assertEqual(identifiers[Q1].n, "n_2") + + self.assertEqual(identifiers[R3].R, "R_3") + self.assertEqual(identifiers[Q2].Y, "Y_4") + self.assertEqual(identifiers[Q2].n, "n_4") + + self.assertEqual(identifiers[R4].R, "R_5") + self.assertEqual(identifiers[Q3].Y, "Y_6") + self.assertEqual(identifiers[Q3].n, "n_6") + + # Make sure the identifiers can be accessed + self.assertEqual(identifiers[R1].R, identifiers[R1]["R"]) + self.assertEqual(identifiers[Q1].Y, identifiers[Q1]["Y"]) + self.assertEqual(identifiers[Q1].n, identifiers[Q1]["n"]) + + with self.assertRaises(AttributeError): + identifiers[R1].Y + + with self.assertRaises(AttributeError): + identifiers[R1]["Y"] + + with self.assertRaises(AttributeError): + identifiers[Q1].N + + with self.assertRaises(AttributeError): + identifiers[Q1]["N"] + + with self.assertRaises(AttributeError): + identifiers[Q1].R + + with self.assertRaises(AttributeError): + identifiers[Q1]["R"] + + fit: FitResult = fit_circuit( + circuit, + data, + method="least_squares", + weight="boukamp", + constraint_expressions={ + identifiers[R3].R: f"{identifiers[R2].R} + alpha", + identifiers[R4].R: f"{identifiers[R3].R} - beta", + identifiers[Q2].Y: f"{identifiers[Q1].Y} + gamma", + identifiers[Q3].Y: f"{identifiers[Q2].Y} + delta", + }, + constraint_variables=dict( + alpha=dict( + value=500, + min=0, + ), + beta=dict( + value=300, + min=0, + ), + gamma=dict( + value=1e-8, + min=0, + ), + delta=dict( + value=2e-7, + min=0, + ), + ), + ) + + R1, R2, Q1, R3, Q2, R4, Q3 = fit.circuit.get_elements() + + # Make sure the values are in the expected order + self.assertLess(R2.get_value("R"), R4.get_value("R")) + self.assertLess(R4.get_value("R"), R3.get_value("R")) + self.assertLess(Q1.get_value("Y"), Q2.get_value("Y")) + self.assertLess(Q2.get_value("Y"), Q3.get_value("Y")) + + # Make sure the values are as expected + self.assertAlmostEqual(R1.get_value("R"), 143, delta=1) + self.assertAlmostEqual(R2.get_value("R"), 229, delta=1) + self.assertAlmostEqual(R3.get_value("R"), 854, delta=1) + self.assertAlmostEqual(R4.get_value("R"), 475, delta=1) + self.assertAlmostEqual(Q1.get_value("Y"), 1.78e-7, delta=5e-9) + self.assertAlmostEqual(Q2.get_value("Y"), 7.83e-7, delta=5e-9) + self.assertAlmostEqual(Q3.get_value("Y"), 1.61e-5, delta=5e-7) + self.assertAlmostEqual(Q1.get_value("n"), 0.913, delta=5e-2) + self.assertAlmostEqual(Q2.get_value("n"), 0.856, delta=5e-2) + self.assertAlmostEqual(Q3.get_value("n"), 0.952, delta=5e-2)