Skip to content

Distributions (polykin.distributions)¤

This module implements methods to create, visualize, fit, combine, integrate, etc. theoretical and experimental chain-length distributions.

DataDistribution ¤

Arbitrary numerical chain-length distribution, defined by chain size and pdf data.

PARAMETER DESCRIPTION
size_data

Chain length or molar mass data.

TYPE: FloatVectorLike

pdf_data

Distribution data.

TYPE: FloatVectorLike

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

M0

Molar mass of the repeating unit, \(M_0\). Unit = kg/mol.

TYPE: float DEFAULT: 0.1

name

Name.

TYPE: str DEFAULT: ''

Source code in src/polykin/distributions/datadistribution.py 
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
class DataDistribution(IndividualDistribution):
    r"""Arbitrary numerical chain-length distribution, defined by chain size
    and pdf data.

    Parameters
    ----------
    size_data : FloatVectorLike
        Chain length or molar mass data.
    pdf_data : FloatVectorLike
        Distribution data.
    kind : Literal['number', 'mass', 'gpc']
        Kind of distribution.
    sizeasmass : bool
        Switch size input between chain-length (if `False`) or molar
        mass (if `True`).
    M0 : float
        Molar mass of the repeating unit, $M_0$. Unit = kg/mol.
    name : str
        Name.
    """
    _continuous = True

    def __init__(self,
                 size_data: FloatVectorLike,
                 pdf_data: FloatVectorLike,
                 kind: Literal['number', 'mass', 'gpc'] = 'mass',
                 sizeasmass: bool = False,
                 M0: float = 0.1,
                 name: str = ''
                 ) -> None:

        # Check and clean input
        self.M0 = M0
        self.name = name
        size_data = np.array(size_data)
        pdf_data = np.array(pdf_data)
        if self._verify_sizeasmass(sizeasmass):
            size_data /= self.M0
        idx_valid = np.logical_and(pdf_data > 0., size_data >= 1.)
        if not idx_valid.all():
            print("Warning: Found and removed inconsistent values.")
        self._pdf_data = pdf_data[idx_valid]
        self._length_data = size_data[idx_valid]
        self._pdf_order = self.kind_order[self._verify_kind(kind)]

        # Base-line correction
        # y = self._pdf_data
        # x = self._length_data
        # baseline = y[0] + (y[-1] - y[0])/(np.log(x[-1]/x[0]))*np.log(x/x[0])
        # self._pdf_data -= baseline

        # Compute spline
        self._pdf_spline = \
            interpolate.UnivariateSpline(self._length_data,
                                         self._pdf_data,
                                         k=3,
                                         s=0,
                                         ext=1)

    @vectorize
    def _moment_quadrature(self,
                           xa: float,
                           xb: float,
                           order: int
                           ) -> float:
        xrange = (max(xa, self._length_data[0]),
                  min(xb, self._length_data[-1]))
        if order == self._pdf_order:
            result = self._pdf_spline.integral(*xrange)
        else:
            result, _ = integrate.quad(
                lambda x: x**(order-self._pdf_order)*self._pdf_spline(x),
                *xrange, limit=50)
        return result

    def _pdf0_length(self, x):
        return x**(-self._pdf_order)*self._pdf_spline(x)

    @functools.cached_property
    def _range_length_default(self):
        return self._length_data[(0, -1), ]

    def fit(self,
            dist_class: Union[type[Flory], type[Poisson], type[LogNormal],
                              type[SchulzZimm]],
            dim: int = 1,
            display_table: bool = True
            ) -> Optional[Union[AnalyticalDistribution, MixtureDistribution]]:
        """Fit (deconvolute) a `DataDistribution` into a linear combination of
        `AnalyticalDistribution`(s).

        Parameters
        ----------
        dist_class : type[Flory] | type[Poisson] | type[LogNormal] | type[SchulzZimm]
            Type of distribution to be used in the fit.
        dim : int
            Number of individual components to use in the fit.
        display_table : bool
            Option to display results table with information about individual
            components.

        Returns
        -------
        AnalyticalDistribution | MixtureDistribution | None
            If fit successful, it returns the fitted distribution.
        """

        check_subclass(dist_class, AnalyticalDistribution, 'dist_class')
        isP2 = issubclass(dist_class, AnalyticalDistributionP2)
        check_bounds(dim, 1, 10, 'dim')

        # Init fit distribution
        dfit = MixtureDistribution({})
        weight = np.full(dim, 1/dim)
        DPn = np.empty_like(weight)
        if isP2:
            PDI = np.full(dim, 2.0)
        for i in range(dim):
            DPn[i] = self._length_data[0] * \
                (self._length_data[-1]/self._length_data[0])**((i+1)/(dim+1))
            if isP2:
                args = (DPn[i], PDI[i])
            else:
                args = (DPn[i],)
            d = dist_class(*args, M0=self.M0)
            dfit = dfit + weight[i]*d

        # Define objective function
        xdata = self._length_data
        ydata = self._cdf_length(xdata, 1)

        def objective_fun(x):
            # Assign values
            for i, d in enumerate(dfit.components.keys()):
                dfit.components[d] = x[i]
                d.DPn = 10**x[dim+i]
                if isP2:
                    d.PDI = x[2*dim+i]
            yfit = dfit._cdf(xdata, 1, False)
            return np.sum((yfit - ydata)**2)/ydata.size

        # Initial guess and bounds
        # We do a log transform on DPn to normalize the changes
        x0 = np.concatenate([weight, log10(DPn)])
        bounds = [(0, 1) for _ in range(dim)] + \
            [(log10(3), log10(xdata[-1])) for _ in range(dim)]
        if isP2:
            x0 = np.concatenate([x0, PDI])
            bounds += [(1.01, 5.) for _ in range(dim)]

        # Equality constraint: w(1) + .. + w(N) = 1
        A = np.zeros(x0.size)
        A[:dim] = 1
        constraint = optimize.LinearConstraint(A, 1, 1)
        constraints = []
        constraints.append(constraint)

        # Inequality constraints: DPn(i) - DPn(i+1) <= 0
        for i in range(dim-1):
            A = np.zeros(x0.size)
            A[dim+i] = 1
            A[dim+i+1] = -1
            constraint = optimize.LinearConstraint(A, -np.inf, 0)
            constraints.append(constraint)

        # Call fit method
        # Tried all methods and 'trust-constr' is the most robust
        solution = optimize.minimize(objective_fun,
                                     x0=x0,
                                     method='trust-constr',
                                     bounds=bounds,
                                     constraints=constraints,
                                     options={'verbose': 0})
        if solution.success:
            if display_table:
                print(dfit.components_table)
            if dim == 1:
                dfit = next(iter(dfit.components))
            dfit.name = self.name + "-fit"
            result = dfit
        else:
            print("Failed to fit distribution: ", solution.message)
            result = None

        return result  # type : ignore

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 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 ¤

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 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
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)

