polykin.distributions

MixtureDistribution

Mixture chain-length distribution.

This kind of distributions are instantiated indirectly by doing linear combinations of IndividualDistribution objects.

Examples:

>>> from polykin.distributions import Flory, SchulzZimm
>>> a = Flory(100, M0=0.050, name='A')
>>> b = SchulzZimm(100, PDI=3., M0=0.10, name='B')
>>> c = 0.3*a + 0.7*b # c is now a MixtureDistribution instance
>>> c
type: MixtureDistribution
name: A+B
DPn:  100.0
DPw:  269.7
DPz:  474.9
PDI:  3.12
M0:   0.077 kg/mol
Mn:   7.692 kg/mol
Mw:   23.985 kg/mol
Mz:   45.635 kg/mol

 #   Weight   Distribution        DPn        DPw    PDI
-------------------------------------------------------
 1    0.300          Flory   1.00e+02   1.99e+02   1.99
 2    0.700     SchulzZimm   1.00e+02   3.00e+02   3.00

>>> c.pdf(c.DPn)
0.002802983984583185

>>> c.cdf([c.DPn, c.DPw, c.DPz])
array([0.21950423, 0.61773034, 0.85164309])

Source code in src/polykin/distributions/base.py

class MixtureDistribution(Distribution):
    r"""Mixture chain-length distribution.

    This kind of distributions are instantiated _indirectly_ by doing linear
    combinations of `IndividualDistribution` objects.

    Examples
    --------
    >>> from polykin.distributions import Flory, SchulzZimm
    >>> a = Flory(100, M0=0.050, name='A')
    >>> b = SchulzZimm(100, PDI=3., M0=0.10, name='B')
    >>> c = 0.3*a + 0.7*b # c is now a MixtureDistribution instance
    >>> c
    type: MixtureDistribution
    name: A+B
    DPn:  100.0
    DPw:  269.7
    DPz:  474.9
    PDI:  3.12
    M0:   0.077 kg/mol
    Mn:   7.692 kg/mol
    Mw:   23.985 kg/mol
    Mz:   45.635 kg/mol
    <BLANKLINE>
     #   Weight   Distribution        DPn        DPw    PDI
    -------------------------------------------------------
     1    0.300          Flory   1.00e+02   1.99e+02   1.99
     2    0.700     SchulzZimm   1.00e+02   3.00e+02   3.00

    >>> c.pdf(c.DPn)
    0.002802983984583185

    >>> c.cdf([c.DPn, c.DPw, c.DPz])
    array([0.21950423, 0.61773034, 0.85164309])

    """

    def __init__(self,
                 components: dict[IndividualDistribution, float],
                 name: str = ''
                 ) -> None:

        self.__components = components
        self.name = name

    def __add__(self, other):
        if isinstance(other, MixtureDistribution):
            return MixtureDistribution(add_dicts(self.components,
                                                 other.components),
                                       name=self.name+'+'+other.name)
        elif isinstance(other, IndividualDistribution):
            return MixtureDistribution(add_dicts(self.components, {other: 1}),
                                       name=self.name+'+'+other.name)
        elif isinstance(other, (int, float)):
            return self
        else:
            return NotImplemented

    def __radd__(self, other):
        return self.__add__(other)

    def __iter__(self):
        return self.components

    def __contains__(self, component):
        return (component in self.components)

    def __repr__(self):
        if len(self.components) > 0:
            return super().__repr__() + "\n\n" + self.components_table
        else:
            return 'empty'

    def __bool__(self):
        return bool(self.components)

    @property
    def components(self) -> dict[IndividualDistribution, float]:
        r"""Individual components of the mixture distribution.

        Returns
        -------
        dict[IndividualDistribution, float]
            Dictionary of individual distributions and corresponding mass
            weight.
        """
        return self.__components

    @property
    def components_table(self) -> str:
        r"""Table of individual components of the mixture distribution.

        Returns
        -------
        str
            Table with key properties of each component.
        """
        result = 'empty'
        if self.components:
            spacer = ' '*3
            header = [f"{'#':>2}", f"{'Weight':>6}", f"{'Distribution':>12}",
                      f"{'DPn':>8}", f"{'DPw':>8}", f"{'PDI':>4}"]
            header = (spacer).join(header)
            ruler = f"{'-'*len(header)}"
            table = [header, ruler]
            for i, (d, w) in enumerate(self.components.items()):
                row = [f"{i+1:2}", f"{w:>6.3f}",
                       f"{d.__class__.__name__:>12}",
                       f"{d.DPn:>4.2e}", f"{d.DPw:>4.2e}", f"{d.PDI:>4.2f}"]
                table.append((spacer).join(row))
            result = ("\n").join(table)
        return result

    @property
    def _molefrac(self) -> np.ndarray:
        """Mole fraction of each individual distribution."""
        xn = np.empty(len(self.components))
        for i, (d, w) in enumerate(self.__iter__().items()):
            xn[i] = w/d.Mn
        xn[:] /= xn.sum()
        return xn

    def _moment_mass(self, order, shift=0):
        xn = self._molefrac
        result = 0
        for i, d in enumerate(self.__iter__()):
            result += xn[i]*d._moment_mass(order, shift)
        return result

    def _pdf(self, size, order, sizeasmass):
        xn = self._molefrac
        numerator = 0
        denominator = 0
        for i, d in enumerate(self.__iter__()):
            term1 = xn[i]*d._moment_mass(order)
            term2 = term1*d._pdf(size, order, sizeasmass)
            denominator += term1
            numerator += term2
        return numerator/denominator

    def _cdf(self, size, order, sizeasmass):
        xn = self._molefrac
        numerator = 0
        denominator = 0
        for i, d in enumerate(self.__iter__()):
            term1 = xn[i]*d._moment_mass(order)
            term2 = term1*d._cdf(size, order, sizeasmass)
            denominator += term1
            numerator += term2
        return numerator/denominator

    def _xrange_plot(self, sizeasmass):
        xrange = np.empty(2)
        xrange[0] = min([d._xrange_plot(sizeasmass)[0]
                        for d in self.__iter__()])
        xrange[1] = max([d._xrange_plot(sizeasmass)[1]
                         for d in self.__iter__()])
        return xrange

