pyemma.coordinates.tica¶

pyemma.coordinates.tica(data=None, lag=10, dim=-1, var_cutoff=0.95, kinetic_map=True, commute_map=False, weights='empirical', stride=1, remove_mean=True, skip=0, reversible=True, ncov_max=inf)¶

Time-lagged independent component analysis (TICA).

TICA is a linear transformation method. In contrast to PCA, which finds coordinates of maximal variance, TICA finds coordinates of maximal autocorrelation at the given lag time. Therefore, TICA is useful in order to find the slow components in a dataset and thus an excellent choice to transform molecular dynamics data before clustering data for the construction of a Markov model. When the input data is the result of a Markov process (such as thermostatted molecular dynamics), TICA finds in fact an approximation to the eigenfunctions and eigenvalues of the underlying Markov operator [1].

It estimates a TICA transformation from data. When input data is given as an argument, the estimation will be carried out straight away, and the resulting object can be used to obtain eigenvalues, eigenvectors or project input data onto the slowest TICA components. If no data is given, this object is an empty estimator and can be put into a pipeline() in order to use TICA in the streaming mode.

Parameters:	data (ndarray (T, d) or list of ndarray (T_i, d) or a reader created by) – source function array with the data, if available. When given, the TICA transformation is immediately computed and can be used to transform data. lag (int, optional, default = 10) – the lag time, in multiples of the input time step dim (int, optional, default -1) – the number of dimensions (independent components) to project onto. A call to the `map` function reduces the d-dimensional input to only dim dimensions such that the data preserves the maximum possible autocorrelation amongst dim-dimensional linear projections. -1 means all numerically available dimensions will be used unless reduced by var_cutoff. Setting dim to a positive value is exclusive with var_cutoff. var_cutoff (float in the range [0,1], optional, default 0.95) – Determines the number of output dimensions by including dimensions until their cumulative kinetic variance exceeds the fraction subspace_variance. var_cutoff=1.0 means all numerically available dimensions (see epsilon) will be used, unless set by dim. Setting var_cutoff smaller than 1.0 is exclusive with dim kinetic_map (bool, optional, default True) – Eigenvectors will be scaled by eigenvalues. As a result, Euclidean distances in the transformed data approximate kinetic distances [4]. This is a good choice when the data is further processed by clustering. commute_map (bool, optional, default False) – Eigenvector_i will be scaled by sqrt(timescale_i / 2). As a result, Euclidean distances in the transformed data will approximate commute distances [5]. stride (int, optional, default = 1) – If set to 1, all input data will be used for estimation. Note that this could cause this calculation to be very slow for large data sets. Since molecular dynamics data is usually correlated at short timescales, it is often sufficient to estimate transformations at a longer stride. Note that the stride option in the get_output() function of the returned object is independent, so you can parametrize at a long stride, and still map all frames through the transformer. weights (optional, default="empirical") – Re-weighting strategy to be used in order to compute equilibrium covariances from non-equilibrium data. “empirical”: no re-weighting “koopman”: use re-weighting procedure from [6] weights: An object that allows to compute re-weighting factors. It must possess a method weights(X) that accepts a trajectory X (np.ndarray(T, n)) and returns a vector of re-weighting factors (np.ndarray(T,)). remove_mean (bool, optional, default True) – remove mean during covariance estimation. Should not be turned off. skip (int, default=0) – skip the first initial n frames per trajectory. reversible (bool, default=True) – symmetrize correlation matrices C_0, C_{tau}. ncov_max (int, default=infinity) – limit the memory usage of the algorithm from [7] to an amount that corresponds to ncov_max additional copies of each correlation matrix
Returns:	tica – Object for time-lagged independent component (TICA) analysis. it contains TICA eigenvalues and eigenvectors, and the projection of input data to the dominant TICA
Return type:	a `TICA` transformation object

Notes

Given a sequence of multivariate data \(X_t\), it computes the mean-free covariance and time-lagged covariance matrix:

\[\begin{split}C_0 &= (X_t - \mu)^T \mathrm{diag}(w) (X_t - \mu) \\ C_{\tau} &= (X_t - \mu)^T \mathrm{diag}(w) (X_t + \tau - \mu)\end{split}\]

where w is a vector of weights for each time step. By default, these weights are all equal to one, but different weights are possible, like the re-weighting to equilibrium described in [6]. Subsequently, the eigenvalue problem

\[C_{\tau} r_i = C_0 \lambda_i r_i,\]