fit ¤

fit(
    dist_class: Union[
        type[Flory],
        type[Poisson],
        type[LogNormal],
        type[SchulzZimm],
    ],
    dim: int = 1,
    display_table: bool = True,
) -> Optional[
    Union[AnalyticalDistribution, MixtureDistribution]
]

Fit (deconvolute) a DataDistribution into a linear combination of AnalyticalDistribution(s).

PARAMETER DESCRIPTION
dist_class

Type of distribution to be used in the fit.

TYPE: type[Flory] | type[Poisson] | type[LogNormal] | type[SchulzZimm]

dim

Number of individual components to use in the fit.

TYPE: int DEFAULT: 1

display_table

Option to display results table with information about individual components.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
AnalyticalDistribution | MixtureDistribution | None

If fit successful, it returns the fitted distribution.
Source code in src/polykin/distributions/datadistribution.py 
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
def fit(self,
        dist_class: Union[type[Flory], type[Poisson], type[LogNormal],
                          type[SchulzZimm]],
        dim: int = 1,
        display_table: bool = True
        ) -> Optional[Union[AnalyticalDistribution, MixtureDistribution]]:
    """Fit (deconvolute) a `DataDistribution` into a linear combination of
    `AnalyticalDistribution`(s).

    Parameters
    ----------
    dist_class : type[Flory] | type[Poisson] | type[LogNormal] | type[SchulzZimm]
        Type of distribution to be used in the fit.
    dim : int
        Number of individual components to use in the fit.
    display_table : bool
        Option to display results table with information about individual
        components.

    Returns
    -------
    AnalyticalDistribution | MixtureDistribution | None
        If fit successful, it returns the fitted distribution.
    """

    check_subclass(dist_class, AnalyticalDistribution, 'dist_class')
    isP2 = issubclass(dist_class, AnalyticalDistributionP2)
    check_bounds(dim, 1, 10, 'dim')

    # Init fit distribution
    dfit = MixtureDistribution({})
    weight = np.full(dim, 1/dim)
    DPn = np.empty_like(weight)
    if isP2:
        PDI = np.full(dim, 2.0)
    for i in range(dim):
        DPn[i] = self._length_data[0] * \
            (self._length_data[-1]/self._length_data[0])**((i+1)/(dim+1))
        if isP2:
            args = (DPn[i], PDI[i])
        else:
            args = (DPn[i],)
        d = dist_class(*args, M0=self.M0)
        dfit = dfit + weight[i]*d

    # Define objective function
    xdata = self._length_data
    ydata = self._cdf_length(xdata, 1)

    def objective_fun(x):
        # Assign values
        for i, d in enumerate(dfit.components.keys()):
            dfit.components[d] = x[i]
            d.DPn = 10**x[dim+i]
            if isP2:
                d.PDI = x[2*dim+i]
        yfit = dfit._cdf(xdata, 1, False)
        return np.sum((yfit - ydata)**2)/ydata.size

    # Initial guess and bounds
    # We do a log transform on DPn to normalize the changes
    x0 = np.concatenate([weight, log10(DPn)])
    bounds = [(0, 1) for _ in range(dim)] + \
        [(log10(3), log10(xdata[-1])) for _ in range(dim)]
    if isP2:
        x0 = np.concatenate([x0, PDI])
        bounds += [(1.01, 5.) for _ in range(dim)]

    # Equality constraint: w(1) + .. + w(N) = 1
    A = np.zeros(x0.size)
    A[:dim] = 1
    constraint = optimize.LinearConstraint(A, 1, 1)
    constraints = []
    constraints.append(constraint)

    # Inequality constraints: DPn(i) - DPn(i+1) <= 0
    for i in range(dim-1):
        A = np.zeros(x0.size)
        A[dim+i] = 1
        A[dim+i+1] = -1
        constraint = optimize.LinearConstraint(A, -np.inf, 0)
        constraints.append(constraint)

    # Call fit method
    # Tried all methods and 'trust-constr' is the most robust
    solution = optimize.minimize(objective_fun,
                                 x0=x0,
                                 method='trust-constr',
                                 bounds=bounds,
                                 constraints=constraints,
                                 options={'verbose': 0})
    if solution.success:
        if display_table:
            print(dfit.components_table)
        if dim == 1:
            dfit = next(iter(dfit.components))
        dfit.name = self.name + "-fit"
        result = dfit
    else:
        print("Failed to fit distribution: ", solution.message)
        result = None

    return result  # type : ignore

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 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
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])

