pyemma.coordinates.pca¶

pyemma.coordinates.pca(data=None, dim=-1, var_cutoff=0.95, stride=1, mean=None, skip=0, chunksize=None, **kwargs)¶

Principal Component Analysis (PCA).

PCA is a linear transformation method that finds coordinates of maximal variance. A linear projection onto the principal components thus makes a minimal error in terms of variation in the data. Note, however, that this method is not optimal for Markov model construction because for that purpose the main objective is to preserve the slow processes which can sometimes be associated with small variance.

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

Parameters

data (ndarray (T, d) or list of ndarray (T_i, d) or a reader created by) – source function data array or list of data arrays. T or T_i are the number of time steps in a trajectory. When data is given, the PCA is immediately parametrized by estimating the covariance matrix and computing its eigenvectors.
dim (int, optional, default -1) – the number of dimensions (principal 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 variance 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
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.
mean (ndarray, optional, default None) – Optionally pass pre-calculated means to avoid their re-computation. The shape has to match the input dimension.
skip (int, default=0) – skip the first initial n frames per trajectory.
chunksize (int, default=None) – Number of data frames to process at once. Choose a higher value here, to optimize thread usage and gain processing speed. If None is passed, use the default value of the underlying reader/data source. Choose zero to disable chunking at all.

Returns

pca – Object for Principle component analysis (PCA) analysis. It contains PCA eigenvalues and eigenvectors, and the projection of input data to the dominant PCA

Return type

a PCA transformation object

Notes

Given a sequence of multivariate data \(X_t\), computes the mean-free covariance matrix.

\[C = (X - \mu)^T (X - \mu)\]

and solves the eigenvalue problem

\[C r_i = \sigma_i r_i,\]

where \(r_i\) are the principal components and \(\sigma_i\) are their respective variances.

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

See Wiki page for more theory and references. for more theory and references.

See also

PCA : pca object

tica : for time-lagged independent component analysis

class pyemma.coordinates.transform.pca.PCA(dim=-1, var_cutoff=0.95, mean=None, stride=1, skip=0)¶

Principal component analysis.

Methods

`describe`()	Get a descriptive string representation of this class.
`dimension`()	output dimension
`estimate`(X, **kwargs)	Estimates the model given the data X
`fit`(X[, y])	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.
`load`(file_name[, model_name])	Loads a previously saved PyEMMA object from disk.
`n_chunks`(chunksize[, stride, skip])	how many chunks an iterator of this sourcde will output, starting (eg.
`n_frames_total`([stride, skip])	Returns total number of frames.
`number_of_trajectories`([stride])	Returns the number of trajectories.
`output_type`()	By default transformers return single precision floats.
`partial_fit`(X)
`save`(file_name[, model_name, overwrite, …])	saves the current state of this object to given file and name.
`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
`write_to_hdf5`(filename[, group, …])	writes all data of this Iterable to a given HDF5 file.

Attributes

`chunksize`	chunksize defines how much data is being processed at once.
`cumvar`
`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`
`eigenvectors`
`feature_PC_correlation`	Instantaneous correlation matrix between input features and PCs
`filenames`	list of file names the data is originally being read 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.
`logger`	The logger for this class instance
`mean`
`model`	The model estimated by this Estimator
`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.

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

cumvar¶

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.

Notes

This variable respects your setting for maximum memory in pyemma.config.default_chunksize

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

dimension()¶: output dimension

eigenvalues¶

eigenvectors¶

estimate(X, **kwargs)¶

Estimates the model given the data X

Parameters

X (object) – A reference to the data from which the model will be estimated
params (dict) – New estimation parameter values. The parameters must that have been announced in the __init__ method of this estimator. The present settings will overwrite the settings of parameters given in the __init__ method, i.e. the parameter values after this call will be those that have been used for this estimation. Use this option if only one or a few parameters change with respect to the __init__ settings for this run, and if you don’t need to remember the original settings of these changed parameters.

Returns

estimator – The estimated estimator with the model being available.

Return type

object

feature_PC_correlation¶

Instantaneous correlation matrix between input features and PCs

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

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

The matrix \(\mathbf{U}\) is the matrix containing, as column vectors, the eigenvectors of the input-feature covariance-maxtrix.

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

filenames¶

list of file names the data is originally being read from.

Returns: names – list of file names at the beginning of the input chain.
Return type: list of str

fit(X, y=None)¶

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

classmethod load(file_name, model_name='default')¶

Loads a previously saved PyEMMA object from disk.

Parameters

file_name (str or file like object (has to provide read method)) – The file like object tried to be read for a serialized object.
model_name (str, default='default') – if multiple models are contained in the file, these can be accessed by their name. Use pyemma.list_models() to get a representation of all stored models.

Returns

obj

Return type

the de-serialized object

logger¶: The logger for this class instance

mean¶

model¶: The model estimated by this Estimator

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(stride=1)¶

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

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.

save(file_name, model_name='default', overwrite=False, save_streaming_chain=False)¶

saves the current state of this object to given file and name.

Parameters

file_name (str) – path to desired output file
model_name (str, default='default') – creates a group named ‘model_name’ in the given file, which will contain all of the data. If the name already exists, and overwrite is False (default) will raise a RuntimeError.
overwrite (bool, default=False) – Should overwrite existing model names?
save_streaming_chain (boolean, default=False) – if True, the data_producer(s) of this object will also be saved in the given file.

Examples

>>> import pyemma, numpy as np
>>> from pyemma.util.contexts import named_temporary_file
>>> m = pyemma.msm.MSM(P=np.array([[0.1, 0.9], [0.9, 0.1]]))

>>> with named_temporary_file() as file: # doctest: +SKIP
...    m.save(file, 'simple') # doctest: +SKIP
...    inst_restored = pyemma.load(file, 'simple') # doctest: +SKIP
>>> np.testing.assert_equal(m.P, inst_restored.P) # doctest: +SKIP

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

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

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=None, **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, default=None) – how many frames to process at once
kw (dict, optional) – 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']

write_to_hdf5(filename, group='/', data_set_prefix='', overwrite=False, stride=1, chunksize=None, h5_opt=None)¶

writes all data of this Iterable to a given HDF5 file. This is equivalent of writing the result of func:pyemma.coordinates.data._base.DataSource.get_output to a file.

Parameters

filename (str) – file name of output HDF5 file
group (str, default='/') – write all trajectories to this HDF5 group. The group name may not already exist in the file.
data_set_prefix (str, default=None) – data set name prefix, will postfixed with the index of the trajectory.
overwrite (bool, default=False) – if group and data sets already exist, shall we overwrite data?
stride (int, default=1) – stride argument to iterator
chunksize (int, default=None) – how many frames to process at once
h5_opt (dict) – optional parameters for h5py.create_dataset

Notes

You can pass the following via h5_opt to enable compression/filters/shuffling etc:

chunks: (Tuple) Chunk shape, or True to enable auto-chunking.
maxshape: (Tuple) Make the dataset resizable up to this shape. Use None for axes you want to be unlimited.
compression: (String or int) Compression strategy. Legal values are ‘gzip’, ‘szip’, ‘lzf’. If an integer in range(10), this indicates gzip compression level. Otherwise, an integer indicates the number of a dynamically loaded compression filter.
compression_opts: Compression settings. This is an integer for gzip, 2-tuple for szip, etc. If specifying a dynamically loaded compression filter number, this must be a tuple of values.
scaleoffset: (Integer) Enable scale/offset filter for (usually) lossy compression of integer or floating-point data. For integer data, the value of scaleoffset is the number of bits to retain (pass 0 to let HDF5 determine the minimum number of bits necessary for lossless compression). For floating point data, scaleoffset is the number of digits after the decimal place to retain; stored values thus have absolute error less than 0.5*10**(-scaleoffset).
shuffle: (T/F) Enable shuffle filter. Only effective in combination with chunks.
fletcher32: (T/F) Enable fletcher32 error detection. Not permitted in conjunction with the scale/offset filter.
fillvalue: (Scalar) Use this value for uninitialized parts of the dataset.
track_times: (T/F) Enable dataset creation timestamps.

References

1: Pearson, K. 1901 On Lines and Planes of Closest Fit to Systems of Points in Space Phil. Mag. 2, 559–572
2: Hotelling, H. 1933. Analysis of a complex of statistical variables into principal components. J. Edu. Psych. 24, 417-441 and 498-520.