DPn property

DPn: float

Number-average degree of polymerization, \(DP_n\).

DPw property

DPw: float

Mass-average degree of polymerization, \(DP_w\).

DPz property

DPz: float

z-average degree of polymerization, \(DP_z\).

M0 property

M0: float

Number-average molar mass of the repeating units, \(M_0=M_n/DP_n\).

Mn property

Mn: float

Number-average molar mass, \(M_n\).

Mw property

Mw: float

Weight-average molar mass, \(M_w\).

Mz property

Mz: float

z-average molar mass, \(M_z\).

PDI property

PDI: float

Polydispersity index, \(M_w/M_n\).

cdf

cdf(
    size: Union[float, FloatArrayLike],
    kind: Literal["number", "mass"] = "mass",
    sizeasmass: bool = False,
) -> Union[float, FloatArray]

Evaluate the cumulative distribution function:

\[ F(s) = \frac{\sum_{k=1}^{s}k^m\,p(k)}{\sum_{k=1}^{\infty}k^m\,p(k)} \]

or

\[ F(s) = \frac{\int_{0}^{s}x^m\,p(x)\mathrm{d}x} {\int_{0}^{\infty}x^m\,p(x)\mathrm{d}x} \]

where \(m\) is the order (0: number, 1: mass).

PARAMETER	DESCRIPTION
size	Chain length or molar mass. TYPE: float | FloatArrayLike
kind	Kind of distribution. TYPE: Literal['number', 'mass'] DEFAULT: 'mass'
sizeasmass	Switch size input between chain-length (if False) or molar mass (if True). TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
float | FloatArray	Cumulative probability.

Source code in src/polykin/distributions/base.py

def cdf(self,
        size: Union[float, FloatArrayLike],
        kind: Literal['number', 'mass'] = 'mass',
        sizeasmass: bool = False,
        ) -> Union[float, FloatArray]:
    r"""Evaluate the cumulative distribution function:

    $$
    F(s) = \frac{\sum_{k=1}^{s}k^m\,p(k)}{\sum_{k=1}^{\infty}k^m\,p(k)}
    $$

    or

    $$
    F(s) = \frac{\int_{0}^{s}x^m\,p(x)\mathrm{d}x}
           {\int_{0}^{\infty}x^m\,p(x)\mathrm{d}x}
    $$

    where $m$ is the order (0: number, 1: mass).

    Parameters
    ----------
    size : float | FloatArrayLike
        Chain length or molar mass.
    kind : Literal['number', 'mass']
        Kind of distribution.
    sizeasmass : bool
        Switch size input between chain-length (if `False`) or molar
        mass (if `True`).

    Returns
    -------
    float | FloatArray
        Cumulative probability.
    """
    # Check inputs
    kind = self._verify_kind(kind)
    if kind == 'gpc':
        custom_error('kind', kind, ValueError,
                     "Please use `mass` instead.")
    order = self.kind_order[kind]
    self._verify_sizeasmass(sizeasmass)
    # Convert list to ndarray
    if isinstance(size, (list, tuple)):
        size = np.array(size)
    # Math is done by the corresponding subclass method
    return self._cdf(size, order, sizeasmass)