Flory ¤

Flory-Schulz (aka most-probable) chain-length distribution.

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

\[ p(k) = (1-a)a^{k-1} \]

where \(a=1-1/DP_n\). Mathematically speaking, this is a geometric distribution.

PARAMETER DESCRIPTION
DPn

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

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 Flory distribution and evaluate the corresponding probability density function and cumulative distribution function for representative chain lengths.

>>> from polykin.distributions import Flory
>>> a = Flory(100, M0=0.050, name='A')
>>> a
type: Flory
name: A
DPn:  100.0
DPw:  199.0
DPz:  298.5
PDI:  1.99
M0:   0.050 kg/mol
Mn:   5.000 kg/mol
Mw:   9.950 kg/mol
Mz:   14.925 kg/mol

>>> a.pdf(a.DPn)
0.003697296376497271

>>> a.cdf([a.DPn, a.DPw, a.DPz])
array([0.26793532, 0.59535432, 0.80159978])
Source code in src/polykin/distributions/analyticaldistributions.py 
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
class Flory(AnalyticalDistributionP1):
    r"""Flory-Schulz (aka most-probable) chain-length distribution.

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

    $$ p(k) = (1-a)a^{k-1} $$

    where $a=1-1/DP_n$. Mathematically speaking, this is a [geometric
    distribution](https://en.wikipedia.org/wiki/Geometric_distribution).

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


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

    >>> from polykin.distributions import Flory
    >>> a = Flory(100, M0=0.050, name='A')
    >>> a
    type: Flory
    name: A
    DPn:  100.0
    DPw:  199.0
    DPz:  298.5
    PDI:  1.99
    M0:   0.050 kg/mol
    Mn:   5.000 kg/mol
    Mw:   9.950 kg/mol
    Mz:   14.925 kg/mol

    >>> a.pdf(a.DPn)
    0.003697296376497271

    >>> a.cdf([a.DPn, a.DPw, a.DPz])
    array([0.26793532, 0.59535432, 0.80159978])

    """
    _continuous = False

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

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

    def _update_internal_parameters(self):
        super()._update_internal_parameters()
        self._a = 1 - 1/self.DPn

    def _pdf0_length(self, k):
        a = self._a
        return (1 - a) * a ** (k - 1)

    @functools.cache
    def _moment_length(self, order):
        a = self._a
        if order == 0:
            result = 1.0
        elif order == 1:
            result = 1/(1 - a)
        elif order == 2:
            result = 2/(1 - a)**2 - 1/(1 - a)
        elif order == 3:
            result = 6/(1 - a)**3 - 6/(1 - a)**2 + 1/(1 - a)
        else:
            raise ValueError("Not defined for order>3.")
        return result

    def _cdf_length(self, k, order):
        a = self._a
        if order == 0:
            result = 1 - a**k
        elif order == 1:
            result = (a**k*(-a*k + k + 1) - 1)/(a - 1)
        else:
            raise ValueError
        return result / self._moment_length(order)

    def _random_length(self, size):
        return self._rng.geometric((1-self._a), size)  # type: ignore

    @functools.cached_property
    def _range_length_default(self):
        return st.geom.ppf(self._ppf_bounds, p=(1-self._a))

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 ¤

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 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
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 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
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 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
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)

