From da306fafa88b001850eaae42e70bfb6f9ad4837f Mon Sep 17 00:00:00 2001 From: Andreas Copan Date: Thu, 13 Feb 2025 13:05:35 -0500 Subject: [PATCH] New: Simplifies plotting --- autochem/rate/_01const.py | 89 +++++++++++++++++++++++++++++++++--- autochem/rate/_02rate.py | 36 ++++++++++++++- autochem/tests/test__rate.py | 16 +++---- autochem/unit_/__init__.py | 3 +- autochem/unit_/_unit.py | 6 ++- autochem/util/__init__.py | 4 +- autochem/util/plot.py | 66 -------------------------- 7 files changed, 134 insertions(+), 86 deletions(-) delete mode 100644 autochem/util/plot.py diff --git a/autochem/rate/_01const.py b/autochem/rate/_01const.py index 78494399..00eb5268 100644 --- a/autochem/rate/_01const.py +++ b/autochem/rate/_01const.py @@ -4,8 +4,10 @@ from collections.abc import Mapping from typing import Annotated, ClassVar +import altair import more_itertools as mit import numpy +import pandas import pint import pydantic import xarray @@ -16,7 +18,7 @@ from .. import unit_ from ..unit_ import UNITS, Dimension, UnitManager, Units, UnitsData from ..util import chemkin -from ..util.type_ import Frozen, NDArray_, Scalable, Scalers, SubclassTyped +from ..util.type_ import Frozen, NDArray_, Number, Scalable, Scalers, SubclassTyped from ._00func import ( BlendingFunction_, extract_blending_function_from_chemkin_parse_results, @@ -79,6 +81,79 @@ def process_input( return t, p + def process_output(self, ktp: ArrayLike) -> NDArray[numpy.float64]: + """Normalize rate constant output, clipping unphyiscal negative values. + + :param ktp: Rate constant values + :return: Rate constant values + """ + return numpy.where(numpy.less_equal(ktp, 0), numpy.nan, ktp) + + def display( + self, + others: "Mapping[str, RateConstant] | None" = None, + t_range: tuple[Number, Number] = (400, 1250), + p: Number = 1, + units: UnitsData | None = None, + label: str = "This work", + x_label: str = "1000/T", + y_label: str = "k", + ) -> altair.Chart: + """Display as an Arrhenius plot. + + :param others: Other rate constants by label + :param t_range: Temperature range + :param p: Pressure + :param units: Units + :param x_label: X-axis label + :param y_label: Y-axis label + :return: Chart + """ + # Process units + units = UNITS if units is None else Units.model_validate(units) + x_unit = unit_.pretty_string(units.temperature**-1) + y_unit = unit_.pretty_string(units.rate_constant(self.order)) + + # Add units to labels + x_label = f"{x_label} ({x_unit})" + y_label = f"{y_label} ({y_unit})" + + # Gather functions + others = {} if others is None else others + funcs = {label: self, **others} + + # Gether data from functons + t = numpy.linspace(*t_range, num=500) + data_dct = {lab: func(t, p) for lab, func in funcs.items()} + data = pandas.DataFrame({"x": numpy.divide(1000, t), **data_dct}) + + # Determine exponent range + vals_arr = numpy.array(list(data_dct.values())) + is_nan = numpy.isnan(vals_arr) + exp_arr = numpy.log10(vals_arr, where=~is_nan) + exp_arr[is_nan] = 0.0 + exp_arr = numpy.rint(exp_arr).astype(int) + exp_max = numpy.max(exp_arr).item() + exp_min = numpy.min(exp_arr).item() + y_vals = [10**x for x in range(exp_min, exp_max + 2)] + + # Prepare encoding parameters + x = altair.X("x", title=x_label) + y = ( + altair.Y("value:Q", title=y_label) + .scale(type="log") + .axis(format=".1e", values=y_vals) + ) + color = "key:N" if others else altair.Undefined + + # Create chart + return ( + altair.Chart(data) + .mark_line() + .transform_fold(fold=list(data_dct.keys())) + .encode(x=x, y=y, color=color) + ) + class RawRateConstant(RateConstant): ts: list[float] @@ -111,7 +186,7 @@ def __call__( """Evaluate rate constant.""" t, p = self.process_input(t, p) ktp: NDArray[numpy.float64] = self.ktp.sel(t=t, p=p, method="ffill").data - return ktp + return self.process_output(ktp) class ParamRateConstant(RateConstant): @@ -160,7 +235,7 @@ def __call__( t, p = self.process_input(t, p) r = unit_.system.gas_constant_value(UNITS) ktp = self.A * (t**self.b) * numpy.exp(-self.E / (r * t)) - return ktp + return self.process_output(ktp) class BlendedRateConstant(ParamRateConstant, abc.ABC): # type: ignore[misc] @@ -245,7 +320,7 @@ def __call__( _, k_high = self.arrhenius_functions p_r = self.effective_reduced_pressure(t, p) ktp = k_high(t) * p_r / (1 + p_r) * self.function(t, p_r) - return ktp + return self.process_output(ktp) class ActivatedRateConstant(BlendedRateConstant): @@ -264,7 +339,7 @@ def __call__( k_low, _ = self.arrhenius_functions p_r = self.effective_reduced_pressure(t, p) ktp = k_low(t) / (1 + p_r) * self.function(t, p_r) - return ktp + return self.process_output(ktp) class PlogRateConstant(ParamRateConstant): @@ -305,7 +380,7 @@ def __call__( elif p == p0: ktp = kt0 - return ktp + return self.process_output(ktp) @property def arrhenius_functions(self) -> list[ArrheniusRateConstant]: @@ -406,7 +481,7 @@ def __call__( # AVC: I don't understand why I need to transpose the coefficient matrix here ktp = chebyshev.chebgrid2d(t_, p_, self.coeffs.T) - return ktp + return self.process_output(ktp) RateConstant_ = Annotated[ diff --git a/autochem/rate/_02rate.py b/autochem/rate/_02rate.py index 6776b12e..35bda011 100644 --- a/autochem/rate/_02rate.py +++ b/autochem/rate/_02rate.py @@ -5,6 +5,7 @@ from collections.abc import Mapping, Sequence from typing import ClassVar +import altair import numpy import pint import pydantic @@ -12,7 +13,7 @@ from ..unit_ import UnitsData from ..util import chemkin -from ..util.type_ import Scalable, Scalers +from ..util.type_ import Number, Scalable, Scalers from ._01const import ( ArrheniusRateConstant, ParamRateConstant, @@ -77,6 +78,39 @@ def __call__( """ return self.rate_constant(t=t, p=p, units=units) + def display( + self, + others: "Mapping[str, Rate] | None" = None, + t_range: tuple[Number, Number] = (400, 1250), + p: Number = 1, + units: UnitsData | None = None, + label: str = "This work", + x_label: str = "1000/T", + y_label: str = "k", + ) -> altair.Chart: + """Display as an Arrhenius plot. + + :param others: Other rate constants by label + :param t_range: Temperature range + :param p: Pressure + :param units: Units + :param x_label: X-axis label + :param y_label: Y-axis label + :return: Chart + """ + if others is not None: + others = {lab: obj.rate_constant for lab, obj in others.items()} + + return self.rate_constant.display( + others=others, + t_range=t_range, + p=p, + units=units, + label=label, + x_label=x_label, + y_label=y_label, + ) + # Constructors def from_chemkin_string( diff --git a/autochem/tests/test__rate.py b/autochem/tests/test__rate.py index 694558af..8a91484c 100644 --- a/autochem/tests/test__rate.py +++ b/autochem/tests/test__rate.py @@ -1,10 +1,8 @@ """Test autochem.rate.""" -import numpy import pytest from autochem import rate -from autochem.util import plot # Elementary SIMPLE = { @@ -115,19 +113,21 @@ def test__from_chemkin_string(name, data, check_roundtrip: bool): # Read k = rate.from_chemkin_string(chem_str, units=units) - # Plot - T = numpy.linspace(400, 1250, 500) - P = 10 - k(T, P) - plot.arrhenius_plot(T, k(T, P), y_unit=k.unit) - # Write chem_str_ = rate.chemkin_string(k) print(chem_str_) k_ = rate.from_chemkin_string(chem_str_) + + # Check roundtrip if check_roundtrip: assert k == k_, f"\n {k}\n!= {k_}" + # Plot against another rate + other_units = SIMPLE.get("units") + other_chem_str = SIMPLE.get("chemkin") + other_k = rate.from_chemkin_string(other_chem_str, units=other_units) + k.display({"Other": other_k}) + if __name__ == "__main__": # test__from_chemkin_string("FALLOFF", FALLOFF) diff --git a/autochem/unit_/__init__.py b/autochem/unit_/__init__.py index 8970c2f6..2660e7c4 100644 --- a/autochem/unit_/__init__.py +++ b/autochem/unit_/__init__.py @@ -2,7 +2,7 @@ from . import system from ._manager import UnitManager, manage_units -from ._unit import string +from ._unit import pretty_string, string from .system import UNITS, Dimension, Units, UnitsData __all__ = [ @@ -12,6 +12,7 @@ "Units", "UnitsData", "manage_units", + "pretty_string", "string", "system", ] diff --git a/autochem/unit_/_unit.py b/autochem/unit_/_unit.py index c0998617..eb009eef 100644 --- a/autochem/unit_/_unit.py +++ b/autochem/unit_/_unit.py @@ -4,4 +4,8 @@ def string(unit: pint.Unit) -> str: - return format(unit, "~") + return format(unit, "~C") + + +def pretty_string(unit: pint.Unit) -> str: + return format(unit, "~P") diff --git a/autochem/util/__init__.py b/autochem/util/__init__.py index e2ddf82a..e173770e 100644 --- a/autochem/util/__init__.py +++ b/autochem/util/__init__.py @@ -1,5 +1,5 @@ """Utilities.""" -from . import chemkin, plot, type_ +from . import chemkin, type_ -__all__ = ["chemkin", "plot", "type_"] +__all__ = ["chemkin", "type_"] diff --git a/autochem/util/plot.py b/autochem/util/plot.py deleted file mode 100644 index 9799dcef..00000000 --- a/autochem/util/plot.py +++ /dev/null @@ -1,66 +0,0 @@ -"""Plotting functions.""" - -from collections.abc import Mapping - -import altair -import numpy -import pandas -import pint -from numpy.typing import ArrayLike - - -def arrhenius_plot( - T: ArrayLike, # noqa: N803 - kT_: ArrayLike | Mapping[str, ArrayLike], # noqa: N803 - x_label: str = "1000/T", - x_unit: str | pint.Unit = "1/K", - y_label: str = "k", - y_unit: str | pint.Unit | None = None, -) -> altair.Chart: - """Generate an Arrhenius plot. - - :param T: Temperature array - :param kT_: Rate array(s), as a dictionary by legend label if multiple - :param x_label: X-axis label - :param x_unit: X-axis unit - :param y_label: Y-axis label - :param y_unit: Y-axis unit - :return: Chart - """ - x_unit = format(pint.Unit(x_unit), "~P") - y_unit = None if y_unit is None else format(pint.Unit(y_unit), "~P") - x_label = f"{x_label} ({x_unit})" - y_label = y_label if y_unit is None else f"{y_label} ({y_unit})" - - # Add data to DataFrame, dropping negative values for log plot - has_series_labels = isinstance(kT_, Mapping) - kTs = kT_ if isinstance(kT_, Mapping) else {"y": kT_} - kTs = {lb: numpy.where(numpy.less(a, 0), numpy.nan, a) for lb, a in kTs.items()} - data = pandas.DataFrame({"x": numpy.divide(1000, T), **kTs}) - - # Determine the exponent range - k_array = numpy.array(list(kTs.values())) - is_nan = numpy.isnan(k_array) - e_array = numpy.log10(k_array, where=~is_nan) - e_array[is_nan] = 0.0 - e_array = numpy.rint(e_array).astype(int) - e_max = numpy.max(e_array).item() - e_min = numpy.min(e_array).item() - values = [10**x for x in range(e_min, e_max + 2)] - - # Prepare encoding parameters - x = altair.X("x", title=x_label) - y = ( - altair.Y("value:Q", title=y_label) - .scale(type="log") - .axis(format=".1e", values=values) - ) - color = "key:N" if has_series_labels else altair.Undefined - - # Create chart - return ( - altair.Chart(data) - .mark_line() - .transform_fold(fold=list(kTs.keys())) - .encode(x=x, y=y, color=color) - )