components property

components: dict[IndividualDistribution, float]

Individual components of the mixture distribution.

RETURNS	DESCRIPTION
dict[IndividualDistribution, float]	Dictionary of individual distributions and corresponding mass weight.

components_table property

components_table: str

Table of individual components of the mixture distribution.

RETURNS	DESCRIPTION
str	Table with key properties of each component.

pdf

pdf(
    size: Union[float, FloatArrayLike],
    kind: Literal["number", "mass", "gpc"] = "mass",
    sizeasmass: bool = False,
) -> Union[float, FloatArray]

Evaluate the probability density function, \(p(k)\).

PARAMETER	DESCRIPTION
size	Chain length or molar mass. TYPE: float | FloatArrayLike
kind	Kind of distribution. TYPE: Literal['number', 'mass', 'gpc'] DEFAULT: 'mass'
sizeasmass	Switch size input between chain-length (if False) or molar mass (if True). TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
float | FloatArray	Probability density.

Source code in src/polykin/distributions/base.py

def pdf(self,
        size: Union[float, FloatArrayLike],
        kind: Literal['number', 'mass', 'gpc'] = 'mass',
        sizeasmass: bool = False,
        ) -> Union[float, FloatArray]:
    r"""Evaluate the probability density function, $p(k)$.

    Parameters
    ----------
    size : float | FloatArrayLike
        Chain length or molar mass.
    kind : Literal['number', 'mass', 'gpc']
        Kind of distribution.
    sizeasmass : bool
        Switch size input between chain-length (if `False`) or molar
        mass (if `True`).

    Returns
    -------
    float | FloatArray
        Probability density.
    """
    # Check inputs
    self._verify_sizeasmass(sizeasmass)
    order = self.kind_order[self._verify_kind(kind)]
    # Convert list to ndarray
    if isinstance(size, (list, tuple)):
        size = np.array(size)
    # Math is done by the corresponding subclass method
    return self._pdf(size, order, sizeasmass)

plot

plot(
    kind: Union[
        Literal["number", "mass", "gpc"],
        list[Literal["number", "mass", "gpc"]],
    ] = "mass",
    sizeasmass: bool = False,
    xscale: Literal["auto", "linear", "log"] = "auto",
    xrange: Union[tuple[float, float], None] = None,
    cdf: Literal[0, 1, 2] = 0,
    title: Optional[str] = None,
    axes: Optional[list[Axes]] = None,
    return_objects: bool = False,
) -> Optional[tuple[Optional[Figure], list[Axes]]]

Plot the chain-length distribution.

PARAMETER	DESCRIPTION
kind	Kind(s) of distribution. TYPE: Literal['number', 'mass', 'gpc'] DEFAULT: 'mass'
sizeasmass	Switch size input between chain-length (if False) or molar mass (if True). TYPE: bool DEFAULT: False
xscale	x-axis scale. TYPE: Literal['auto', 'linear', 'log'] DEFAULT: 'auto'
xrange	x-axis range. TYPE: tuple[float, float] | None DEFAULT: None
cdf	y-axis where cdf is displayed. If 0 the cdf is not displayed; if 1 the cdf is displayed on the primary y-axis; if 2 the cdf is displayed on the secondary axis. TYPE: Literal[0, 1, 2] DEFAULT: 0
title	Title of plot. If None, the object name will be used. TYPE: str | None DEFAULT: None
axes	Matplotlib Axes object. TYPE: list[Axes] | None DEFAULT: None
return_objects	If True, the Figure and Axes objects are returned (for saving or further manipulations). TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
tuple[Figure | None, list[Axes]] | None	Figure and Axes objects if return_objects is True.

Source code in src/polykin/distributions/base.py