LogNormal ¤

Log-normal chain-length distribution.

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

\[ p(x) = \frac{1}{x \sigma \sqrt{2 \pi}} \exp\left (- \frac{(\ln{x}-\mu)^2}{2\sigma^2} \right ) \]

where \(\mu = \ln{(DP_n/\sqrt{PDI})}\) and \(\sigma=\sqrt{\ln(PDI)}\).

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 log-normal distribution and evaluate the corresponding probability density function and cumulative distribution function for representative chain lengths.

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

>>> a.pdf(a.DPn)
0.003317780747597256

>>> a.cdf([a.DPn, a.DPw, a.DPz])
array([0.3001137, 0.6998863, 0.9420503])
Source code in src/polykin/distributions/analyticaldistributions.py 
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
class LogNormal(AnalyticalDistributionP2):
    r"""Log-normal chain-length distribution.

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

    $$ p(x) = \frac{1}{x \sigma \sqrt{2 \pi}}
    \exp\left (- \frac{(\ln{x}-\mu)^2}{2\sigma^2} \right ) $$

    where $\mu = \ln{(DP_n/\sqrt{PDI})}$ and $\sigma=\sqrt{\ln(PDI)}$.

    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 log-normal distribution and evaluate the corresponding probability
    density function and cumulative distribution function for representative
    chain lengths.

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

    >>> a.pdf(a.DPn)
    0.003317780747597256

    >>> a.cdf([a.DPn, a.DPw, a.DPz])
    array([0.3001137, 0.6998863, 0.9420503])

    """
    # https://reference.wolfram.com/language/ref/LogNormalDistribution.html
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html

    _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._sigma = sqrt(log(PDI))
            self._mu = log(DPn/sqrt(PDI))
        except AttributeError:
            pass

    def _pdf0_length(self, x):
        mu = self._mu
        sigma = self._sigma
        return exp(-(log(x) - mu)**2/(2*sigma**2))/(x*sigma*sqrt(2*pi))

    @functools.cache
    def _moment_length(self, order):
        mu = self._mu
        sigma = self._sigma
        return exp(order*mu + 0.5*order**2*sigma**2)

    def _cdf_length(self, x, order):
        mu = self._mu
        sigma = self._sigma
        if order == 0:
            result = (1 + sp.erf((log(x) - mu)/(sigma*sqrt(2))))/2
        elif order == 1:
            result = sp.erfc((mu + sigma**2 - log(x))/(sigma*sqrt(2)))/2
        else:
            raise ValueError
        return result

    def _random_length(self, size):
        mu = self._mu
        sigma = self._sigma
        return np.rint(self._rng.lognormal(mu, sigma, size))  # type: ignore

    @functools.cached_property
    def _range_length_default(self):
        mu = self._mu
        sigma = self._sigma
        return st.lognorm.ppf(self._ppf_bounds, s=sigma, scale=exp(mu), 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 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
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 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
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 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
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)

Poisson ¤

Poisson chain-length distribution.

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

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

where \(a=DP_n-1\).

PARAMETER DESCRIPTION
DPn

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

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 Poisson distribution and evaluate the corresponding probability density function and cumulative distribution function for representative chain lengths.

>>> from polykin.distributions import Poisson
>>> a = Poisson(100, M0=0.050, name='A')
>>> a
type: Poisson
name: A
DPn:  100.0
DPw:  101.0
DPz:  102.0
PDI:  1.01
M0:   0.050 kg/mol
Mn:   5.000 kg/mol
Mw:   5.050 kg/mol
Mz:   5.099 kg/mol

>>> a.pdf(a.DPn)
0.04006147193133002