is solved,where \(r_i\) are the independent components and \(\lambda_i\) are their respective normalized time-autocorrelations. The eigenvalues are related to the relaxation timescale by

\[t_i = -\frac{\tau}{\ln |\lambda_i|}.\]

When used as a dimension reduction method, the input data is projected onto the dominant independent components.

TICA was originally introduced for signal processing in [2]. It was introduced to molecular dynamics and as a method for the construction of Markov models in [1] and [3]. It was shown in [1] that when applied to molecular dynamics data, TICA is an approximation to the eigenvalues and eigenvectors of the true underlying dynamics.

Examples

Invoke TICA transformation with a given lag time and output dimension:

>>> import numpy as np
>>> from pyemma.coordinates import tica
>>> data = np.random.random((100,3))
>>> projected_data = tica(data, lag=2, dim=1).get_output()[0]

For a brief explaination why TICA outperforms PCA to extract a good reaction coordinate have a look here.

See also

TICA : tica object

pca : for principal component analysis

class pyemma.coordinates.transform.tica.TICA(lag, dim=-1, var_cutoff=0.95, kinetic_map=True, commute_map=False, epsilon=1e-06, stride=1, skip=0, reversible=True, weights=None, ncov_max=inf)¶

Time-lagged independent component analysis (TICA)

Methods

`describe`()	Get a descriptive string representation of this class.
`dimension`()	output dimension
`estimate`(X, **kwargs)	Chunk-based parameterization of TICA.
`fit`(X)	Estimates parameters - for compatibility with sklearn.
`fit_transform`(X[, y])	Fit to data, then transform it.
`get_output`([dimensions, stride, skip, chunk])	Maps all input data of this transformer and returns it as an array or list of arrays
`get_params`([deep])	Get parameters for this estimator.
`iterator`([stride, lag, chunk, ...])	creates an iterator to stream over the (transformed) data.
`n_chunks`(chunksize[, stride, skip])	how many chunks an iterator of this sourcde will output, starting (eg. after calling reset())
`n_frames_total`([stride, skip])	Returns total number of frames.
`number_of_trajectories`()	Returns the number of trajectories.
`output_type`()	By default transformers return single precision floats.
`partial_fit`(X)	incrementally update the covariances and mean.
`register_progress_callback`(call_back[, stage])	Registers the progress reporter.
`set_params`(**params)	Set the parameters of this estimator.
`trajectory_length`(itraj[, stride, skip])	Returns the length of trajectory of the requested index.
`trajectory_lengths`([stride, skip])	Returns the length of each trajectory.
`transform`(X)	Maps the input data through the transformer to correspondingly shaped output data array/list.
`write_to_csv`([filename, extension, ...])	write all data to csv with numpy.savetxt

Attributes