def plot(self,
         kind: Union[Literal['number', 'mass', 'gpc'],
                     list[Literal['number', 'mass', 'gpc']]] = 'mass',
         sizeasmass: bool = False,
         xscale: Literal['auto', 'linear', 'log'] = 'auto',
         xrange: Union[tuple[float, float], None] = None,
         cdf: Literal[0, 1, 2] = 0,
         title: Optional[str] = None,
         axes: Optional[list[Axes]] = None,
         return_objects: bool = False
         ) -> Optional[tuple[Optional[Figure], list[Axes]]]:
    """Plot the chain-length distribution.

    Parameters
    ----------
    kind : Literal['number', 'mass', 'gpc']
        Kind(s) of distribution.
    sizeasmass : bool
        Switch size input between chain-length (if `False`) or molar
        mass (if `True`).
    xscale : Literal['auto', 'linear', 'log']
        x-axis scale.
    xrange : tuple[float, float] | None
        x-axis range.
    cdf : Literal[0, 1, 2]
        y-axis where cdf is displayed. If `0` the cdf is not displayed; if
        `1` the cdf is displayed on the primary y-axis; if `2` the cdf is
        displayed on the secondary axis.
    title : str | None
        Title of plot. If `None`, the object name will be used.
    axes : list[Axes] | None
        Matplotlib Axes object.
    return_objects : bool
        If `True`, the Figure and Axes objects are returned (for saving or
        further manipulations).

    Returns
    -------
    tuple[Figure | None, list[Axes]] | None
        Figure and Axes objects if return_objects is `True`.
    """
    # Check inputs
    kind = self._verify_kind(kind, accept_list=True)
    self._verify_sizeasmass(sizeasmass)
    check_in_set(xscale, {'linear', 'log', 'auto'}, 'xscale')
    check_in_set(cdf, {0, 1, 2}, 'cdf')
    if isinstance(kind, str):
        kind = [kind]

    # x-axis scale
    if xscale == 'auto' and set(kind) == {'gpc'}:
        xscale = 'log'
    elif xscale == 'log':
        pass
    else:
        xscale = 'linear'

    # x-axis range
    if xrange is not None:
        check_valid_range(xrange, 0., np.inf, 'xrange')
    else:
        vrange = self._xrange_plot(sizeasmass)  # type : ignore
        if xscale == 'log' and log10(vrange[1]/vrange[0]) > 3 and \
                isinstance(self, (AnalyticalDistribution,
                                  MixtureDistribution)):
            vrange[1] *= 10
        xrange = tuple(vrange)  # type: ignore

    # x-axis vector
    npoints = 200
    if isinstance(self, MixtureDistribution):
        npoints += 100*(len(self.components)-1)
    if xscale == 'log':
        x = np.geomspace(*xrange, npoints)  # type: ignore
    else:
        x = np.linspace(*xrange, npoints)  # type: ignore

    # x-axis label
    if sizeasmass:
        label_x = f"Molar mass [{self.units['molar_mass']}]"
    else:
        label_x = "Chain length"

    # Create axis if none is provided
    ax2 = None
    if axes is None:
        ext_mode = False
        fig, ax = plt.subplots(1, 1)
        if title is None:
            title = f"Distribution: {self.name}"
        fig.suptitle(title)
        if cdf == 2:
            ax2 = ax.twinx()
    else:
        ext_mode = True
        fig = None
        ax = axes[0]
        if cdf == 2:
            ax2 = axes[1]

    # y-values
    for mykind in kind:
        if cdf != 1:
            y1 = self.pdf(x, kind=mykind, sizeasmass=sizeasmass)
        if cdf > 0:
            if mykind == 'gpc':
                _mykind = 'mass'
            else:
                _mykind = mykind
            y2 = self.cdf(x, kind=_mykind, sizeasmass=sizeasmass)
        if cdf == 1:
            y1 = y2
        if ext_mode:
            label = self.name
            if label == '':
                label = '?'
        else:
            label = mykind
        ax.plot(x, y1, label=label)
        if cdf == 2:
            ax2.plot(x, y2, linestyle='--')

    # y-axis and labels
    label_y_pdf = 'Relative abundance'
    label_y_cdf = 'Cumulative probability'
    bbox_to_anchor = (1.05, 1.0)
    if cdf == 0:
        label_y1 = label_y_pdf
    elif cdf == 1:
        label_y1 = label_y_cdf
    elif cdf == 2:
        label_y1 = label_y_pdf
        label_y2 = label_y_cdf
        ax.set_ylabel(label_y1)
        ax2.set_ylabel(label_y2)
        bbox_to_anchor = (1.1, 1.0)
    else:
        raise ValueError
    ax.set_xlabel(label_x)
    ax.set_ylabel(label_y1)
    ax.set_xscale(xscale)
    ax.grid(True)
    ax.legend(bbox_to_anchor=bbox_to_anchor, loc="upper left")

    if return_objects:
        return (fig, [ax, ax2])