>>> a.cdf([a.DPn, a.DPw, a.DPz])
array([0.48703481, 0.52669305, 0.56558077])
Source code in src/polykin/distributions/analyticaldistributions.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
class Poisson(AnalyticalDistributionP1):
    r"""Poisson chain-length distribution.

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

    $$ p(k) = \frac{a^{k-1} e^{-a}}{\Gamma(k)} $$

    where $a=DP_n-1$.

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

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

    >>> from polykin.distributions import Poisson
    >>> a = Poisson(100, M0=0.050, name='A')
    >>> a
    type: Poisson
    name: A
    DPn:  100.0
    DPw:  101.0
    DPz:  102.0
    PDI:  1.01
    M0:   0.050 kg/mol
    Mn:   5.000 kg/mol
    Mw:   5.050 kg/mol
    Mz:   5.099 kg/mol

    >>> a.pdf(a.DPn)
    0.04006147193133002

    >>> a.cdf([a.DPn, a.DPw, a.DPz])
    array([0.48703481, 0.52669305, 0.56558077])

    """

    _continuous = False

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

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

    def _update_internal_parameters(self):
        super()._update_internal_parameters()
        self._a = self.DPn - 1

    def _pdf0_length(self, k):
        a = self._a
        return poisson(k, a, False)

    @functools.cache
    def _moment_length(self, order):
        a = self._a
        if order == 0:
            result = 1.0
        elif order == 1:
            result = a + 1
        elif order == 2:
            result = a**2 + 3*a + 1
        elif order == 3:
            result = a**3 + 6*a**2 + 7*a + 1
        else:
            raise ValueError("Not defined for order>3.")
        return result

    def _cdf_length(self, k, order):
        a = self._a
        if order == 0:
            result = sp.gammaincc(k, a)
        elif order == 1:
            result = (a + 1)*sp.gammaincc(k, a) - \
                exp(k*log(a) - a - sp.gammaln(k))
        else:
            raise ValueError
        return result/self._moment_length(order)

    def _random_length(self, size):
        return self._rng.poisson(self._a, size) + 1  # type: ignore

    @functools.cached_property
    def _range_length_default(self):
        return st.poisson.ppf(self._ppf_bounds, mu=self._a, 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 ¤

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 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
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 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
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 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
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)

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 
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
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 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
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 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
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 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
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)

WeibullNycanderGold_pdf ¤

WeibullNycanderGold_pdf(
    k: Union[int, IntArrayLike], v: float, r: float
) -> Union[float, FloatArray]

Weibull, Nycander and Golds's analytical chain-length distribution for living polymerization with different initiation and polymerization rate coefficients.

For a living polymerization with only initiation and propagation (i.e., constant number of chains), the number fraction of chains of length \(k\) can computed in two steps. First, the number fraction of unreacted initiator molecules, \(p_0=p(0)\), is found by solving the equation:

\[ v + r \ln(p_0) + (r - 1)(1 - p_0) = 0 \]

where \(v\) denotes the number-average degree of polymerization of all chains, including unreacted initiator molecules, and \(r=k_p/k_i\) is the ratio of the polymerization and initiation rate coefficients. Then, the number fraction of chains with \(k \ge 1\) monomer units can be evaluated by:

\[\begin{aligned} p(k) & = p_0 \frac{r^{k-1}}{(r-1)^k} \left[1- \Gamma(k,(1-r)\ln{p_0}) \right] & \text{if } r>1 \\ p(k) & = p_0 \frac{r^{k-1}}{(1-r)^k} (-1)^k\left[1- \Gamma(k,(1-r)\ln{p_0}) \right] & \text{if } r<1 \end{aligned}\]

where \(\Gamma\) is the regularized upper incomplete gamma function. For \(r>1\), the argument of \(\Gamma\) is always positive, while for \(r<1\) it is negative. This analytical solution has an obvious singularity at \(r=1\); in that case, the solution reduces to the well-known Poisson distribution:

\[ p(k) = \frac{v^k}{k!}e^{-v} \]

valid for \(k \ge 0\).

Note

  • The solution is numerically unstable in certain domains, namely for \(r\) close to 1, and also for \(k>>v\). This is an intrinsic feature of the equation.
  • For \(|r-1|<10^{-2}\), the algorithm automatically switches to the Poisson distribution. Some numerical discontinuity at this boundary is to be expected.
  • For \(r<1\), no solution is currently computed, because the incomplete gamma function algorithm available in SciPy is restricted to positive arguments.

References

  • Weibull, B.; Nycander, E.. "The Distribution of Compounds Formed in the Reaction." Acta Chemica Scandinavica 49 (1995): 207-216.
  • Gold, L. "Statistics of polymer molecular size distribution for an invariant number of propagating chains." The Journal of Chemical Physics 28.1 (1958): 91-99.
PARAMETER DESCRIPTION
k

Chain length (>=0).

TYPE: int | IntArrayLike

v

Number-average degree of polymerization considering chains with zero length.

TYPE: float

r

Ratio of propagation and initiation rate coefficients.

TYPE: float

RETURNS DESCRIPTION
float | FloatArray

Number probability density.

Examples:

Compute the fraction of chains with lengths 0 to 2 for a system with \(r=5\) and \(v=1\).

