Mathematical helper functions#
Helper functions for mathematical and physical operations.
- maicos.lib.math.FT(t, x, indvar=True)[source]#
Discrete fast fourier transform.
Takes the time series and the function as arguments. By default, returns the FT and the frequency: setting indvar=False means the function returns only the FT.
- maicos.lib.math.center_cluster(ag, weights)[source]#
Calculate the center of the atomgroup with respect to some weights.
- Parameters:
ag (MDAnalysis.core.groups.AtomGroup) – Group of atoms to calculate the center for.
weights (numpy.ndarray) – Weights in the shape of ag.
- Returns:
com – The center with respect to the weights.
- Return type:
Without proper treatment of periodic boundrary conditions most algorithms will result in wrong center of mass calculations where molecules or clusters of particles are broken over the boundrary.
Example
+-----------+ | | | 1 x 2 | | | +-----------+
Following
Linge Bai & David Breen (2008) Calculating Center of Mass in an Unbounded 2D Environment, Journal of Graphics Tools, 13:4, 53-60, DOI: 10.1080/2151237X.2008.10129266
the coordinates of the particles are projected on a circle and weighted by their mass in this two dimensional space. The center of mass is obtained by transforming this point back to the corresponding point in the real system. This is done seperately for each dimension.
Reasons for doing this include the analysis of clusters in periodic boundrary conditions and consistent center of mass calculation across box boundraries. This procedure results in the right center of mass as seen below
+-----------+ | | x 1 2 | | | +-----------+
- maicos.lib.math.compute_form_factor(q, atom_type)[source]#
Calculate the form factor for the given element for given q (1/Å).
Handles united atom types like CH4 etc …
- maicos.lib.math.correlation(a, b=None, subtract_mean=False)[source]#
Calculate correlation or autocorrelation.
Uses fast fourier transforms to give the correlation function of two arrays, or, if only one array is given, the autocorrelation. Setting subtract_mean=True causes the mean to be subtracted from the input data.
- maicos.lib.math.correlation_time(x_n, method='sokal', c=8, mintime=3)[source]#
Compute the integrated correlation time of a timeseries.
- Parameters:
x_n (numpy.ndarray, float) – timeseries
method (str, {'sokal', 'chodera'}) – Method to choose integration cutoff.
c (float) – cut-off factor for calculation of correlation time \(\tau\) for the Sokal method. The cut-off \(T\) for integration is determined to be \(T >= c \cdot tau\).
mintime (int) – minimum possible value for cut-off
- Returns:
tau – integrated correlation time
- Return type:
- Raises:
ValueError – If method is not one of ‘Sokal’ or ‘Chodera’
- maicos.lib.math.iFT(k, xf, indvar=True)[source]#
Inverse discrete fast fourier transform.
Takes the frequency series and the function as arguments. By default, returns the iFT and the time series: setting indvar=False means the function returns only the iFT.
- maicos.lib.math.new_mean(old_mean, data, length)[source]#
Compute the arithmetic mean of a series iteratively.
Compute the arithmetic mean of n samples based on an existing mean of n-1 and the n-th value.
Given the mean of a data series
\[\bar x_N = \frac{1}{N} \sum_{n=1}^N x_n\]we seperate the last value
\[\bar x_N = \frac{1}{N} \sum_{n=1}^{N-1} x_n + \frac{x_N}{N}\]and multiply 1 = (N - 1)/(N - 1)
\[\begin{split}\bar x_N = \frac{N-1}{N} \frac{1}{N-1} \\ \sum_{n=1}^{N-1} x_n + \frac{x_N}{N}\end{split}\]The first term can be identified as the mean of the first N - 1 values and we arrive at
\[\bar x_N = \frac{N-1}{N} \bar x_{N-1} + \frac{x_N}{N}\]- Parameters:
- Returns:
new_mean – Updated mean of the series of n values.
- Return type:
Examples
The mean of a data set can easily be calculated from the data points. However this requires one to keep all data points on hand until the end of the calculation.
>>> np.mean([1,3,5,7]) 4.0
Alternatively, one can update an existing mean, this requires only knowledge of the total number of samples.
>>> maicos.utils.new_mean(np.mean([1, 3, 5]), 7, 4) 4.0
- maicos.lib.math.new_variance(old_variance, old_mean, new_mean, data, length)[source]#
Calculate the variance of a timeseries iteratively.
The variance of a timeseries \(x_n\) can be calculated iteratively by using the following formula:
\[S_n = S_n-1 + (n-1) * (x_n - \bar{x}_n-1)^2 / (n-1)\]Here, \(\bar{x}_n\) is the mean of the timeseries up to the \(n\)-th value.
Floating point imprecision can lead to slight negative variances leading non defined standard deviations. Therefore a negetaive variance is set to 0.
- Parameters:
- Returns:
new_variance – Updated variance of the series of n values.
- Return type:
Examples
The data set
[1,5,5,1]has a variance of4.0>>> np.var([1,5,5,1]) 4.0
Knowing the total number of data points, this operation can be performed iteratively.
>>> maicos.utils.new_variance(np.var([1,5,5]), 1, 4) 4.0
- maicos.lib.math.scalar_prod_corr(a, b=None, subtract_mean=False)[source]#
Give the corr. function of the scalar product of two vector timeseries.
Arguments should be given in the form a[t, i], where t is the time variable along which the correlation is calculated, and i indexes the vector components.
- maicos.lib.math.symmetrize(m, axis=None, inplace=False)[source]#
Symmeterize an array.
The shape of the array is preserved, but the elements are symmetrized with respect to the given axis.
- Parameters:
m (array_like) – Input array to symmetrize
axis (None or int or tuple of ints) – Axis or axes along which to symmetrize over. The default, axis=None, will symmetrize over all of the axes of the input array. If axis is negative it counts from the last to the first axis. If axis is a tuple of ints, symmetrizing is performed on all of the axes specified in the tuple.
inplace (bool) – Do symmetrizations inplace. If False a new array is returned.
- Returns:
out – the symmetrized array
- Return type:
array_like
Notes
symmetrize uses
np.flip()for flipping the indices.Examples
>>> A = np.arange(10).astype(float) >>> A array([0., 1., 2., 3., 4., 5., 6., 7., 8., 9.]) >>> maicos.utils.symmetrize(A) array([4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5]) >>> maicos.utils.symmetrize(A, inplace=True) array([4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5]) >>> A array([4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5])
It also works for arrays with more than 1 dimensions in a general dimension.
>>> A = np.arange(20).astype(float).reshape(2,10).T >>> A array([[ 0., 10.], [ 1., 11.], [ 2., 12.], [ 3., 13.], [ 4., 14.], [ 5., 15.], [ 6., 16.], [ 7., 17.], [ 8., 18.], [ 9., 19.]]) >>> maicos.utils.symmetrize(A) array([[9.5, 9.5], [9.5, 9.5], [9.5, 9.5], [9.5, 9.5], [9.5, 9.5], [9.5, 9.5], [9.5, 9.5], [9.5, 9.5], [9.5, 9.5], [9.5, 9.5]]) >>> maicos.utils.symmetrize(A, axis=0) array([[ 4.5, 14.5], [ 4.5, 14.5], [ 4.5, 14.5], [ 4.5, 14.5], [ 4.5, 14.5], [ 4.5, 14.5], [ 4.5, 14.5], [ 4.5, 14.5], [ 4.5, 14.5], [ 4.5, 14.5]])
- maicos.lib._cutil.compute_structure_factor()#
Calculates \(S(\vert q \vert)\) for all possible \(q\) values.
Returns the \(q\) values as well as the scattering factor.
Use via from maicos.lib.math import compute_structure_factor
- Parameters:
positions (numpy.ndarray) – position array.
boxdimensions (numpy.ndarray) – dimensions of the cell.
startq (float) – Starting q (1/Å).
endq (float) – Ending q (1/Å).
mintheta (float) – Minimal angle (°) between the q vectors and the z-axis.
maxtheta (float) – Maximal angle (°) between the q vectors and the z-axis.