polykin.distributions

SchulzZimm

Schulz-Zimm chain-length distribution.

This distribution is based on the following number probability density function:

\[ p(x) = \frac{x^{k-1} e^{-x/\theta}}{\Gamma(k) \theta^k} \]

where \(k = 1/(DP_n-1)\) and \(\theta = DP_n(PDI-1)\). Mathematically speaking, this is a Gamma distribution.

PARAMETER	DESCRIPTION
DPn	Number-average degree of polymerization, \(DP_n\). TYPE: float
PDI	Polydispersity index, \(PDI\). TYPE: float
M0	Molar mass of the repeating unit, \(M_0\). Unit = kg/mol. TYPE: float DEFAULT: 0.1
name	Name. TYPE: str DEFAULT: ''

Examples:

Define a Schulz-Zimm distribution and evaluate the corresponding probability density function and cumulative distribution function for representative chain lengths.

>>> from polykin.distributions import SchulzZimm
>>> a = SchulzZimm(100, PDI=3., M0=0.050, name='A')
>>> a
type: SchulzZimm
name: A
DPn:  100.0
DPw:  300.0
DPz:  500.0
PDI:  3.00
M0:   0.050 kg/mol
Mn:   5.000 kg/mol
Mw:   15.000 kg/mol
Mz:   25.000 kg/mol

>>> a.pdf(a.DPn)
0.0024197072451914337

>>> a.cdf([a.DPn, a.DPw, a.DPz])
array([0.19874804, 0.60837482, 0.82820286])

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

class SchulzZimm(AnalyticalDistributionP2):
    r"""Schulz-Zimm chain-length distribution.

    This distribution is based on the following number probability density
    function:

    $$ p(x) = \frac{x^{k-1} e^{-x/\theta}}{\Gamma(k) \theta^k} $$

    where $k = 1/(DP_n-1)$ and $\theta = DP_n(PDI-1)$. Mathematically speaking,
    this is a
    [Gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution).

    Parameters
    ----------
    DPn : float
        Number-average degree of polymerization, $DP_n$.
    PDI : float
        Polydispersity index, $PDI$.
    M0 : float
        Molar mass of the repeating unit, $M_0$. Unit = kg/mol.
    name : str
        Name.

    Examples
    --------
    Define a Schulz-Zimm distribution and evaluate the corresponding
    probability density function and cumulative distribution function for
    representative chain lengths.

    >>> from polykin.distributions import SchulzZimm
    >>> a = SchulzZimm(100, PDI=3., M0=0.050, name='A')
    >>> a
    type: SchulzZimm
    name: A
    DPn:  100.0
    DPw:  300.0
    DPz:  500.0
    PDI:  3.00
    M0:   0.050 kg/mol
    Mn:   5.000 kg/mol
    Mw:   15.000 kg/mol
    Mz:   25.000 kg/mol

    >>> a.pdf(a.DPn)
    0.0024197072451914337

    >>> a.cdf([a.DPn, a.DPw, a.DPz])
    array([0.19874804, 0.60837482, 0.82820286])

    """
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html
    # https://goldbook.iupac.org/terms/view/S05502

    _continuous = True

    def __init__(self,
                 DPn: float,
                 PDI: float,
                 M0: float = 0.1,
                 name: str = ''
                 ) -> None:

        super().__init__(DPn, PDI, M0, name)

    def _update_internal_parameters(self):
        super()._update_internal_parameters()
        try:
            PDI = self.PDI
            DPn = self.DPn
            self._k = 1/(PDI - 1)
            self._theta = DPn*(PDI - 1)
        except AttributeError:
            pass

    def _pdf0_length(self, x):
        k = self._k
        theta = self._theta
        return x**(k-1)*exp(-x/theta)/(theta**k*sp.gamma(k))

    @functools.cache
    def _moment_length(self, order):
        k = self._k
        theta = self._theta
        return sp.poch(k, order)*theta**order

    def _cdf_length(self, x, order):
        k = self._k
        theta = self._theta
        if order == 0:
            result = sp.gammainc(k, x/theta)
        elif order == 1:
            result = 1 - sp.gammaincc(1+k, x/theta)
        else:
            raise ValueError
        return result

    def _random_length(self, size):
        k = self._k
        theta = self._theta
        return np.rint(self._rng.gamma(k, theta, size))  # type: ignore

    @functools.cached_property
    def _range_length_default(self):
        k = self._k
        theta = self._theta
        return st.gamma.ppf(self._ppf_bounds, a=k, scale=theta, loc=1)

DPn property writable

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 instance-attribute

M0 = M0

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 writable

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)

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])

random

random(
    shape: Optional[Union[int, tuple[int, ...]]] = None
) -> Union[int, IntArray]

Generate random sample of chain lengths according to the corresponding number probability density/mass function.

PARAMETER	DESCRIPTION
shape	Sample shape. TYPE: int | tuple[int, ...] | None DEFAULT: None

RETURNS	DESCRIPTION
int | IntArray	Random sample of chain lengths.

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

def random(self,
           shape: Optional[Union[int, tuple[int, ...]]] = None
           ) -> Union[int, IntArray]:
    r"""Generate random sample of chain lengths according to the
    corresponding number probability density/mass function.

    Parameters
    ----------
    shape : int | tuple[int, ...] | None
        Sample shape.

    Returns
    -------
    int | IntArray
        Random sample of chain lengths.
    """
    if self._rng is None:
        self._rng = np.random.default_rng()
    return self._random_length(shape)