>>> from polykin.distributions import WeibullNycanderGold_pdf
>>> WeibullNycanderGold_pdf([0, 1, 2], 1., 5)
array([0.58958989, 0.1295864 , 0.11493254])
Source code in src/polykin/distributions/analyticaldistributions.py 
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
def WeibullNycanderGold_pdf(k: Union[int, IntArrayLike],
                            v: float,
                            r: float
                            ) -> Union[float, FloatArray]:
    r"""Weibull, Nycander and Golds's analytical chain-length distribution for
    living polymerization with different initiation and polymerization rate
    coefficients.

    For a living polymerization with only initiation and propagation (i.e.,
    constant number of chains), the number fraction of chains of length
    $k$ can computed in two steps. First, the number fraction of unreacted 
    initiator molecules, $p_0=p(0)$, is found by solving the equation:

    $$ v + r \ln(p_0) + (r - 1)(1 - p_0) = 0 $$

    where $v$ denotes the number-average degree of polymerization of all chains,
    including unreacted initiator molecules, and $r=k_p/k_i$ is the ratio of the
    polymerization and initiation rate coefficients. Then, the number fraction
    of chains with $k \ge 1$ monomer units can be evaluated by:

    \begin{aligned}
    p(k) & = p_0 \frac{r^{k-1}}{(r-1)^k} \left[1- \Gamma(k,(1-r)\ln{p_0}) \right] & \text{if } r>1 \\
    p(k) & = p_0 \frac{r^{k-1}}{(1-r)^k} (-1)^k\left[1- \Gamma(k,(1-r)\ln{p_0}) \right] & \text{if } r<1
    \end{aligned}

    where $\Gamma$ is the regularized upper incomplete gamma function. For $r>1$,
    the argument of $\Gamma$ is always positive, while for $r<1$ it is negative. 
    This analytical solution has an obvious singularity at $r=1$; in that case,
    the solution reduces to the well-known Poisson distribution:

    $$ p(k) = \frac{v^k}{k!}e^{-v} $$

    valid for $k \ge 0$.

    !!! note

        * The solution is numerically unstable in certain domains, namely for
        $r$ close to 1, and also for $k>>v$. This is an intrinsic feature  of
        the equation.
        * For $|r-1|<10^{-2}$, the algorithm automatically switches to the Poisson
        distribution. Some numerical discontinuity at this boundary is to be
        expected. 
        * For $r<1$, no solution is currently computed, because the incomplete
        gamma function algorithm available in SciPy is restricted to positive
        arguments.    

    **References**

    * Weibull, B.; Nycander, E.. "The Distribution of Compounds Formed in the
    Reaction." Acta Chemica Scandinavica 49 (1995): 207-216.
    * Gold, L. "Statistics of polymer molecular size distribution for an
    invariant number of propagating chains." The Journal of Chemical Physics
    28.1 (1958): 91-99.

    Parameters
    ----------
    k : int | IntArrayLike
        Chain length (>=0).
    v : float
        Number-average degree of polymerization considering chains with zero
        length.
    r : float
        Ratio of propagation and initiation rate coefficients.

    Returns
    -------
    float | FloatArray
        Number probability density.

    Examples
    --------
    Compute the fraction of chains with lengths 0 to 2 for a system with $r=5$
    and $v=1$.

    >>> from polykin.distributions import WeibullNycanderGold_pdf
    >>> WeibullNycanderGold_pdf([0, 1, 2], 1., 5)
    array([0.58958989, 0.1295864 , 0.11493254])
    """

    if isinstance(k, (list, tuple)):
        k = np.asarray(k)

    # Special case where it reduces to Poisson(kmin=0)
    if abs(r - 1.) < 1e-2:
        return poisson(k, v, True)

    # Find p0=p(k=0)
    def find_p0(p0):
        return v + r*log(p0) + (r - 1.)*(1. - p0)

    sol = root_scalar(find_p0, method='brentq', bracket=(0. + eps, 1. - eps))
    p0 = sol.root

    # Avoid issue 'Integers to negative integer powers'
    r = float(r)

    # Weibull-Nycander-Gold solution p(k>=1)
    if r > 1.:
        A = sp.gammaincc(k, (1. - r)*log(p0))
        result = p0 * r**(k - 1) / (r - 1.)**k * (1. - A)
    else:
        # gammaincc in scipy is restricted to positive arguments
        # A = sp.gammaincc(k, (1. - r)*log(p0))
        # result = p0 * r**(k - 1) / (1. - r)**k * (-1)**k * (1. - A)
        result = NotImplemented

    # Overwite p(k=0)
    result[k == 0] = p0

    return result

convolve_moments ¤