`chunksize`	chunksize defines how much data is being processed at once.
`cov`	covariance matrix of input data.
`cov_tau`	covariance matrix of time-lagged input data.
`cumvar`	Cumulative sum of the the TICA eigenvalues
`data_producer`	The data producer for this data source object (can be another data source object).
`default_chunksize`	How much data will be processed at once, in case no chunksize has been provided.
`eigenvalues`	Eigenvalues of the TICA problem (usually denoted \(\lambda\)
`eigenvectors`	Eigenvectors of the TICA problem, columnwise
`feature_TIC_correlation`	Instantaneous correlation matrix between mean-free input features and TICs
`filenames`	Property which returns a list of filenames the data is originally from.
`in_memory`	are results stored in memory?
`is_random_accessible`	Check if self._is_random_accessible is set to true and if all the random access strategies are implemented.
`is_reader`	Property telling if this data source is a reader or not.
`lag`	lag time of correlation matrix \(C_{ au}\)
`logger`	The logger for this class instance
`mean`	mean of input features
`model`	The model estimated by this Estimator
`mu`	DEPRECATED – please use the “mean” property
`name`	The name of this instance
`ndim`
`ntraj`
`ra_itraj_cuboid`	Implementation of random access with slicing that can be up to 3-dimensional, where the first dimension corresponds to the trajectory index, the second dimension corresponds to the frames and the third dimension corresponds to the dimensions of the frames.
`ra_itraj_jagged`	Behaves like ra_itraj_cuboid just that the trajectories are not truncated and returned as a list.
`ra_itraj_linear`	Implementation of random access that takes arguments as the default random access (i.e., up to three dimensions with trajs, frames and dims, respectively), but which considers the frame indexing to be contiguous.
`ra_linear`	Implementation of random access that takes a (maximal) two-dimensional slice where the first component corresponds to the frames and the second component corresponds to the dimensions.
`show_progress`	whether to show the progress of heavy calculations on this object.
`timescales`	Implied timescales of the TICA transformation

chunksize¶: chunksize defines how much data is being processed at once.

cov¶: covariance matrix of input data.

cov_tau¶: covariance matrix of time-lagged input data.

cumvar¶

Cumulative sum of the the TICA eigenvalues

Returns:	cumvar
Return type:	1D np.array

data_producer¶: The data producer for this data source object (can be another data source object). :returns: :rtype: This data source’s data producer.

default_chunksize¶: How much data will be processed at once, in case no chunksize has been provided.

describe()¶: Get a descriptive string representation of this class.

dimension()¶: output dimension

eigenvalues¶

Eigenvalues of the TICA problem (usually denoted \(\lambda\)

Returns:	eigenvalues
Return type:	1D np.array

eigenvectors¶

Eigenvectors of the TICA problem, columnwise

Returns:	eigenvectors
Return type:	(N,M) ndarray

estimate(X, **kwargs)¶: Chunk-based parameterization of TICA. Iterates over all data and estimates the mean, covariance and time lagged covariance. Finally, the generalized eigenvalue problem is solved to determine the independent components.

feature_TIC_correlation¶

Instantaneous correlation matrix between mean-free input features and TICs

Denoting the input features as \(X_i\) and the TICs as \(\theta_j\), the instantaneous, linear correlation between them can be written as

\[\mathbf{Corr}(X_i - \mu_i, \mathbf{\theta}_j) = \frac{1}{\sigma_{X_i - \mu_i}}\sum_l \sigma_{(X_i - \mu_i)(X_l - \mu_l} \mathbf{U}_{li}\]

The matrix \(\mathbf{U}\) is the matrix containing, as column vectors, the eigenvectors of the TICA generalized eigenvalue problem .

Returns:	feature_TIC_correlation – correlation matrix between input features and TICs. There is a row for each feature and a column for each TIC.
Return type:	ndarray(n,m)

filenames¶: Property which returns a list of filenames the data is originally from. :returns: list of str :rtype: list of filenames if data is originating from a file based reader

fit(X)¶

Estimates parameters - for compatibility with sklearn.

Parameters:	X (object) – A reference to the data from which the model will be estimated
Returns:	estimator – The estimator (self) with estimated model.
Return type:	object

fit_transform(X, y=None, **fit_params)¶

Fit to data, then transform it. Fits transformer to X and y with optional parameters fit_params and returns a transformed version of X. :param X: Training set. :type X: numpy array of shape [n_samples, n_features] :param y: Target values. :type y: numpy array of shape [n_samples]

Returns:	X_new – Transformed array.
Return type:	numpy array of shape [n_samples, n_features_new]

get_output(dimensions=slice(0, None, None), stride=1, skip=0, chunk=None)¶

Maps all input data of this transformer and returns it as an array or list of arrays

Parameters:	dimensions (list-like of indexes or slice, default=all) – indices of dimensions you like to keep. stride (int, default=1) – only take every n’th frame. skip (int, default=0) – initially skip n frames of each file. chunk (int, default=None) – How many frames to process at once. If not given obtain the chunk size from the source.
Returns:	output – the mapped data, where T is the number of time steps of the input data, or if stride > 1, floor(T_in / stride). d is the output dimension of this transformer. If the input consists of a list of trajectories, Y will also be a corresponding list of trajectories
Return type:	list of ndarray(T_i, d)

get_params(deep=True)¶

Get parameters for this estimator.

Parameters:	deep (boolean, optional) – If True, will return the parameters for this estimator and contained subobjects that are estimators.
Returns:	params – Parameter names mapped to their values.
Return type:	mapping of string to any

in_memory¶: are results stored in memory?

is_random_accessible¶: Check if self._is_random_accessible is set to true and if all the random access strategies are implemented. :returns: bool :rtype: Returns True if random accessible via strategies and False otherwise.

is_reader¶: Property telling if this data source is a reader or not. :returns: bool :rtype: True if this data source is a reader and False otherwise

iterator(stride=1, lag=0, chunk=None, return_trajindex=True, cols=None, skip=0)¶

creates an iterator to stream over the (transformed) data.

If your data is too large to fit into memory and you want to incrementally compute some quantities on it, you can create an iterator on a reader or transformer (eg. TICA) to avoid memory overflows.

Parameters:	stride (int, default=1) – Take only every stride’th frame. lag (int, default=0) – how many frame to omit for each file. chunk (int, default=None) – How many frames to process at once. If not given obtain the chunk size from the source. return_trajindex (boolean, default=True) – a chunk of data if return_trajindex is False, otherwise a tuple of (trajindex, data). cols (array like, default=None) – return only the given columns. skip (int, default=0) – skip ‘n’ first frames of each trajectory.
Returns:	iter – a implementation of a DataSourceIterator to stream over the data
Return type:	instance of DataSourceIterator

Examples

>>> from pyemma.coordinates import source; import numpy as np
>>> data = [np.arange(3), np.arange(4, 7)]
>>> reader = source(data)
>>> iterator = reader.iterator(chunk=1)
>>> for array_index, chunk in iterator:
...     print(array_index, chunk)
0 [[0]]
0 [[1]]
0 [[2]]
1 [[4]]
1 [[5]]
1 [[6]]

lag¶: lag time of correlation matrix \(C_{ au}\)

logger¶: The logger for this class instance

mean¶: mean of input features

model¶: The model estimated by this Estimator

mu¶: DEPRECATED – please use the “mean” property

n_chunks(chunksize, stride=1, skip=0)¶

how many chunks an iterator of this sourcde will output, starting (eg. after calling reset())

Parameters:	chunksize – stride – skip –

n_frames_total(stride=1, skip=0)¶

Returns total number of frames.

Parameters:	stride (int) – return value is the number of frames in trajectories when running through them with a step size of stride. skip (int, default=0) – skip the first initial n frames per trajectory.
Returns:	n_frames_total – total number of frames.
Return type:	int

name¶: The name of this instance

ndim¶

ntraj¶

number_of_trajectories()¶

Returns the number of trajectories.

Parameters:	stride (None (default) or np.ndarray) –
Returns:	int
Return type:	number of trajectories

output_type()¶: By default transformers return single precision floats.

partial_fit(X)¶

incrementally update the covariances and mean.

Parameters:	X (array, list of arrays, PyEMMA reader) – input data.

Notes

The projection matrix is first being calculated upon its first access.

ra_itraj_cuboid¶

Implementation of random access with slicing that can be up to 3-dimensional, where the first dimension corresponds to the trajectory index, the second dimension corresponds to the frames and the third dimension corresponds to the dimensions of the frames.

The with the frame slice selected frames will be loaded from each in the trajectory-slice selected trajectories and then sliced with the dimension slice. For example: The data consists out of three trajectories with length 10, 20, 10, respectively. The slice data[:, :15, :3] returns a 3D array of shape (3, 10, 3), where the first component corresponds to the three trajectories, the second component corresponds to 10 frames (note that the last 5 frames are being truncated as the other two trajectories only have 10 frames) and the third component corresponds to the selected first three dimensions.

Returns:	Returns an object that allows access by slices in the described manner.

ra_itraj_jagged¶

Behaves like ra_itraj_cuboid just that the trajectories are not truncated and returned as a list.

Returns:	Returns an object that allows access by slices in the described manner.

ra_itraj_linear¶

Implementation of random access that takes arguments as the default random access (i.e., up to three dimensions with trajs, frames and dims, respectively), but which considers the frame indexing to be contiguous. Therefore, it returns a simple 2D array.

Returns:	A 2D array of the sliced data containing [frames, dims].

ra_linear¶

Implementation of random access that takes a (maximal) two-dimensional slice where the first component corresponds to the frames and the second component corresponds to the dimensions. Here it is assumed that the frame indexing is contiguous, i.e., the first frame of the second trajectory has the index of the last frame of the first trajectory plus one.

Returns:	Returns an object that allows access by slices in the described manner.

register_progress_callback(call_back, stage=0)¶

Registers the progress reporter.

Parameters:	call_back (function) – This function will be called with the following arguments: stage (int) instance of pyemma.utils.progressbar.ProgressBar optional args and named keywords (kw), for future changes stage* (int, optional, default=0) – The stage you want the given call back function to be fired.

set_params(**params)¶: Set the parameters of this estimator. The method works on simple estimators as well as on nested objects (such as pipelines). The former have parameters of the form <component>__<parameter> so that it’s possible to update each component of a nested object. :returns: :rtype: self

show_progress¶: whether to show the progress of heavy calculations on this object.

timescales¶

Implied timescales of the TICA transformation

For each \(i\)-th eigenvalue, this returns

\[t_i = -\frac{\tau}{\log(|\lambda_i|)}\]

where \(\tau\) is the lag of the TICA object and \(\lambda_i\) is the i-th eigenvalue of the TICA object.

Returns:	timescales – numpy array with the implied timescales. In principle, one should expect as many timescales as input coordinates were available. However, less eigenvalues will be returned if the TICA matrices were not full rank or `var_cutoff` was parsed
Return type:	1D np.array

trajectory_length(itraj, stride=1, skip=None)¶

Returns the length of trajectory of the requested index.

Parameters:	itraj (int) – trajectory index stride (int) – return value is the number of frames in the trajectory when running through it with a step size of stride. skip (int or None) – skip n frames.
Returns:	int
Return type:	length of trajectory

trajectory_lengths(stride=1, skip=0)¶

Returns the length of each trajectory.

Parameters:	stride (int) – return value is the number of frames of the trajectories when running through them with a step size of stride. skip (int) – skip parameter
Returns:	array(dtype=int)
Return type:	containing length of each trajectory

transform(X)¶

Maps the input data through the transformer to correspondingly shaped output data array/list.

Parameters:	X (ndarray(T, n) or list of ndarray(T_i, n)) – The input data, where T is the number of time steps and n is the number of dimensions. If a list is provided, the number of time steps is allowed to vary, but the number of dimensions are required to be to be consistent.
Returns:	Y – The mapped data, where T is the number of time steps of the input data and d is the output dimension of this transformer. If called with a list of trajectories, Y will also be a corresponding list of trajectories
Return type:	ndarray(T, d) or list of ndarray(T_i, d)

write_to_csv(filename=None, extension='.dat', overwrite=False, stride=1, chunksize=100, **kw)¶

write all data to csv with numpy.savetxt

Parameters:

filename (str, optional) –
filename string, which may contain placeholders {itraj} and {stride}:
- itraj will be replaced by trajetory index
- stride is stride argument of this method
If filename is not given, it is being tried to obtain the filenames from the data source of this iterator.
extension (str, optional, default='.dat') – filename extension of created files
overwrite (bool, optional, default=False) – shall existing files be overwritten? If a file exists, this method will raise.
stride (int) – omit every n’th frame
chunksize (int) – how many frames to process at once
kw (dict) – named arguments passed into numpy.savetxt (header, seperator etc.)

Example

Assume you want to save features calculated by some FeatureReader to ASCII:

>>> import numpy as np, pyemma
>>> import os
>>> from pyemma.util.files import TemporaryDirectory
>>> from pyemma.util.contexts import settings
>>> data = [np.random.random((10,3))] * 3
>>> reader = pyemma.coordinates.source(data)
>>> filename = "distances_{itraj}.dat"
>>> with TemporaryDirectory() as td, settings(show_progress_bars=False):
...    out = os.path.join(td, filename)
...    reader.write_to_csv(out, header='', delimiter=';')
...    print(sorted(os.listdir(td)))
['distances_0.dat', 'distances_1.dat', 'distances_2.dat']

References

[1]	(1, 2, 3) Perez-Hernandez G, F Paul, T Giorgino, G De Fabritiis and F Noe. 2013. Identification of slow molecular order parameters for Markov model construction J. Chem. Phys. 139, 015102. doi:10.1063/1.4811489

[2]	L. Molgedey and H. G. Schuster. 1994. Separation of a mixture of independent signals using time delayed correlations Phys. Rev. Lett. 72, 3634.

[3]	Schwantes C, V S Pande. 2013. Improvements in Markov State Model Construction Reveal Many Non-Native Interactions in the Folding of NTL9 J. Chem. Theory. Comput. 9, 2000-2009. doi:10.1021/ct300878a

[4]	Noe, F. and Clementi, C. 2015. Kinetic distance and kinetic maps from molecular dynamics simulation. J. Chem. Theory. Comput. doi:10.1021/acs.jctc.5b00553

[5]	Noe, F., Banisch, R., Clementi, C. 2016. Commute maps: separating slowly-mixing molecular configurations for kinetic modeling. J. Chem. Theory. Comput. doi:10.1021/acs.jctc.6b00762

[6]	(1, 2) Wu, H., Nueske, F., Paul, F., Klus, S., Koltai, P., and Noe, F. 2016. Bias reduced variational approximation of molecular kinetics from short off-equilibrium simulations. J. Chem. Phys. (submitted), https://arxiv.org/abs/1610.06773.

[7]	Chan, T. F., Golub G. H., LeVeque R. J. 1979. Updating formulae and pairwiese algorithms for computing sample variances. Technical Report STAN-CS-79-773, Department of Computer Science, Stanford University.