convolve_moments(
    q0: float,
    q1: float,
    q2: float,
    r0: float,
    r1: float,
    r2: float,
) -> tuple[float, float, float]

Compute the first three moments of the convolution of two distributions.

If \(P = Q * R\) is the convolution of \(Q\) and \(R\), defined as:

\[ P_n = \sum_{i=0}^{n} Q_{n-i}R_{i} \]

then the first three moments of \(P\) are related to the moments of \(Q\) and \(R\) by:

\[\begin{aligned} p_0 &= q_0 r_0 \\ p_1 &= q_1 r_0 + q_0 r_1 \\ p_2 &= q_2 r_0 + 2 q_1 r_1 + q_0 r_2 \end{aligned}\]

where \(p_i\), \(q_i\) and \(r_i\) denote the \(i\)-th moments of \(P\), \(Q\) and \(R\), respectively.

PARAMETER DESCRIPTION
q0

0-th moment of \(Q\).

TYPE: float

q1

1-st moment of \(Q\).

TYPE: float

q2

2-nd moment of \(Q\).

TYPE: float

r0

0-th moment of \(R\).

TYPE: float

r1

1-st moment of \(R\).

TYPE: float

r2

2-nd moment of \(R\).

TYPE: float

RETURNS DESCRIPTION
tuple[float, float, float]

0-th, 1-st and 2-nd moments of \(P=Q*R\).

Examples:

>>> from polykin.distributions import convolve_moments
>>> convolve_moments(1., 1e2, 2e4, 1., 50., 5e4)
(1.0, 150.0, 80000.0)
Source code in src/polykin/distributions/base.py 
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
def convolve_moments(q0: float,
                     q1: float,
                     q2: float,
                     r0: float,
                     r1: float,
                     r2: float
                     ) -> tuple[float, float, float]:
    r"""Compute the first three moments of the convolution of two distributions.

    If $P = Q * R$ is the convolution of $Q$ and $R$, defined as:

    $$ P_n = \sum_{i=0}^{n} Q_{n-i}R_{i} $$

    then the first three moments of $P$ are related to the moments of $Q$ and
    $R$ by:

    \begin{aligned}
    p_0 &= q_0 r_0 \\
    p_1 &= q_1 r_0 + q_0 r_1 \\
    p_2 &= q_2 r_0 + 2 q_1 r_1 + q_0 r_2
    \end{aligned}

    where $p_i$, $q_i$ and $r_i$ denote the $i$-th moments of $P$, $Q$ and $R$,
    respectively.    

    Parameters
    ----------
    q0 : float
        0-th moment of $Q$.
    q1 : float
        1-st moment of $Q$.
    q2 : float
        2-nd moment of $Q$.
    r0 : float
        0-th moment of $R$.
    r1 : float
        1-st moment of $R$.
    r2 : float
        2-nd moment of $R$.

    Returns
    -------
    tuple[float, float, float]
        0-th, 1-st and 2-nd moments of $P=Q*R$.

    Examples
    --------
    >>> from polykin.distributions import convolve_moments
    >>> convolve_moments(1., 1e2, 2e4, 1., 50., 5e4)
    (1.0, 150.0, 80000.0)
    """
    p0 = q0*r0
    p1 = q1*r0 + q0*r1
    p2 = q2*r0 + 2*q1*r1 + q0*r2
    return p0, p1, p2

convolve_moments_self ¤

convolve_moments_self(
    q0: float, q1: float, q2: float, order: int = 1
) -> tuple[float, float, float]

Compute the first three moments of the k-th order convolution of a distribution with itself.

If \(P^k\) is the \(k\)-th order convolution of \(Q\) with itself, defined as:

\[\begin{aligned} P^1_n &= Q*Q = \sum_{i=0}^{n} Q_{n-i} Q_{i} \\ P^2_n &= (Q*Q)*Q = \sum_{i=0}^{n} Q_{n-i} P^1_{i} \\ P^3_n &= ((Q*Q)*Q)*Q = \sum_{i=0}^{n} Q_{n-i} P^2_{i} \\ ... \end{aligned}\]

then the first three moments of \(P^k\) are related to the moments of \(Q\) by:

\[\begin{aligned} p_0 &= q_0^{k+1} \\ p_1 &= (k+1) q_0^k q_1 \\ p_2 &= (k+1) q_0^{k-1} (k q_1^2 +q_0 q_2) \end{aligned}\]

where \(p_i\) and \(q_i\) denote the \(i\)-th moments of \(P^k\) and \(Q\), respectively.

PARAMETER DESCRIPTION
q0

0-th moment of \(Q\).

TYPE: float

q1

1-st moment of \(Q\).

TYPE: float

q2

2-nd moment of \(Q\).

TYPE: float

RETURNS DESCRIPTION
tuple[float, float, float]

0-th, 1-st and 2-nd moments of \(P^k=(Q*Q)*...\).

Examples:

>>> from polykin.distributions import convolve_moments_self
>>> convolve_moments_self(1., 1e2, 2e4, 2)
(1.0, 300.0, 120000.0)
Source code in src/polykin/distributions/base.py 
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
def convolve_moments_self(q0: float,
                          q1: float,
                          q2: float,
                          order: int = 1
                          ) -> tuple[float, float, float]:
    r"""Compute the first three moments of the k-th order convolution of a
    distribution with itself.

    If $P^k$ is the $k$-th order convolution of $Q$ with itself, defined as:

    \begin{aligned}
    P^1_n &= Q*Q = \sum_{i=0}^{n} Q_{n-i} Q_{i} \\
    P^2_n &= (Q*Q)*Q = \sum_{i=0}^{n} Q_{n-i} P^1_{i} \\
    P^3_n &= ((Q*Q)*Q)*Q = \sum_{i=0}^{n} Q_{n-i} P^2_{i} \\
    ...
    \end{aligned}

    then the first three moments of $P^k$ are related to the moments of $Q$ by:

    \begin{aligned}
    p_0 &= q_0^{k+1}  \\
    p_1 &= (k+1) q_0^k q_1 \\
    p_2 &= (k+1) q_0^{k-1} (k q_1^2 +q_0 q_2)
    \end{aligned}

    where $p_i$ and $q_i$ denote the $i$-th moments of $P^k$ and $Q$,
    respectively.    

    Parameters
    ----------
    q0 : float
        0-th moment of $Q$.
    q1 : float
        1-st moment of $Q$.
    q2 : float
        2-nd moment of $Q$.

    Returns
    -------
    tuple[float, float, float]
        0-th, 1-st and 2-nd moments of $P^k=(Q*Q)*...$.

    Examples
    --------
    >>> from polykin.distributions import convolve_moments_self
    >>> convolve_moments_self(1., 1e2, 2e4, 2)
    (1.0, 300.0, 120000.0)
    """
    p0 = q0**(order+1)
    p1 = (order+1) * q0**order * q1
    p2 = (order+1) * q0**(order-1) * (order*q1**2 + q0*q2)
    return p0, p1, p2

plotdists ¤

plotdists(
    dists: list[Distribution],
    kind: Literal["number", "mass", "gpc"],
    title: Optional[str] = None,
    **kwargs
) -> Figure

Plot a list of distributions in a joint plot.

PARAMETER DESCRIPTION
dists

List of distributions to be ploted together.

TYPE: list[Distribution]

kind

Kind of distribution.

TYPE: Literal['number', 'mass', 'gpc']

title

Title of plot.

TYPE: str | None DEFAULT: None

RETURNS DESCRIPTION
Figure

Matplotlib Figure object holding the joint plot.

Examples:

>>> from polykin.distributions import Flory, LogNormal, plotdists
>>> a = Flory(100, M0=0.050, name='A')
>>> b = LogNormal(100, PDI=3., M0=0.050, name='B')
>>> fig = plotdists([a, b], kind='gpc', xrange=(1, 1e4), cdf=2)
>>> fig.show()
Source code in src/polykin/distributions/base.py 
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
def plotdists(dists: list[Distribution],
              kind: Literal['number', 'mass', 'gpc'],
              title: Optional[str] = None,
              **kwargs
              ) -> Figure:
    """Plot a list of distributions in a joint plot.

    Parameters
    ----------
    dists : list[Distribution]
        List of distributions to be ploted together.
    kind : Literal['number', 'mass', 'gpc']
        Kind of distribution.
    title : str | None
        Title of plot.

    Returns
    -------
    Figure
        Matplotlib Figure object holding the joint plot.

    Examples
    --------
    >>> from polykin.distributions import Flory, LogNormal, plotdists
    >>> a = Flory(100, M0=0.050, name='A')
    >>> b = LogNormal(100, PDI=3., M0=0.050, name='B')
    >>> fig = plotdists([a, b], kind='gpc', xrange=(1, 1e4), cdf=2)
    >>> fig.show()
    """

    # Check input
    kind = Distribution._verify_kind(kind)

    # Create matplotlib objects
    fig, ax = plt.subplots(1, 1)
    if kwargs.get('cdf', 1) == 2:
        ax.twinx()

    # Title
    if title is None:
        titles = {'number': 'Number', 'mass': 'Mass', 'gpc': 'GPC'}
        title = f"{titles.get(kind,'')} distributions"
    fig.suptitle(title)

    # Draw plots sequentially
    for d in dists:
        d.plot(kind=kind, axes=fig.axes, **kwargs)

    return fig

Tutorial