From 8ce6291a4388d8a99994fbbcb70d489edd50f885 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 22 Jun 2016 17:26:29 -0700 Subject: [PATCH 01/39] Initial PCA Work --- package/MDAnalysis/analysis/pca.py | 105 +++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 package/MDAnalysis/analysis/pca.py diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py new file mode 100644 index 00000000000..b3de026a19c --- /dev/null +++ b/package/MDAnalysis/analysis/pca.py @@ -0,0 +1,105 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- http://www.MDAnalysis.org +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein +# and contributors (see AUTHORS for the full list) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + +""" +Diffusion map --- :mod:`MDAnalysis.analysis.pca` +===================================================================== + +:Authors: John Detlefs +:Year: 2016 +:Copyright: GNU Public License v3 + +TODO, write documentation + +""" +from six.moves import range +import logging + +import numpy as np +from MDAnalysis import Universe + +from .base import AnalysisBase + +logger = logging.getLogger("MDAnalysis.analysis.pca") + + +class PCA(object): + """Principal component analysis on an MD trajectory + + Attributes + ---------- + components: array, [number of frames, number of atoms] + The principal components of the feature space, + representing the directions of maximum variance in the data. + + explained_variance_ratio : array, [number of frames] + Percentage of variance explained by each of the selected components. + If a subset of components is not chosen then all components are stored + and the sum of explained variances is equal to 1.0 + + Methods + ------- + + """ + + def __init__(self, u, select='All', start=None, stop=None,step=None, + **kwargs): + """ + Parameters + ---------- + u : MDAnalysis Universe + The universe containing the trajectory frames for Principal + Component Analysis. + start : int, optional + First frame of trajectory to analyse, Default: 0 + stop : int, optional + Last frame of trajectory to analyse, Default: -1 + step : int, optional + Step between frames to analyse, Default: 1 + """ + self.setup_frames(start, stop, step) + + + def transform(self, traj=None): + """Apply the dimensionality reduction on a trajectory + + Parameters + ---------- + traj : MDAnalysis Trajectory + Trajectory for PCA transformation + + Returns + ------- + pca_space : array, shape (number of atoms, number of components) + """ + # apply Transform + # + pass + + def inverse_tranform(self, pca_space): + """ Transform PCA-transformed data back to original configuration space. + + Parameters + ---------- + pca_space : array + The space corresponding to the trajectory coordinates after the PCA + transformation. + + Returns + ------- + original_traj : array, shape (number of frames, number of atoms) + """ + pass From 74941884722270a9077a4eda364647aae75426bb Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Mon, 4 Jul 2016 10:16:30 -0700 Subject: [PATCH 02/39] Fleshed out API --- package/MDAnalysis/analysis/pca.py | 38 ++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index b3de026a19c..55106a43c1c 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -35,27 +35,33 @@ logger = logging.getLogger("MDAnalysis.analysis.pca") - class PCA(object): """Principal component analysis on an MD trajectory Attributes ---------- - components: array, [number of frames, number of atoms] + components: array, (n_components, n_atoms) The principal components of the feature space, representing the directions of maximum variance in the data. - explained_variance_ratio : array, [number of frames] + explained_variance_ratio : array, (n_components, ) Percentage of variance explained by each of the selected components. If a subset of components is not chosen then all components are stored and the sum of explained variances is equal to 1.0 + + + Methods ------- + fit(traj) + transform(traj, n_components) + + inverse_tranform(pc_space) """ - def __init__(self, u, select='All', start=None, stop=None,step=None, + def __init__(self, u, n_components=None **kwargs): """ Parameters @@ -63,17 +69,25 @@ def __init__(self, u, select='All', start=None, stop=None,step=None, u : MDAnalysis Universe The universe containing the trajectory frames for Principal Component Analysis. - start : int, optional - First frame of trajectory to analyse, Default: 0 - stop : int, optional - Last frame of trajectory to analyse, Default: -1 - step : int, optional - Step between frames to analyse, Default: 1 """ - self.setup_frames(start, stop, step) + self.u = u + self._calculated = False + + def fit(self, traj=None, n_components = None, start=None, step=None, + stop=None): + """ Use a subset of the trajectory to generate principal components """ + if traj is None: + traj = self.u.trajectory + elif isinstance(traj, u.trajectory): + traj = traj + else: + raise AttributeError("Trajectory can not be found.") + + self._setup_frames(traj, start, stop, step) + - def transform(self, traj=None): + def transform(self, traj=None, n_components=None): """Apply the dimensionality reduction on a trajectory Parameters From dbbef39bc0cbd30bbef82bbd29d0a51f356d9f05 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 6 Jul 2016 11:17:43 -0700 Subject: [PATCH 03/39] Work started on actual details --- package/MDAnalysis/analysis/pca.py | 60 ++++++++++++++++++++---------- 1 file changed, 40 insertions(+), 20 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 55106a43c1c..a2df5d58694 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -32,26 +32,25 @@ from MDAnalysis import Universe from .base import AnalysisBase +from .rms import process_selection logger = logging.getLogger("MDAnalysis.analysis.pca") -class PCA(object): +class PCA(AnalysisBase): """Principal component analysis on an MD trajectory Attributes ---------- - components: array, (n_components, n_atoms) + p_components: array, (n_components, n_atoms) The principal components of the feature space, representing the directions of maximum variance in the data. - explained_variance_ratio : array, (n_components, ) + cumulated_variance : array, (n_components, ) Percentage of variance explained by each of the selected components. If a subset of components is not chosen then all components are stored and the sum of explained variances is equal to 1.0 - - Methods ------- fit(traj) @@ -61,8 +60,7 @@ class PCA(object): inverse_tranform(pc_space) """ - def __init__(self, u, n_components=None - **kwargs): + def __init__(self, u): """ Parameters ---------- @@ -71,21 +69,41 @@ def __init__(self, u, n_components=None Component Analysis. """ self.u = u + self._calculated = False - def fit(self, traj=None, n_components = None, start=None, step=None, - stop=None): + + def fit(self, select='All', n_components=-1, start=None, stop=None, + step=None): """ Use a subset of the trajectory to generate principal components """ - if traj is None: - traj = self.u.trajectory - elif isinstance(traj, u.trajectory): - traj = traj - else: - raise AttributeError("Trajectory can not be found.") + self._setup_frames(self.u.trajectory, start, stop, step) + self.atoms = self.u.select_atoms(select) + self._n_atoms = self.atoms.n_atoms + traj = self.u.trajectory + self.run() + self._calculated = True + return self.cumulated_variance, self.p_components + - self._setup_frames(traj, start, stop, step) + def _prepare(self): + self._xyz = np.zeros((self.nframes, self._n_atoms)) + def _single_frame(self): + self._xyz[self._frame_index] = self.atoms.positions.copy() + + + def _conclude(self): + self._xyz.reshape(self.nframes, self._n_atoms * 3, order='F') + x = self._xyz - self._xyz.mean(0) + cov = np.cov(x, rowvar = 0) + e_vals, e_vects = np.linalg.eig(cov) + sort_idx = np.argsort(e_vals)[::-1] + variance = e_vals[sort_idx] + p_components = e_vects[sort_idx] + self.cumulated_variance = np.cumsum(variance) + self.p_components = p_components[:n_components] + def transform(self, traj=None, n_components=None): """Apply the dimensionality reduction on a trajectory @@ -99,9 +117,10 @@ def transform(self, traj=None, n_components=None): ------- pca_space : array, shape (number of atoms, number of components) """ - # apply Transform - # - pass + if traj is None: + traj = self._original_traj + return np.dot(traj, self.p_components) + def inverse_tranform(self, pca_space): """ Transform PCA-transformed data back to original configuration space. @@ -116,4 +135,5 @@ def inverse_tranform(self, pca_space): ------- original_traj : array, shape (number of frames, number of atoms) """ - pass + inv = np.linalg.inv(self.p_components) + return np.dot(pca_space, inv) From ad250484806e00ebd63a1dfa5446e619bcff5591 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 6 Jul 2016 14:06:24 -0700 Subject: [PATCH 04/39] More work done --- package/MDAnalysis/analysis/pca.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index a2df5d58694..475f0b0b401 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -77,6 +77,7 @@ def fit(self, select='All', n_components=-1, start=None, stop=None, step=None): """ Use a subset of the trajectory to generate principal components """ self._setup_frames(self.u.trajectory, start, stop, step) + self.n_components = n_components self.atoms = self.u.select_atoms(select) self._n_atoms = self.atoms.n_atoms traj = self.u.trajectory @@ -86,7 +87,7 @@ def fit(self, select='All', n_components=-1, start=None, stop=None, def _prepare(self): - self._xyz = np.zeros((self.nframes, self._n_atoms)) + self._xyz = np.zeros((self.nframes, self._n_atoms, 3)) def _single_frame(self): @@ -94,15 +95,16 @@ def _single_frame(self): def _conclude(self): - self._xyz.reshape(self.nframes, self._n_atoms * 3, order='F') + self._xyz = self._xyz.reshape(self.nframes, self._n_atoms * 3, + order='F') x = self._xyz - self._xyz.mean(0) cov = np.cov(x, rowvar = 0) e_vals, e_vects = np.linalg.eig(cov) sort_idx = np.argsort(e_vals)[::-1] - variance = e_vals[sort_idx] + self.variance = e_vals[sort_idx] p_components = e_vects[sort_idx] - self.cumulated_variance = np.cumsum(variance) - self.p_components = p_components[:n_components] + self.cumulated_variance = np.cumsum(variance) / np.sum(variance) + self.p_components = p_components[:self.n_components] def transform(self, traj=None, n_components=None): From c732a922683bf3455ca62086fc56591f9b987b33 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 6 Jul 2016 15:52:50 -0700 Subject: [PATCH 05/39] More work done, Docs update, added variance attribute --- package/MDAnalysis/analysis/pca.py | 47 +++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 13 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 475f0b0b401..050b496c205 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -44,20 +44,24 @@ class PCA(AnalysisBase): p_components: array, (n_components, n_atoms) The principal components of the feature space, representing the directions of maximum variance in the data. - - cumulated_variance : array, (n_components, ) + variance : array (n_components, ) + The raw variance explained by each eigenvector of the covariance + matrix. + explained_variance : array, (n_components, ) Percentage of variance explained by each of the selected components. If a subset of components is not chosen then all components are stored - and the sum of explained variances is equal to 1.0 - + and the sum of explained variances is equal to 1.0. Methods ------- - fit(traj) - - transform(traj, n_components) - - inverse_tranform(pc_space) + fit(traj=None, select='All', start=None, stop=None, step=None) + Find the principal components of selected atoms in a subset of the + trajectory. + transform(traj=None, n_components=None) + Take the original trajectory and project it onto the principal + components. + inverse_tranform(pca_space) + Take a pca_space and map it back onto the trajectory used to create it. """ def __init__(self, u): @@ -69,18 +73,36 @@ def __init__(self, u): Component Analysis. """ self.u = u - + self._original_traj = self.u.trajectory self._calculated = False def fit(self, select='All', n_components=-1, start=None, stop=None, step=None): - """ Use a subset of the trajectory to generate principal components """ + """ Use a subset of the trajectory to generate principal components. + + Parameters + ---------- + select : string, optional + A valid select statement for picking a subset of atoms from an + MDAnalysis Trajectory. + n_components : int, optional + The number of principal components to be saved, default saves + all principal components, Default: -1 + start : int, optional + First frame of trajectory to use for generation + of covariance matrix, Default: 0 + stop : int, optional + Last frame of trajectory to use for generation + of covariance matrix, Default: -1 + step : int, optional + Step between frames of trajectory to use for generation + of covariance matrix, Default: 1 + """ self._setup_frames(self.u.trajectory, start, stop, step) self.n_components = n_components self.atoms = self.u.select_atoms(select) self._n_atoms = self.atoms.n_atoms - traj = self.u.trajectory self.run() self._calculated = True return self.cumulated_variance, self.p_components @@ -114,7 +136,6 @@ def transform(self, traj=None, n_components=None): ---------- traj : MDAnalysis Trajectory Trajectory for PCA transformation - Returns ------- pca_space : array, shape (number of atoms, number of components) From 9fdc3906cfe5123e457d4f1d30202c426ec3c425 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 6 Jul 2016 21:08:33 -0700 Subject: [PATCH 06/39] more docs --- package/MDAnalysis/analysis/pca.py | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 050b496c205..08353087e6f 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -22,7 +22,19 @@ :Year: 2016 :Copyright: GNU Public License v3 -TODO, write documentation +This module contains the linear dimension reduction method Principal +Component Analysis. This module constructs a covariance matrix wherein each +element of the matrix is denoted by (i,j) row-column coordinates. The (i,j) +coordinate reflects the influence of the of the ith frame on the jth frame +of the trajectory. The Principal Components are the eigenvectors of this matrix. + +For each eigenvector, its eigenvalue reflects the variance that the eigenvector +explains. This value is made into a ratio stored in `explained_variance`, which +provides divides accumulated variance of the nth eigenvector and the +n-1 eigenvectors preceding by the total variance in the data. + +From here, we can project a trajectory onto these principal components, + """ from six.moves import range @@ -32,7 +44,6 @@ from MDAnalysis import Universe from .base import AnalysisBase -from .rms import process_selection logger = logging.getLogger("MDAnalysis.analysis.pca") @@ -62,6 +73,7 @@ class PCA(AnalysisBase): components. inverse_tranform(pca_space) Take a pca_space and map it back onto the trajectory used to create it. + """ def __init__(self, u): @@ -125,7 +137,8 @@ def _conclude(self): sort_idx = np.argsort(e_vals)[::-1] self.variance = e_vals[sort_idx] p_components = e_vects[sort_idx] - self.cumulated_variance = np.cumsum(variance) / np.sum(variance) + self.cumulated_variance = (np.cumsum(self.variance) / + np.sum(self.variance)) self.p_components = p_components[:self.n_components] @@ -152,7 +165,8 @@ def inverse_tranform(self, pca_space): ---------- pca_space : array The space corresponding to the trajectory coordinates after the PCA - transformation. + transformation. Assumes current principal components were the + components projected onto to create the pca_space. Returns ------- From c673951d261cb290f7e7c1a2b728539068387c2c Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 7 Jul 2016 11:57:28 -0700 Subject: [PATCH 07/39] Mostly finished with docs --- package/MDAnalysis/analysis/pca.py | 66 ++++++++++++++++++++++++------ 1 file changed, 53 insertions(+), 13 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 08353087e6f..aa456c4d66d 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -15,7 +15,7 @@ # """ -Diffusion map --- :mod:`MDAnalysis.analysis.pca` +Principal Component Analysis (PCA) --- :mod:`MDAnalysis.analysis.pca` ===================================================================== :Authors: John Detlefs @@ -25,17 +25,55 @@ This module contains the linear dimension reduction method Principal Component Analysis. This module constructs a covariance matrix wherein each element of the matrix is denoted by (i,j) row-column coordinates. The (i,j) -coordinate reflects the influence of the of the ith frame on the jth frame -of the trajectory. The Principal Components are the eigenvectors of this matrix. +coordinate reflects the influence of the of the ith frame's coordinates on the +jth frame's coordinates of a given trajectory. The eponymous components are the +eigenvectors of this matrix. For each eigenvector, its eigenvalue reflects the variance that the eigenvector -explains. This value is made into a ratio stored in `explained_variance`, which -provides divides accumulated variance of the nth eigenvector and the -n-1 eigenvectors preceding by the total variance in the data. - -From here, we can project a trajectory onto these principal components, - - +explains. This value is made into a ratio. Stored in +:attribute:`explained_variance`, this ratio divides the accumulated variance +of the nth eigenvector and the n-1 eigenvectors preceding the eigenvector by +the total variance in the data. For most data, :attribute:`explained_variance` +will be approximately equal to one for some n that is significantly smaller +than the total number of components, these are the components of interest given +by Principal Component Analysis. + +From here, we can project a trajectory onto these principal components and +attempt to retrieve some structure from our high dimensional data. We have +provided a [notebook](# TODO edit max's notebook to use the new module) +containing a thorough demonstration of Principal Component Analysis. + +For a basic introduction to the module, the :ref:`PCA-tutorial` shows how +to perform Principal Component Analysis. + + +.. _PCA-tutorial: +The example uses files provided as part of the MDAnalysis test suite +(in the variables :data:`~MDAnalysis.tests.datafiles.PSF` and +:data:`~MDAnalysis.tests.datafiles.DCD`). This tutorial shows how to use the +PCA class. + + +First load all modules and test data :: + >>> import MDAnalysis as mda + >>> import numpy as np + >>> import MDAnalysis.analysis.pca as pca + >>> from MDAnalysis.tests.datafiles import PSF, DCD + +Given a universe containing trajectory data we can perform PCA using +:class:`PCA`:: and retrieve the principal components. + >>> u = mda.Universe(PSF,DCD) + >>> PSF_pca = pca.PCA(u) + >>> cumulated_variance, principal_components = PSF_pca.fit() + +Inspect the components to determine the principal components you would like +to retain. The choice is arbitrary, but I will stop when 95 percent of the +variance is explained by the components. + >>> n_pcs = next(x[0] for x in enumerate(cumulated_variance) if x[1] > 0.95) + >>> pca_space = PSF_pca.transform(n_components=n_pcs) + +From here, inspection of the pca_space and conclusions to be drawn from the +data are left to the user. """ from six.moves import range import logging @@ -147,15 +185,17 @@ def transform(self, traj=None, n_components=None): Parameters ---------- - traj : MDAnalysis Trajectory - Trajectory for PCA transformation + traj : MDAnalysis Universe + Universe containing trajectory for PCA transformation Returns ------- pca_space : array, shape (number of atoms, number of components) """ if traj is None: traj = self._original_traj - return np.dot(traj, self.p_components) + + ca_xyz = np.array([ca.positions.copy() for ts in traj]) + return np.dot(traj, self.p_components[:n_components]) def inverse_tranform(self, pca_space): From 7501982d445aa10dab54f7ff3626c3e4f41e6ef8 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 7 Jul 2016 13:09:29 -0700 Subject: [PATCH 08/39] Transform function done --- package/MDAnalysis/analysis/pca.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index aa456c4d66d..b428fb0af73 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -111,7 +111,6 @@ class PCA(AnalysisBase): components. inverse_tranform(pca_space) Take a pca_space and map it back onto the trajectory used to create it. - """ def __init__(self, u): @@ -123,7 +122,6 @@ def __init__(self, u): Component Analysis. """ self.u = u - self._original_traj = self.u.trajectory self._calculated = False @@ -180,22 +178,28 @@ def _conclude(self): self.p_components = p_components[:self.n_components] - def transform(self, traj=None, n_components=None): + def transform(self, select='All', u=None, n_components=None): """Apply the dimensionality reduction on a trajectory Parameters ---------- - traj : MDAnalysis Universe + select : string, optional + A valid select statement for picking a subset of atoms from an + MDAnalysis Universe. + u : MDAnalysis Universe Universe containing trajectory for PCA transformation Returns ------- - pca_space : array, shape (number of atoms, number of components) + pca_space : array, shape (number of atoms*3, number of components) """ - if traj is None: - traj = self._original_traj + if u is None: + u = self.u + + atoms = u.select_atoms(select) - ca_xyz = np.array([ca.positions.copy() for ts in traj]) - return np.dot(traj, self.p_components[:n_components]) + xyz = np.array([atoms.positions.copy() for ts in u.trajectory]) + xyz.reshape(u.trajectory.nframes, atoms.n_atoms*3, order='F') + return np.dot(xyz, self.p_components[:n_components]) def inverse_tranform(self, pca_space): @@ -210,7 +214,7 @@ def inverse_tranform(self, pca_space): Returns ------- - original_traj : array, shape (number of frames, number of atoms) + original_traj : array, shape (number of frames, number of atoms*3) """ inv = np.linalg.inv(self.p_components) return np.dot(pca_space, inv) From f59c1ba354ece6eeae1fc3e5f88fc87e907bddbe Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 8 Jul 2016 13:26:09 -0700 Subject: [PATCH 09/39] Test driven bug fixing Updated tests --- package/MDAnalysis/analysis/pca.py | 68 +++++++++++-------- .../MDAnalysisTests/analysis/test_pca.py | 48 +++++++++++++ 2 files changed, 86 insertions(+), 30 deletions(-) create mode 100644 testsuite/MDAnalysisTests/analysis/test_pca.py diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index b428fb0af73..83686a7e958 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -113,7 +113,7 @@ class PCA(AnalysisBase): Take a pca_space and map it back onto the trajectory used to create it. """ - def __init__(self, u): + def __init__(self, u, select='all'): """ Parameters ---------- @@ -121,19 +121,20 @@ def __init__(self, u): The universe containing the trajectory frames for Principal Component Analysis. """ - self.u = u + self._u = u + self._atoms = self._u.select_atoms(select) + self._n_atoms = self._atoms.n_atoms self._calculated = False - def fit(self, select='All', n_components=-1, start=None, stop=None, + def fit(self, n_components=None, start=None, stop=None, step=None): - """ Use a subset of the trajectory to generate principal components. + """ Use a subset of frames from the trajectory to generate the + principal components. Parameters ---------- - select : string, optional - A valid select statement for picking a subset of atoms from an - MDAnalysis Trajectory. + n_components : int, optional The number of principal components to be saved, default saves all principal components, Default: -1 @@ -146,60 +147,67 @@ def fit(self, select='All', n_components=-1, start=None, stop=None, step : int, optional Step between frames of trajectory to use for generation of covariance matrix, Default: 1 + + Return + ------ + cumulated_variance: array, (n_components, ) + The amount of variance explained by the nth component and the n-1 + components preceding it. + p_components: array, (n_components, n_atoms * 3) + """ - self._setup_frames(self.u.trajectory, start, stop, step) + self._setup_frames(self._u.trajectory, start, stop, step) self.n_components = n_components - self.atoms = self.u.select_atoms(select) - self._n_atoms = self.atoms.n_atoms self.run() self._calculated = True - return self.cumulated_variance, self.p_components + if n_components is None: + return self.cumulated_variance, self.p_components + else: + return (self.cumulated_variance[:n_components], + self.p_components[:n_components]) def _prepare(self): - self._xyz = np.zeros((self.nframes, self._n_atoms, 3)) + self._xyz = np.zeros((self.n_frames, self._n_atoms, 3)) def _single_frame(self): - self._xyz[self._frame_index] = self.atoms.positions.copy() + self._xyz[self._frame_index] = self._atoms.positions.copy() def _conclude(self): - self._xyz = self._xyz.reshape(self.nframes, self._n_atoms * 3, + self._xyz = self._xyz.reshape(self.n_frames, self._n_atoms * 3, order='F') x = self._xyz - self._xyz.mean(0) cov = np.cov(x, rowvar = 0) e_vals, e_vects = np.linalg.eig(cov) sort_idx = np.argsort(e_vals)[::-1] self.variance = e_vals[sort_idx] - p_components = e_vects[sort_idx] + self.p_components = e_vects[sort_idx] self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance)) - self.p_components = p_components[:self.n_components] - def transform(self, select='All', u=None, n_components=None): + def transform(self, n_components=None): """Apply the dimensionality reduction on a trajectory Parameters ---------- - select : string, optional - A valid select statement for picking a subset of atoms from an - MDAnalysis Universe. - u : MDAnalysis Universe - Universe containing trajectory for PCA transformation + n_components: int, optional + The number of components to be projected onto, default maps onto + all components. Returns ------- - pca_space : array, shape (number of atoms*3, number of components) + pca_space : array, shape (number of frames, number of components) """ - if u is None: - u = self.u - - atoms = u.select_atoms(select) - xyz = np.array([atoms.positions.copy() for ts in u.trajectory]) - xyz.reshape(u.trajectory.nframes, atoms.n_atoms*3, order='F') - return np.dot(xyz, self.p_components[:n_components]) + xyz = np.array([self._atoms.positions.copy() for ts in self._u.trajectory]) + xyz = xyz.reshape(self._u.trajectory.n_frames, + self._n_atoms*3, order='F') + if n_components is None: + return np.dot(xyz, self.p_components.T) + else: + return np.dot(xyz, self.p_components[:n_components].T) def inverse_tranform(self, pca_space): diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py new file mode 100644 index 00000000000..daf2a0cfd37 --- /dev/null +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -0,0 +1,48 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- http://www.MDAnalysis.org +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# +from __future__ import print_function +import numpy as np +import MDAnalysis +import MDAnalysis.analysis.pca as pca +from numpy.testing import (assert_almost_equal, assert_equal, + assert_array_almost_equal,raises) + + +from MDAnalysisTests.datafiles import PDB, XTC + + +class TestPCA(object): + def setUp(self): + self.u = MDAnalysis.Universe(PDB, XTC) + self.pca = pca.PCA(self.u, select='backbone and name CA') + self.cumulated_variance, self.pcs = self.pca.fit() + self.n_atoms = self.u.select_atoms('backbone and name CA').n_atoms + + def test_cum_var(self): + assert_almost_equal(self.cumulated_variance[-1], 1) + # makes no sense to test values here, no physical meaning + + def test_pcs(self): + assert_equal(self.pcs.shape, (self.n_atoms*3, self.n_atoms*3)) + + def test_different_steps(self): + cum_var, pcs = self.pca.fit(step=3) + + def test_transform(self): + self.n_components = 1 + pca_space = self.pca.transform(n_components=self.n_components) + assert_equal(pca_space.shape, + (self.u.trajectory.n_frames, self.n_components)) From bdb3848934bc9b0e4bd40d3337598a7a5c9c8139 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 21 Jul 2016 16:04:43 -0700 Subject: [PATCH 10/39] update to Covariance algorithm --- package/MDAnalysis/analysis/pca.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 83686a7e958..e04a9584fbd 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -83,7 +83,7 @@ from .base import AnalysisBase -logger = logging.getLogger("MDAnalysis.analysis.pca") +logger = logging.getLogger(__name__) class PCA(AnalysisBase): """Principal component analysis on an MD trajectory @@ -126,7 +126,6 @@ def __init__(self, u, select='all'): self._n_atoms = self._atoms.n_atoms self._calculated = False - def fit(self, n_components=None, start=None, stop=None, step=None): """ Use a subset of frames from the trajectory to generate the @@ -157,6 +156,10 @@ def fit(self, n_components=None, start=None, stop=None, """ self._setup_frames(self._u.trajectory, start, stop, step) + self.start = start + self.stop = stop + self.step = step + self.n_components = n_components self.run() self._calculated = True @@ -168,26 +171,28 @@ def fit(self, n_components=None, start=None, stop=None, self.p_components[:n_components]) def _prepare(self): - self._xyz = np.zeros((self.n_frames, self._n_atoms, 3)) - + n_dim = self._n_atoms * 3 + self.cov = np.zeros((n_dim, n_dim)) + self.mean = np.zeros(n_dim,) + for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): + x = self._atoms.positions.ravel() + self.mean += x + self.mean /= self.n_frames def _single_frame(self): - self._xyz[self._frame_index] = self._atoms.positions.copy() - + x = self._atoms.positions.ravel() + x -= self.mean + self.cov += np.dot(x[:, np.newaxis], x[:, np.newaxis].T) def _conclude(self): - self._xyz = self._xyz.reshape(self.n_frames, self._n_atoms * 3, - order='F') - x = self._xyz - self._xyz.mean(0) - cov = np.cov(x, rowvar = 0) - e_vals, e_vects = np.linalg.eig(cov) + self.cov /= self.n_frames - 1 + e_vals, e_vects = np.linalg.eig(self.cov) sort_idx = np.argsort(e_vals)[::-1] self.variance = e_vals[sort_idx] self.p_components = e_vects[sort_idx] self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance)) - def transform(self, n_components=None): """Apply the dimensionality reduction on a trajectory From b945237147fcf5a47dcf3fd31f08f64251e642e8 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 21 Jul 2016 16:21:00 -0700 Subject: [PATCH 11/39] Transform update --- package/MDAnalysis/analysis/pca.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index e04a9584fbd..c2078842e06 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -193,7 +193,7 @@ def _conclude(self): self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance)) - def transform(self, n_components=None): + def transform(self, atomgroup=None, n_components=None): """Apply the dimensionality reduction on a trajectory Parameters @@ -205,8 +205,15 @@ def transform(self, n_components=None): ------- pca_space : array, shape (number of frames, number of components) """ + if atomgroup is not None: + atoms = atomgroup + else: + atoms = self._atoms + + if self._atoms != atoms: + warnings.warn('This is a transform for different atom types.') - xyz = np.array([self._atoms.positions.copy() for ts in self._u.trajectory]) + xyz = np.array([atoms.positions.copy() for ts in self._u.trajectory]) xyz = xyz.reshape(self._u.trajectory.n_frames, self._n_atoms*3, order='F') if n_components is None: From 1678d58122170bea80f0ed24da311e56e830ca2f Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sat, 23 Jul 2016 17:26:31 -0700 Subject: [PATCH 12/39] Switched to new base api --- package/MDAnalysis/analysis/pca.py | 56 ++++++++++++------------------ 1 file changed, 22 insertions(+), 34 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index c2078842e06..c3b1ba73ea6 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -113,27 +113,13 @@ class PCA(AnalysisBase): Take a pca_space and map it back onto the trajectory used to create it. """ - def __init__(self, u, select='all'): + def __init__(self, atomgroup, select='All', n_components=None, + **kwargs): """ Parameters ---------- - u : MDAnalysis Universe - The universe containing the trajectory frames for Principal - Component Analysis. - """ - self._u = u - self._atoms = self._u.select_atoms(select) - self._n_atoms = self._atoms.n_atoms - self._calculated = False - - def fit(self, n_components=None, start=None, stop=None, - step=None): - """ Use a subset of frames from the trajectory to generate the - principal components. - - Parameters - ---------- - + atomgroup: MDAnalysis atomgroup + AtomGroup to be used for fitting. n_components : int, optional The number of principal components to be saved, default saves all principal components, Default: -1 @@ -146,7 +132,21 @@ def fit(self, n_components=None, start=None, stop=None, step : int, optional Step between frames of trajectory to use for generation of covariance matrix, Default: 1 + """ + super(PCA, self).__init__(atomgroup.universe.trajectory, + **kwargs) + + self._atoms = atomgroup.select_atoms(select) + self.n_components = n_components + self._n_atoms = self._atoms.n_atoms + self._calculated = False + + def fit(self): + """ Use a subset of frames from the trajectory to generate the + principal components. + Parameters + ---------- Return ------ cumulated_variance: array, (n_components, ) @@ -155,20 +155,8 @@ def fit(self, n_components=None, start=None, stop=None, p_components: array, (n_components, n_atoms * 3) """ - self._setup_frames(self._u.trajectory, start, stop, step) - self.start = start - self.stop = stop - self.step = step - - self.n_components = n_components - self.run() - self._calculated = True - - if n_components is None: - return self.cumulated_variance, self.p_components - else: - return (self.cumulated_variance[:n_components], - self.p_components[:n_components]) + return (self.cumulated_variance[:n_components], + self.p_components[:n_components]) def _prepare(self): n_dim = self._n_atoms * 3 @@ -192,8 +180,9 @@ def _conclude(self): self.p_components = e_vects[sort_idx] self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance)) + self._calculated = True - def transform(self, atomgroup=None, n_components=None): + def transform(self, atomgroup, n_components=None): """Apply the dimensionality reduction on a trajectory Parameters @@ -221,7 +210,6 @@ def transform(self, atomgroup=None, n_components=None): else: return np.dot(xyz, self.p_components[:n_components].T) - def inverse_tranform(self, pca_space): """ Transform PCA-transformed data back to original configuration space. From 5dc007668d310227ca6b3691e61666fff9c5a510 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sat, 23 Jul 2016 18:32:07 -0700 Subject: [PATCH 13/39] Tests working after refactor --- package/MDAnalysis/analysis/pca.py | 27 +++++++++---------- .../MDAnalysisTests/analysis/test_pca.py | 10 ++++--- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index c3b1ba73ea6..0e532f32e23 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -77,6 +77,7 @@ """ from six.moves import range import logging +import warnings import numpy as np from MDAnalysis import Universe @@ -113,7 +114,7 @@ class PCA(AnalysisBase): Take a pca_space and map it back onto the trajectory used to create it. """ - def __init__(self, atomgroup, select='All', n_components=None, + def __init__(self, atomgroup, select='all', n_components=None, **kwargs): """ Parameters @@ -135,7 +136,10 @@ def __init__(self, atomgroup, select='All', n_components=None, """ super(PCA, self).__init__(atomgroup.universe.trajectory, **kwargs) - + self._u = atomgroup.universe + self.start = kwargs.get('start') + self.stop = kwargs.get('stop') + self.step = kwargs.get('step') self._atoms = atomgroup.select_atoms(select) self.n_components = n_components self._n_atoms = self._atoms.n_atoms @@ -155,8 +159,9 @@ def fit(self): p_components: array, (n_components, n_atoms * 3) """ - return (self.cumulated_variance[:n_components], - self.p_components[:n_components]) + self.run() + return (self.cumulated_variance[:self.n_components], + self.p_components[:self.n_components]) def _prepare(self): n_dim = self._n_atoms * 3 @@ -194,21 +199,15 @@ def transform(self, atomgroup, n_components=None): ------- pca_space : array, shape (number of frames, number of components) """ - if atomgroup is not None: - atoms = atomgroup - else: - atoms = self._atoms - if self._atoms != atoms: + if self._atoms != atomgroup: warnings.warn('This is a transform for different atom types.') - xyz = np.array([atoms.positions.copy() for ts in self._u.trajectory]) + xyz = np.array([atomgroup.positions.copy() for ts in self._u.trajectory]) xyz = xyz.reshape(self._u.trajectory.n_frames, self._n_atoms*3, order='F') - if n_components is None: - return np.dot(xyz, self.p_components.T) - else: - return np.dot(xyz, self.p_components[:n_components].T) + + return np.dot(xyz, self.p_components[:n_components].T) def inverse_tranform(self, pca_space): """ Transform PCA-transformed data back to original configuration space. diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index daf2a0cfd37..4c5ba06e387 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -27,7 +27,8 @@ class TestPCA(object): def setUp(self): self.u = MDAnalysis.Universe(PDB, XTC) - self.pca = pca.PCA(self.u, select='backbone and name CA') + self.pca = pca.PCA(self.u.atoms, select='backbone and name CA') + self.pca.run() self.cumulated_variance, self.pcs = self.pca.fit() self.n_atoms = self.u.select_atoms('backbone and name CA').n_atoms @@ -39,10 +40,13 @@ def test_pcs(self): assert_equal(self.pcs.shape, (self.n_atoms*3, self.n_atoms*3)) def test_different_steps(self): - cum_var, pcs = self.pca.fit(step=3) + self.pca.run() + cum_var = self.pca.cumulated_variance + pcs = self.pca.p_components def test_transform(self): self.n_components = 1 - pca_space = self.pca.transform(n_components=self.n_components) + self.ag = self.u.select_atoms('backbone and name CA') + pca_space = self.pca.transform(self.ag, n_components=self.n_components) assert_equal(pca_space.shape, (self.u.trajectory.n_frames, self.n_components)) From 9aa086d140824d1f491160f3a079b60d61790528 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Mon, 25 Jul 2016 12:25:09 -0700 Subject: [PATCH 14/39] Updated transform --- package/MDAnalysis/analysis/pca.py | 57 +++++++++++++++++------------- 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 0e532f32e23..78368b21f7d 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -56,7 +56,6 @@ First load all modules and test data :: >>> import MDAnalysis as mda - >>> import numpy as np >>> import MDAnalysis.analysis.pca as pca >>> from MDAnalysis.tests.datafiles import PSF, DCD @@ -86,6 +85,7 @@ logger = logging.getLogger(__name__) + class PCA(AnalysisBase): """Principal component analysis on an MD trajectory @@ -114,8 +114,8 @@ class PCA(AnalysisBase): Take a pca_space and map it back onto the trajectory used to create it. """ - def __init__(self, atomgroup, select='all', n_components=None, - **kwargs): + def __init__(self, atomgroup, select='all', mean_free=True, reference=None, + n_components=None, **kwargs): """ Parameters ---------- @@ -140,33 +140,21 @@ def __init__(self, atomgroup, select='all', n_components=None, self.start = kwargs.get('start') self.stop = kwargs.get('stop') self.step = kwargs.get('step') + if mean_free: + if reference is None: + self.reference = self.atomgroup + else: + self.reference = reference self._atoms = atomgroup.select_atoms(select) self.n_components = n_components self._n_atoms = self._atoms.n_atoms self._calculated = False - def fit(self): - """ Use a subset of frames from the trajectory to generate the - principal components. - - Parameters - ---------- - Return - ------ - cumulated_variance: array, (n_components, ) - The amount of variance explained by the nth component and the n-1 - components preceding it. - p_components: array, (n_components, n_atoms * 3) - - """ - self.run() - return (self.cumulated_variance[:self.n_components], - self.p_components[:self.n_components]) - def _prepare(self): n_dim = self._n_atoms * 3 self.cov = np.zeros((n_dim, n_dim)) self.mean = np.zeros(n_dim,) + for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): x = self._atoms.positions.ravel() self.mean += x @@ -187,6 +175,23 @@ def _conclude(self): np.sum(self.variance)) self._calculated = True + def fit(self): + """ Use a subset of frames from the trajectory to generate the + principal components. + + Parameters + ---------- + Return + ------ + cumulated_variance: array, (n_components, ) + The amount of variance explained by the nth component and the n-1 + components preceding it. + p_components: array, (n_components, n_atoms * 3) + """ + self.run() + return (self.cumulated_variance[:self.n_components], + self.p_components[:self.n_components]) + def transform(self, atomgroup, n_components=None): """Apply the dimensionality reduction on a trajectory @@ -203,11 +208,13 @@ def transform(self, atomgroup, n_components=None): if self._atoms != atomgroup: warnings.warn('This is a transform for different atom types.') - xyz = np.array([atomgroup.positions.copy() for ts in self._u.trajectory]) - xyz = xyz.reshape(self._u.trajectory.n_frames, - self._n_atoms*3, order='F') + dot = np.zeros((u.trajectory.n_frames, ca.n_atoms*3)) + + for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): + xyz = atomgroup.positions.ravel() + dot[i] = np.dot(xyz, self.p_components[:n_components].T) - return np.dot(xyz, self.p_components[:n_components].T) + return dot def inverse_tranform(self, pca_space): """ Transform PCA-transformed data back to original configuration space. From 740fc89c9460123b7feac672939d7fbdf6d6afe1 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Mon, 25 Jul 2016 15:17:12 -0700 Subject: [PATCH 15/39] Some demeaning fixing, work in progress --- package/MDAnalysis/analysis/pca.py | 48 ++++++++++++------------------ 1 file changed, 19 insertions(+), 29 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 78368b21f7d..46d14dcfc69 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -114,25 +114,25 @@ class PCA(AnalysisBase): Take a pca_space and map it back onto the trajectory used to create it. """ - def __init__(self, atomgroup, select='all', mean_free=True, reference=None, + def __init__(self, atomgroup, demean=True, reference=None, n_components=None, **kwargs): """ Parameters ---------- atomgroup: MDAnalysis atomgroup - AtomGroup to be used for fitting. + AtomGroup to be used for PCA. n_components : int, optional The number of principal components to be saved, default saves all principal components, Default: -1 start : int, optional First frame of trajectory to use for generation - of covariance matrix, Default: 0 + of covariance matrix, Default: None stop : int, optional Last frame of trajectory to use for generation - of covariance matrix, Default: -1 + of covariance matrix, Default: None step : int, optional Step between frames of trajectory to use for generation - of covariance matrix, Default: 1 + of covariance matrix, Default: None """ super(PCA, self).__init__(atomgroup.universe.trajectory, **kwargs) @@ -140,15 +140,22 @@ def __init__(self, atomgroup, select='all', mean_free=True, reference=None, self.start = kwargs.get('start') self.stop = kwargs.get('stop') self.step = kwargs.get('step') - if mean_free: - if reference is None: - self.reference = self.atomgroup + + if demean: + if frame is None: + self.reference = self._u else: - self.reference = reference - self._atoms = atomgroup.select_atoms(select) + # change u to traj number for fitting + self._u.trajectory[frame] + self.reference = self._u + aligned = align.AlignTraj(self.reference, self._u).run() + # need help here + self._u = mda.Universe(aligned.filename,) + atomgroup.universe = self._u + + self._atoms = atomgroup self.n_components = n_components self._n_atoms = self._atoms.n_atoms - self._calculated = False def _prepare(self): n_dim = self._n_atoms * 3 @@ -158,6 +165,7 @@ def _prepare(self): for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): x = self._atoms.positions.ravel() self.mean += x + self.mean /= self.n_frames def _single_frame(self): @@ -173,24 +181,6 @@ def _conclude(self): self.p_components = e_vects[sort_idx] self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance)) - self._calculated = True - - def fit(self): - """ Use a subset of frames from the trajectory to generate the - principal components. - - Parameters - ---------- - Return - ------ - cumulated_variance: array, (n_components, ) - The amount of variance explained by the nth component and the n-1 - components preceding it. - p_components: array, (n_components, n_atoms * 3) - """ - self.run() - return (self.cumulated_variance[:self.n_components], - self.p_components[:self.n_components]) def transform(self, atomgroup, n_components=None): """Apply the dimensionality reduction on a trajectory From bb24f2b14348804b0269f78271334fca184a043d Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 15:34:05 -0700 Subject: [PATCH 16/39] Crappy code to jump off from. Demeaning using _fit_to. --- package/MDAnalysis/analysis/pca.py | 42 ++++++++++++++---------------- 1 file changed, 19 insertions(+), 23 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 46d14dcfc69..ab75a1b19cf 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -80,6 +80,7 @@ import numpy as np from MDAnalysis import Universe +from MDAnalysis.analysis.align import _fit_to from .base import AnalysisBase @@ -137,22 +138,13 @@ def __init__(self, atomgroup, demean=True, reference=None, super(PCA, self).__init__(atomgroup.universe.trajectory, **kwargs) self._u = atomgroup.universe + self.reference = reference + # for transform function self.start = kwargs.get('start') self.stop = kwargs.get('stop') self.step = kwargs.get('step') - if demean: - if frame is None: - self.reference = self._u - else: - # change u to traj number for fitting - self._u.trajectory[frame] - self.reference = self._u - aligned = align.AlignTraj(self.reference, self._u).run() - # need help here - self._u = mda.Universe(aligned.filename,) - atomgroup.universe = self._u - + self.demean = demean self._atoms = atomgroup self.n_components = n_components self._n_atoms = self._atoms.n_atoms @@ -160,17 +152,21 @@ def __init__(self, atomgroup, demean=True, reference=None, def _prepare(self): n_dim = self._n_atoms * 3 self.cov = np.zeros((n_dim, n_dim)) - self.mean = np.zeros(n_dim,) - - for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): - x = self._atoms.positions.ravel() - self.mean += x - - self.mean /= self.n_frames + if self.demean: + self._ref_atoms = self.reference.atoms + self.ref_cog = self.ref_atoms.center_of_geometry def _single_frame(self): - x = self._atoms.positions.ravel() - x -= self.mean + if self.demean: + # TODO make this more functional, initial proof of concept + mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, + self._ref_atoms.positons, + self._atoms, self._ref_cog, + self._atoms.center_of_geometry) + # now all structures are aligned to reference + x = mobile_atoms.positions.ravel() + else: + x = self._atoms.positions.ravel() self.cov += np.dot(x[:, np.newaxis], x[:, np.newaxis].T) def _conclude(self): @@ -194,11 +190,11 @@ def transform(self, atomgroup, n_components=None): ------- pca_space : array, shape (number of frames, number of components) """ - + # TODO has to accept universe or trajectory slice here if self._atoms != atomgroup: warnings.warn('This is a transform for different atom types.') - dot = np.zeros((u.trajectory.n_frames, ca.n_atoms*3)) + dot = np.zeros((self._u.trajectory.n_frames, ca.n_atoms*3)) for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): xyz = atomgroup.positions.ravel() From 4ed2070b9ca89ecf89c34827fd8fc325077d4a16 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 15:57:11 -0700 Subject: [PATCH 17/39] Iteration over crappy code to get things working. --- package/MDAnalysis/analysis/pca.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index ab75a1b19cf..17947bfb16a 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -115,7 +115,7 @@ class PCA(AnalysisBase): Take a pca_space and map it back onto the trajectory used to create it. """ - def __init__(self, atomgroup, demean=True, reference=None, + def __init__(self, atomgroup, select='All', demean_ref=None, n_components=None, **kwargs): """ Parameters @@ -138,31 +138,35 @@ def __init__(self, atomgroup, demean=True, reference=None, super(PCA, self).__init__(atomgroup.universe.trajectory, **kwargs) self._u = atomgroup.universe - self.reference = reference # for transform function self.start = kwargs.get('start') self.stop = kwargs.get('stop') self.step = kwargs.get('step') - self.demean = demean - self._atoms = atomgroup + if demean_ref: + self._demean = True + self._reference = demean_ref + else: + self._demean = False + + self._atoms = atomgroup.select_atoms(select) self.n_components = n_components self._n_atoms = self._atoms.n_atoms def _prepare(self): n_dim = self._n_atoms * 3 self.cov = np.zeros((n_dim, n_dim)) - if self.demean: - self._ref_atoms = self.reference.atoms - self.ref_cog = self.ref_atoms.center_of_geometry + if self._demean: + self._ref_atoms = self._reference + self._ref_cog = self._ref_atoms.center_of_geometry() def _single_frame(self): - if self.demean: + if self._demean: # TODO make this more functional, initial proof of concept mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, - self._ref_atoms.positons, + self._ref_atoms.positions, self._atoms, self._ref_cog, - self._atoms.center_of_geometry) + self._atoms.center_of_geometry()) # now all structures are aligned to reference x = mobile_atoms.positions.ravel() else: From dae191ae29ee9ece7f76294b3a0105e4e89f62bf Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 16:07:01 -0700 Subject: [PATCH 18/39] Get transform working, add calculated protected attribute. --- package/MDAnalysis/analysis/pca.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 17947bfb16a..b30490e09e5 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -139,10 +139,6 @@ def __init__(self, atomgroup, select='All', demean_ref=None, **kwargs) self._u = atomgroup.universe # for transform function - self.start = kwargs.get('start') - self.stop = kwargs.get('stop') - self.step = kwargs.get('step') - if demean_ref: self._demean = True self._reference = demean_ref @@ -152,6 +148,7 @@ def __init__(self, atomgroup, select='All', demean_ref=None, self._atoms = atomgroup.select_atoms(select) self.n_components = n_components self._n_atoms = self._atoms.n_atoms + self._calculated = False def _prepare(self): n_dim = self._n_atoms * 3 @@ -181,6 +178,7 @@ def _conclude(self): self.p_components = e_vects[sort_idx] self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance)) + self._calculated = True def transform(self, atomgroup, n_components=None): """Apply the dimensionality reduction on a trajectory @@ -194,11 +192,13 @@ def transform(self, atomgroup, n_components=None): ------- pca_space : array, shape (number of frames, number of components) """ + if not self._calculated: + self.run() # TODO has to accept universe or trajectory slice here if self._atoms != atomgroup: warnings.warn('This is a transform for different atom types.') - dot = np.zeros((self._u.trajectory.n_frames, ca.n_atoms*3)) + dot = np.zeros((self._u.trajectory.n_frames, atomgroup.n_atoms*3)) for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): xyz = atomgroup.positions.ravel() From c15d339c6ae0929ed077922d0597861a31133b51 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 28 Jul 2016 11:29:22 -0700 Subject: [PATCH 19/39] updated docs --- package/MDAnalysis/analysis/pca.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index b30490e09e5..af47cab9a55 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -102,7 +102,14 @@ class PCA(AnalysisBase): Percentage of variance explained by each of the selected components. If a subset of components is not chosen then all components are stored and the sum of explained variances is equal to 1.0. - + start: int + The index of the first frame to be used for the creation of the + covariance matrix. + stop: int + The index to stop before in the creation of the covariance matrix. + step: int + The amount of frames stepped between in the creation of the covariance + matrix. Methods ------- fit(traj=None, select='All', start=None, stop=None, step=None) From b22b2dfb886c376cee9de00124adcc1e5af8e733 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 28 Jul 2016 13:12:38 -0700 Subject: [PATCH 20/39] Reset previous commit, force pushing the appropriate changes --- package/MDAnalysis/analysis/pca.py | 19 +++++++++++++------ .../MDAnalysisTests/analysis/test_pca.py | 13 ++++++------- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index af47cab9a55..c9fbf3a820b 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -122,13 +122,19 @@ class PCA(AnalysisBase): Take a pca_space and map it back onto the trajectory used to create it. """ - def __init__(self, atomgroup, select='All', demean_ref=None, + def __init__(self, atomgroup, select='all', demean=True, n_components=None, **kwargs): """ Parameters ---------- atomgroup: MDAnalysis atomgroup AtomGroup to be used for PCA. + select: string, optional + A valid selection statement for choosing a subset of atoms from + the atomgroup. + demean: boolean, optional + If False, the trajectory will not be aligned to a reference + and there will be no demeaning of covariance data. n_components : int, optional The number of principal components to be saved, default saves all principal components, Default: -1 @@ -146,11 +152,12 @@ def __init__(self, atomgroup, select='All', demean_ref=None, **kwargs) self._u = atomgroup.universe # for transform function - if demean_ref: - self._demean = True - self._reference = demean_ref + if demean: + self._demean = demean + self._u.trajectory[0] + self._reference = self._u.atoms.select_atoms(select) else: - self._demean = False + self._demean = demean self._atoms = atomgroup.select_atoms(select) self.n_components = n_components @@ -205,7 +212,7 @@ def transform(self, atomgroup, n_components=None): if self._atoms != atomgroup: warnings.warn('This is a transform for different atom types.') - dot = np.zeros((self._u.trajectory.n_frames, atomgroup.n_atoms*3)) + dot = np.zeros((self._u.trajectory.n_frames, n_components)) for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): xyz = atomgroup.positions.ravel() diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 4c5ba06e387..98898f57da5 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -27,17 +27,17 @@ class TestPCA(object): def setUp(self): self.u = MDAnalysis.Universe(PDB, XTC) - self.pca = pca.PCA(self.u.atoms, select='backbone and name CA') + self.pca = pca.PCA(self.u.atoms, select='backbone and name CA', + demean=True) self.pca.run() - self.cumulated_variance, self.pcs = self.pca.fit() self.n_atoms = self.u.select_atoms('backbone and name CA').n_atoms def test_cum_var(self): - assert_almost_equal(self.cumulated_variance[-1], 1) + assert_almost_equal(self.pca.cumulated_variance[-1], 1) # makes no sense to test values here, no physical meaning def test_pcs(self): - assert_equal(self.pcs.shape, (self.n_atoms*3, self.n_atoms*3)) + assert_equal(self.pca.p_components.shape, (self.n_atoms*3, self.n_atoms*3)) def test_different_steps(self): self.pca.run() @@ -45,8 +45,7 @@ def test_different_steps(self): pcs = self.pca.p_components def test_transform(self): - self.n_components = 1 self.ag = self.u.select_atoms('backbone and name CA') - pca_space = self.pca.transform(self.ag, n_components=self.n_components) + pca_space = self.pca.transform(self.ag, n_components=1) assert_equal(pca_space.shape, - (self.u.trajectory.n_frames, self.n_components)) + (self.u.trajectory.n_frames, 1)) From 736d95a68eacd79cd8e9005d9a27ee63a925aae6 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 28 Jul 2016 14:15:31 -0700 Subject: [PATCH 21/39] Fixed transform function per Maxs comments --- package/MDAnalysis/analysis/pca.py | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index c9fbf3a820b..96db6ea98cf 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -194,14 +194,27 @@ def _conclude(self): np.sum(self.variance)) self._calculated = True - def transform(self, atomgroup, n_components=None): + def transform(self, atomgroup, n_components=None, start=None, stop=None, + step=None): """Apply the dimensionality reduction on a trajectory Parameters ---------- + atomgroup: MDAnalysis atomgroup + The atoms to be PCA transformed. n_components: int, optional - The number of components to be projected onto, default maps onto - all components. + The number of components to be projected onto, Default none: maps + onto all components. + start: int, optional + The frame to start on for the PCA transform. Default: None becomes + 0, the first frame index. + stop: int, optional + Frame index to stop PCA transform. Default: None becomes n_frames. + Iteration stops *before* this frame number, which means that the + trajectory would be read until the end. + step: int, optional + Number of frames to skip over for PCA transform. Default: None + becomes 1. Returns ------- pca_space : array, shape (number of frames, number of components) @@ -212,9 +225,11 @@ def transform(self, atomgroup, n_components=None): if self._atoms != atomgroup: warnings.warn('This is a transform for different atom types.') - dot = np.zeros((self._u.trajectory.n_frames, n_components)) + traj = atomgroup.universe.trajectory + start, stop, step = traj.check_slice_indices(start, stop, step) + dot = np.zeros((traj.n_frames, n_components)) - for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): + for i, ts in enumerate(traj[start:stop:step]): xyz = atomgroup.positions.ravel() dot[i] = np.dot(xyz, self.p_components[:n_components].T) From f318a3e370c34852bd02e72844723e7ecd41ac22 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 29 Jul 2016 13:14:38 -0700 Subject: [PATCH 22/39] Transform changes --- package/MDAnalysis/analysis/pca.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 96db6ea98cf..0718e1eeaa8 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -152,13 +152,11 @@ def __init__(self, atomgroup, select='all', demean=True, **kwargs) self._u = atomgroup.universe # for transform function - if demean: - self._demean = demean - self._u.trajectory[0] - self._reference = self._u.atoms.select_atoms(select) - else: - self._demean = demean - + self._demean = demean + # access 0th index + self._u.trajectory[0] + # reference will be 0th index + self._reference = self._u.atoms.select_atoms(select) self._atoms = atomgroup.select_atoms(select) self.n_components = n_components self._n_atoms = self._atoms.n_atoms @@ -235,7 +233,7 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, return dot - def inverse_tranform(self, pca_space): + def inverse_transform(self, pca_space): """ Transform PCA-transformed data back to original configuration space. Parameters @@ -249,5 +247,6 @@ def inverse_tranform(self, pca_space): ------- original_traj : array, shape (number of frames, number of atoms*3) """ - inv = np.linalg.inv(self.p_components) - return np.dot(pca_space, inv) + inv = np.linalg.pinv(self.p_components[:pca_space.shape[1]]) + original = np.dot(pca_space, inv.T) + return original From f0e6501c5cb9626c83775e9bcd0512f9549f9399 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 29 Jul 2016 14:29:55 -0700 Subject: [PATCH 23/39] Removed inverse_transform function, fixed rather large issue in transform function --- package/MDAnalysis/analysis/pca.py | 20 +------------------ .../MDAnalysisTests/analysis/test_pca.py | 7 +++++-- 2 files changed, 6 insertions(+), 21 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 0718e1eeaa8..e9b4a3d77b7 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -229,24 +229,6 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, for i, ts in enumerate(traj[start:stop:step]): xyz = atomgroup.positions.ravel() - dot[i] = np.dot(xyz, self.p_components[:n_components].T) + dot[i] = np.dot(xyz, self.p_components[:, :n_components]) return dot - - def inverse_transform(self, pca_space): - """ Transform PCA-transformed data back to original configuration space. - - Parameters - ---------- - pca_space : array - The space corresponding to the trajectory coordinates after the PCA - transformation. Assumes current principal components were the - components projected onto to create the pca_space. - - Returns - ------- - original_traj : array, shape (number of frames, number of atoms*3) - """ - inv = np.linalg.pinv(self.p_components[:pca_space.shape[1]]) - original = np.dot(pca_space, inv.T) - return original diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 98898f57da5..9b96da19346 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -46,6 +46,9 @@ def test_different_steps(self): def test_transform(self): self.ag = self.u.select_atoms('backbone and name CA') - pca_space = self.pca.transform(self.ag, n_components=1) - assert_equal(pca_space.shape, + self.pca_space = self.pca.transform(self.ag, n_components=1) + assert_equal(self.pca_space.shape, (self.u.trajectory.n_frames, 1)) + + # TODO test different inputs + # TODO test cosine content From e0ec9ff9dd16e39c02f43c07d93ad426e5ce3a30 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Mon, 1 Aug 2016 13:01:07 -0700 Subject: [PATCH 24/39] PCA covariance test --- .../MDAnalysisTests/analysis/test_pca.py | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 9b96da19346..3f584aad5c6 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -17,6 +17,8 @@ import numpy as np import MDAnalysis import MDAnalysis.analysis.pca as pca +from MDAnalysis.analysis.align import _fit_to + from numpy.testing import (assert_almost_equal, assert_equal, assert_array_almost_equal,raises) @@ -32,6 +34,24 @@ def setUp(self): self.pca.run() self.n_atoms = self.u.select_atoms('backbone and name CA').n_atoms + def test_cov(self): + atoms = self.u.select_atoms('backbone and name CA') + reference = atoms.positions + reference -= atoms.center_of_geometry() + ref_cog = atoms.center_of_geometry() + xyz = np.zeros((10, 642)) + for i, ts in enumerate(self.u.trajectory): + mobile_cog = atoms.center_of_geometry() + mobile_atoms, old_rmsd = _fit_to(atoms.positions, + reference, + atoms, + mobile_com=mobile_cog, + ref_com=ref_cog) + xyz[i] = mobile_atoms.positions.ravel() + + cov = np.cov(xyz, rowvar=0) + assert_array_almost_equal(self.pca.cov, cov, 5) + def test_cum_var(self): assert_almost_equal(self.pca.cumulated_variance[-1], 1) # makes no sense to test values here, no physical meaning From 6c11174405bd5d7b172391c04ade2753df36ca7c Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Mon, 1 Aug 2016 13:01:51 -0700 Subject: [PATCH 25/39] small changes --- package/MDAnalysis/analysis/pca.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index e9b4a3d77b7..ed563ce884b 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -168,14 +168,17 @@ def _prepare(self): if self._demean: self._ref_atoms = self._reference self._ref_cog = self._ref_atoms.center_of_geometry() + self._ref_atoms.positions -= self._ref_cog def _single_frame(self): if self._demean: # TODO make this more functional, initial proof of concept + mobile_cog = self._atoms.center_of_geometry() mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, self._ref_atoms.positions, - self._atoms, self._ref_cog, - self._atoms.center_of_geometry()) + self._atoms, + mobile_com=mobile_cog, + ref_com=self._ref_cog) # now all structures are aligned to reference x = mobile_atoms.positions.ravel() else: From 8e3f8abd96a39cbac574ada0f64c8022d0db69a5 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 3 Aug 2016 13:27:20 -0700 Subject: [PATCH 26/39] I came, I saw, I analyzed some principal components. --- package/MDAnalysis/analysis/pca.py | 57 +++++++++++++------ .../MDAnalysisTests/analysis/test_pca.py | 9 +-- 2 files changed, 45 insertions(+), 21 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index ed563ce884b..35f6e7391fe 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -122,7 +122,7 @@ class PCA(AnalysisBase): Take a pca_space and map it back onto the trajectory used to create it. """ - def __init__(self, atomgroup, select='all', demean=True, + def __init__(self, atomgroup, select='all', align=False, mean=None, n_components=None, **kwargs): """ Parameters @@ -132,9 +132,12 @@ def __init__(self, atomgroup, select='all', demean=True, select: string, optional A valid selection statement for choosing a subset of atoms from the atomgroup. - demean: boolean, optional - If False, the trajectory will not be aligned to a reference - and there will be no demeaning of covariance data. + align: boolean, optional + If True, the trajectory will be aligned to a reference + structure. + mean: MDAnalysis atomgroup, option + An optional reference structure to be used as the mean of the + covariance matrix. n_components : int, optional The number of principal components to be saved, default saves all principal components, Default: -1 @@ -152,30 +155,48 @@ def __init__(self, atomgroup, select='all', demean=True, **kwargs) self._u = atomgroup.universe # for transform function - self._demean = demean + self.align = align # access 0th index self._u.trajectory[0] # reference will be 0th index - self._reference = self._u.atoms.select_atoms(select) - self._atoms = atomgroup.select_atoms(select) + self._reference = self._u.select_atoms(select) + self._atoms = self._u.select_atoms(select) self.n_components = n_components self._n_atoms = self._atoms.n_atoms self._calculated = False + if mean is None: + self.mean = np.zeros(642) + else: + self.mean = mean.positions def _prepare(self): n_dim = self._n_atoms * 3 self.cov = np.zeros((n_dim, n_dim)) - if self._demean: - self._ref_atoms = self._reference - self._ref_cog = self._ref_atoms.center_of_geometry() - self._ref_atoms.positions -= self._ref_cog + self._ref_atom_positions = self._reference.positions + self._ref_cog = self._reference.center_of_geometry() + self._ref_atom_positions -= self._ref_cog + + for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): + if self.align: + mobile_cog = self._atoms.center_of_geometry() + mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, + self._ref_atom_positions, + self._atoms, + mobile_com=mobile_cog, + ref_com=self._ref_cog) + self.mean += mobile_atoms.positions.ravel() + else: + self.mean += self._atoms.positions.ravel() + + self.mean /= self.n_frames + # TODO if not quiet, inform user about double iteration + # TODO mean structure should be accessible as an atomgroup object def _single_frame(self): - if self._demean: - # TODO make this more functional, initial proof of concept + if self.align: mobile_cog = self._atoms.center_of_geometry() mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, - self._ref_atoms.positions, + self._ref_atom_positions, self._atoms, mobile_com=mobile_cog, ref_com=self._ref_cog) @@ -183,6 +204,8 @@ def _single_frame(self): x = mobile_atoms.positions.ravel() else: x = self._atoms.positions.ravel() + + x -= self.mean self.cov += np.dot(x[:, np.newaxis], x[:, np.newaxis].T) def _conclude(self): @@ -190,7 +213,7 @@ def _conclude(self): e_vals, e_vects = np.linalg.eig(self.cov) sort_idx = np.argsort(e_vals)[::-1] self.variance = e_vals[sort_idx] - self.p_components = e_vects[sort_idx] + self.p_components = e_vects[:, sort_idx] self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance)) self._calculated = True @@ -223,7 +246,7 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, if not self._calculated: self.run() # TODO has to accept universe or trajectory slice here - if self._atoms != atomgroup: + if self._atoms.atoms != atomgroup.atoms: warnings.warn('This is a transform for different atom types.') traj = atomgroup.universe.trajectory @@ -232,6 +255,6 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, for i, ts in enumerate(traj[start:stop:step]): xyz = atomgroup.positions.ravel() - dot[i] = np.dot(xyz, self.p_components[:, :n_components]) + dot[i] = np.dot(xyz, self.p_components[:,:n_components]) return dot diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 3f584aad5c6..d105db87eb0 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -30,27 +30,28 @@ class TestPCA(object): def setUp(self): self.u = MDAnalysis.Universe(PDB, XTC) self.pca = pca.PCA(self.u.atoms, select='backbone and name CA', - demean=True) + align=True) self.pca.run() self.n_atoms = self.u.select_atoms('backbone and name CA').n_atoms def test_cov(self): atoms = self.u.select_atoms('backbone and name CA') reference = atoms.positions - reference -= atoms.center_of_geometry() + ref_cog = atoms.center_of_geometry() + reference -= ref_cog xyz = np.zeros((10, 642)) for i, ts in enumerate(self.u.trajectory): mobile_cog = atoms.center_of_geometry() mobile_atoms, old_rmsd = _fit_to(atoms.positions, reference, - atoms, + mobile_atoms=atoms, mobile_com=mobile_cog, ref_com=ref_cog) xyz[i] = mobile_atoms.positions.ravel() cov = np.cov(xyz, rowvar=0) - assert_array_almost_equal(self.pca.cov, cov, 5) + assert_array_almost_equal(self.pca.cov, cov, 4) def test_cum_var(self): assert_almost_equal(self.pca.cumulated_variance[-1], 1) From 8938ff0a163b2006925d0e751aeaba133bc9d4bf Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 3 Aug 2016 14:15:03 -0700 Subject: [PATCH 27/39] Warnings added --- package/MDAnalysis/analysis/pca.py | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 35f6e7391fe..1eca19cf90a 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -140,7 +140,7 @@ def __init__(self, atomgroup, select='all', align=False, mean=None, covariance matrix. n_components : int, optional The number of principal components to be saved, default saves - all principal components, Default: -1 + all principal components, Default: None start : int, optional First frame of trajectory to use for generation of covariance matrix, Default: None @@ -154,6 +154,8 @@ def __init__(self, atomgroup, select='all', align=False, mean=None, super(PCA, self).__init__(atomgroup.universe.trajectory, **kwargs) self._u = atomgroup.universe + if self._quiet: + logging.disable(logging.WARN) # for transform function self.align = align # access 0th index @@ -165,9 +167,14 @@ def __init__(self, atomgroup, select='all', align=False, mean=None, self._n_atoms = self._atoms.n_atoms self._calculated = False if mean is None: + logger.warn('In order to demean to generate the covariance matrix\n' + + 'the frames have to be iterated over twice. To avoid\n' + 'this slowdown, provide an atomgroup for demeaning.') self.mean = np.zeros(642) + self._calc_mean = True else: self.mean = mean.positions + self._calc_mean = False def _prepare(self): n_dim = self._n_atoms * 3 @@ -176,20 +183,21 @@ def _prepare(self): self._ref_cog = self._reference.center_of_geometry() self._ref_atom_positions -= self._ref_cog - for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): - if self.align: - mobile_cog = self._atoms.center_of_geometry() - mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, - self._ref_atom_positions, - self._atoms, - mobile_com=mobile_cog, - ref_com=self._ref_cog) + if self._calc_mean: + for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): + if self.align: + mobile_cog = self._atoms.center_of_geometry() + mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, + self._ref_atom_positions, + self._atoms, + mobile_com=mobile_cog, + ref_com=self._ref_cog) self.mean += mobile_atoms.positions.ravel() else: self.mean += self._atoms.positions.ravel() - self.mean /= self.n_frames - # TODO if not quiet, inform user about double iteration + self.mean /= self.n_frames + # TODO mean structure should be accessible as an atomgroup object def _single_frame(self): From 32db2cc6d66faf9d37d7a1ad65c8b9d85697b1a2 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 4 Aug 2016 12:46:48 -0700 Subject: [PATCH 28/39] Fixed broken tests --- package/MDAnalysis/analysis/pca.py | 6 +++--- testsuite/MDAnalysisTests/analysis/test_pca.py | 5 +++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 1eca19cf90a..71c02b58487 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -192,9 +192,9 @@ def _prepare(self): self._atoms, mobile_com=mobile_cog, ref_com=self._ref_cog) - self.mean += mobile_atoms.positions.ravel() - else: - self.mean += self._atoms.positions.ravel() + self.mean += mobile_atoms.positions.ravel() + else: + self.mean += self._atoms.positions.ravel() self.mean /= self.n_frames diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index d105db87eb0..2b291364328 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -18,6 +18,7 @@ import MDAnalysis import MDAnalysis.analysis.pca as pca from MDAnalysis.analysis.align import _fit_to +from scipy.integrate import simps from numpy.testing import (assert_almost_equal, assert_equal, assert_array_almost_equal,raises) @@ -30,7 +31,7 @@ class TestPCA(object): def setUp(self): self.u = MDAnalysis.Universe(PDB, XTC) self.pca = pca.PCA(self.u.atoms, select='backbone and name CA', - align=True) + align=True) self.pca.run() self.n_atoms = self.u.select_atoms('backbone and name CA').n_atoms @@ -71,5 +72,5 @@ def test_transform(self): assert_equal(self.pca_space.shape, (self.u.trajectory.n_frames, 1)) + # TODO test different inputs - # TODO test cosine content From 9812f9e6d6383c39faa7cea89bb51a330f9b78bf Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 4 Aug 2016 13:17:07 -0700 Subject: [PATCH 29/39] Accepts a Universe as a transform argument, added tests that fails --- package/MDAnalysis/analysis/pca.py | 2 ++ testsuite/MDAnalysisTests/analysis/test_pca.py | 8 +++++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 71c02b58487..4f8be3cda61 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -256,6 +256,8 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, # TODO has to accept universe or trajectory slice here if self._atoms.atoms != atomgroup.atoms: warnings.warn('This is a transform for different atom types.') + if isinstance(atomgroup, Universe): + atomgroup = atomgroup.atoms traj = atomgroup.universe.trajectory start, stop, step = traj.check_slice_indices(start, stop, step) diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 2b291364328..4df8863e63b 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -73,4 +73,10 @@ def test_transform(self): (self.u.trajectory.n_frames, 1)) - # TODO test different inputs + # Accepts universe as input, but shapes are not aligned + @raises(ValueError) + def test_transform_inputs(self): + self.ag = self.u.select_atoms('backbone and name CA') + self.pca_space = self.pca.transform(self.u, n_components=1) + assert_equal(self.pca_space.shape, + (self.u.trajectory.n_frames, 1)) From 84fc3ce0f87bdec3118a778d706af54dc34dec02 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Mon, 8 Aug 2016 14:24:35 -0700 Subject: [PATCH 30/39] Made changes to pca to align with Scikit-learn, specifically demeaning in transform. Added random walk trajectory and tested for cosine content. Added some more functionality and testing Docs fixes and reverting some mistakes. Updated docs, added testing for universe as inputs, added mean positions as atom group. Some other style fixes. Updated CHANGELOG Progress meter added, switch to XTC --- package/CHANGELOG | 3 +- package/MDAnalysis/analysis/pca.py | 114 +++++++++++------- .../MDAnalysisTests/analysis/test_pca.py | 54 +++++---- .../MDAnalysisTests/data/RANDOM_WALK_TOPO.pdb | 105 ++++++++++++++++ .../MDAnalysisTests/data/xyz_random_walk.xtc | Bin 0 -> 57380 bytes testsuite/MDAnalysisTests/datafiles.py | 6 +- 6 files changed, 214 insertions(+), 68 deletions(-) create mode 100644 testsuite/MDAnalysisTests/data/RANDOM_WALK_TOPO.pdb create mode 100644 testsuite/MDAnalysisTests/data/xyz_random_walk.xtc diff --git a/package/CHANGELOG b/package/CHANGELOG index fa9d2be518d..7b215b398da 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -26,7 +26,8 @@ Enhancements class itself. * Added boolean flag at lib.distances.USED_OPENMP to reveal if OpenMP was used in compilation (Issue #883) - + * Added Principal Component Analysis module for linear dimensionality + reduction. Fixes * GROWriter resids now truncated properly (Issue #886) * reading/writing lambda value in trr files (Issue #859) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 4f8be3cda61..f4e8e9257f0 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -31,7 +31,7 @@ For each eigenvector, its eigenvalue reflects the variance that the eigenvector explains. This value is made into a ratio. Stored in -:attribute:`explained_variance`, this ratio divides the accumulated variance +:attribute:`cumulat_variance`, this ratio divides the accumulated variance of the nth eigenvector and the n-1 eigenvectors preceding the eigenvector by the total variance in the data. For most data, :attribute:`explained_variance` will be approximately equal to one for some n that is significantly smaller @@ -39,21 +39,17 @@ by Principal Component Analysis. From here, we can project a trajectory onto these principal components and -attempt to retrieve some structure from our high dimensional data. We have -provided a [notebook](# TODO edit max's notebook to use the new module) -containing a thorough demonstration of Principal Component Analysis. +attempt to retrieve some structure from our high dimensional data. For a basic introduction to the module, the :ref:`PCA-tutorial` shows how to perform Principal Component Analysis. - .. _PCA-tutorial: The example uses files provided as part of the MDAnalysis test suite (in the variables :data:`~MDAnalysis.tests.datafiles.PSF` and :data:`~MDAnalysis.tests.datafiles.DCD`). This tutorial shows how to use the PCA class. - First load all modules and test data :: >>> import MDAnalysis as mda >>> import MDAnalysis.analysis.pca as pca @@ -62,14 +58,20 @@ Given a universe containing trajectory data we can perform PCA using :class:`PCA`:: and retrieve the principal components. >>> u = mda.Universe(PSF,DCD) - >>> PSF_pca = pca.PCA(u) - >>> cumulated_variance, principal_components = PSF_pca.fit() + >>> PSF_pca = pca.PCA(u, select='backbone') + >>> PSF_pca.run() Inspect the components to determine the principal components you would like to retain. The choice is arbitrary, but I will stop when 95 percent of the -variance is explained by the components. - >>> n_pcs = next(x[0] for x in enumerate(cumulated_variance) if x[1] > 0.95) - >>> pca_space = PSF_pca.transform(n_components=n_pcs) +variance is explained by the components. This cumulated variance by the +components is conveniently stored in the one-dimensional array attribute +`cumulated_variance`. The value at the ith index of `cumulated_variance` is the +sum of the variances from 0 to i. + + >>> n_pcs = next(x[0] for x in enumerate(PSF_pca.cumulated_variance) + >>> if x[1] > 0.95) + >>> atomgroup = u.select_atoms('backbone') + >>> pca_space = PSF_pca.transform(atomgroup, n_components=n_pcs) From here, inspection of the pca_space and conclusions to be drawn from the data are left to the user. @@ -79,29 +81,37 @@ import warnings import numpy as np + +from MDAnalysis.core import AtomGroup from MDAnalysis import Universe from MDAnalysis.analysis.align import _fit_to +from MDAnalysis.lib.log import ProgressMeter from .base import AnalysisBase -logger = logging.getLogger(__name__) - class PCA(AnalysisBase): """Principal component analysis on an MD trajectory Attributes ---------- - p_components: array, (n_components, n_atoms) + p_components: array, (n_components, n_atoms * 3) The principal components of the feature space, representing the directions of maximum variance in the data. variance : array (n_components, ) The raw variance explained by each eigenvector of the covariance matrix. - explained_variance : array, (n_components, ) - Percentage of variance explained by each of the selected components. - If a subset of components is not chosen then all components are stored - and the sum of explained variances is equal to 1.0. + cumulated_variance : array, (n_components, ) + Percentage of variance explained by the selected components and the sum + of the components preceding it. If a subset of components is not chosen + then all components are stored and the cumulated variance will converge + to 1. + pca_space : array (n_frames, n_components) + After running :method:`pca.transform(atomgroup)` the projection of the + positions onto the principal components will exist here. + mean_atoms: MDAnalyis atomgroup + After running :method:`PCA.run()`, the mean position of all the atoms + used for the creation of the covariance matrix will exist here. start: int The index of the first frame to be used for the creation of the covariance matrix. @@ -112,14 +122,10 @@ class PCA(AnalysisBase): matrix. Methods ------- - fit(traj=None, select='All', start=None, stop=None, step=None) - Find the principal components of selected atoms in a subset of the - trajectory. - transform(traj=None, n_components=None) - Take the original trajectory and project it onto the principal - components. - inverse_tranform(pca_space) - Take a pca_space and map it back onto the trajectory used to create it. + transform(atomgroup, n_components=None) + Take an atomgroup or universe with the same number of atoms as was + used for the calculation in :method:`PCA.run()` and project it onto the + principal components. """ def __init__(self, atomgroup, select='all', align=False, mean=None, @@ -135,7 +141,7 @@ def __init__(self, atomgroup, select='all', align=False, mean=None, align: boolean, optional If True, the trajectory will be aligned to a reference structure. - mean: MDAnalysis atomgroup, option + mean: MDAnalysis atomgroup, optional An optional reference structure to be used as the mean of the covariance matrix. n_components : int, optional @@ -153,7 +159,12 @@ def __init__(self, atomgroup, select='all', align=False, mean=None, """ super(PCA, self).__init__(atomgroup.universe.trajectory, **kwargs) + if self.n_frames == 1: + raise ValueError('No covariance information can be gathered from a' + 'single trajectory frame.\n') + self._u = atomgroup.universe + if self._quiet: logging.disable(logging.WARN) # for transform function @@ -167,10 +178,10 @@ def __init__(self, atomgroup, select='all', align=False, mean=None, self._n_atoms = self._atoms.n_atoms self._calculated = False if mean is None: - logger.warn('In order to demean to generate the covariance matrix\n' - + 'the frames have to be iterated over twice. To avoid\n' - 'this slowdown, provide an atomgroup for demeaning.') - self.mean = np.zeros(642) + logging.warn('In order to demean to generate the covariance matrix\n' + 'the frames have to be iterated over twice. To avoid\n' + 'this slowdown, provide an atomgroup for demeaning.') + self.mean = np.zeros(self._n_atoms*3) self._calc_mean = True else: self.mean = mean.positions @@ -184,6 +195,10 @@ def _prepare(self): self._ref_atom_positions -= self._ref_cog if self._calc_mean: + interval = int(self.n_frames // 100) + interval = interval if interval > 0 else 1 + mean_pm = ProgressMeter(self.n_frames if self.n_frames else 1, + interval=interval, quiet=self._quiet) for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): if self.align: mobile_cog = self._atoms.center_of_geometry() @@ -195,10 +210,12 @@ def _prepare(self): self.mean += mobile_atoms.positions.ravel() else: self.mean += self._atoms.positions.ravel() - + mean_pm.echo(i) self.mean /= self.n_frames - # TODO mean structure should be accessible as an atomgroup object + atom_positions = self.mean.reshape(self._n_atoms, 3) + self.mean_atoms = AtomGroup.AtomGroup(self._atoms) + self.mean_atoms.positions = atom_positions def _single_frame(self): if self.align: @@ -212,7 +229,6 @@ def _single_frame(self): x = mobile_atoms.positions.ravel() else: x = self._atoms.positions.ravel() - x -= self.mean self.cov += np.dot(x[:, np.newaxis], x[:, np.newaxis].T) @@ -221,7 +237,8 @@ def _conclude(self): e_vals, e_vects = np.linalg.eig(self.cov) sort_idx = np.argsort(e_vals)[::-1] self.variance = e_vals[sort_idx] - self.p_components = e_vects[:, sort_idx] + self.variance = self.variance[:self.n_components] + self.p_components = e_vects[:self.n_components, sort_idx] self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance)) self._calculated = True @@ -232,8 +249,8 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, Parameters ---------- - atomgroup: MDAnalysis atomgroup - The atoms to be PCA transformed. + atomgroup: MDAnalysis atomgroup/ Universe + The atomgroup or universe containing atoms to be PCA transformed. n_components: int, optional The number of components to be projected onto, Default none: maps onto all components. @@ -253,18 +270,29 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, """ if not self._calculated: self.run() - # TODO has to accept universe or trajectory slice here - if self._atoms.atoms != atomgroup.atoms: - warnings.warn('This is a transform for different atom types.') + if isinstance(atomgroup, Universe): atomgroup = atomgroup.atoms + if(self._n_atoms != atomgroup.n_atoms): + raise ValueError('PCA has been fit for' + '{} atoms. Your atomgroup' + 'has {} atoms'.format(self._n_atoms, + atomgroup.n_atoms)) + if not (self._atoms.types == atomgroup.types).all(): + warnings.warn('Atom types do not match with types used to fit PCA') + traj = atomgroup.universe.trajectory start, stop, step = traj.check_slice_indices(start, stop, step) - dot = np.zeros((traj.n_frames, n_components)) + n_frames = len(range(start, stop, step)) + + dim = (n_components if n_components is not None else + self.p_components.shape[1]) + + dot = np.zeros((n_frames, dim)) for i, ts in enumerate(traj[start:stop:step]): - xyz = atomgroup.positions.ravel() - dot[i] = np.dot(xyz, self.p_components[:,:n_components]) + xyz = atomgroup.positions.ravel() - self.mean + dot[i] = np.dot(xyz, self.p_components[:, :n_components]) return dot diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 4df8863e63b..3da1bc8bfd7 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -21,50 +21,39 @@ from scipy.integrate import simps from numpy.testing import (assert_almost_equal, assert_equal, - assert_array_almost_equal,raises) + assert_array_almost_equal, raises) -from MDAnalysisTests.datafiles import PDB, XTC - +from MDAnalysisTests.datafiles import PDB, XTC, RANDOM_WALK, waterPSF, waterDCD class TestPCA(object): def setUp(self): self.u = MDAnalysis.Universe(PDB, XTC) self.pca = pca.PCA(self.u.atoms, select='backbone and name CA', - align=True) + align=False) self.pca.run() self.n_atoms = self.u.select_atoms('backbone and name CA').n_atoms def test_cov(self): atoms = self.u.select_atoms('backbone and name CA') - reference = atoms.positions - - ref_cog = atoms.center_of_geometry() - reference -= ref_cog xyz = np.zeros((10, 642)) for i, ts in enumerate(self.u.trajectory): - mobile_cog = atoms.center_of_geometry() - mobile_atoms, old_rmsd = _fit_to(atoms.positions, - reference, - mobile_atoms=atoms, - mobile_com=mobile_cog, - ref_com=ref_cog) - xyz[i] = mobile_atoms.positions.ravel() + xyz[i] = atoms.positions.ravel() cov = np.cov(xyz, rowvar=0) assert_array_almost_equal(self.pca.cov, cov, 4) def test_cum_var(self): assert_almost_equal(self.pca.cumulated_variance[-1], 1) - # makes no sense to test values here, no physical meaning def test_pcs(self): - assert_equal(self.pca.p_components.shape, (self.n_atoms*3, self.n_atoms*3)) + assert_equal(self.pca.p_components.shape, + (self.n_atoms*3, self.n_atoms*3)) def test_different_steps(self): - self.pca.run() - cum_var = self.pca.cumulated_variance - pcs = self.pca.p_components + dot = self.pca.transform(self.u.select_atoms('backbone and name CA'), + start=5, stop=7, step=1) + assert_equal(dot.shape, (2, self.n_atoms*3)) def test_transform(self): self.ag = self.u.select_atoms('backbone and name CA') @@ -72,11 +61,30 @@ def test_transform(self): assert_equal(self.pca_space.shape, (self.u.trajectory.n_frames, 1)) - - # Accepts universe as input, but shapes are not aligned + # Accepts universe as input, but shapes are not aligned due to n_atoms @raises(ValueError) - def test_transform_inputs(self): + def test_transform_mismatch(self): self.ag = self.u.select_atoms('backbone and name CA') self.pca_space = self.pca.transform(self.u, n_components=1) assert_equal(self.pca_space.shape, (self.u.trajectory.n_frames, 1)) + + def test_transform_universe(self): + u1 = MDAnalysis.Universe(waterPSF, waterDCD) + u2 = MDAnalysis.Universe(waterPSF, waterDCD) + pca_test = pca.PCA(u1).run() + pca_test.transform(u2) + + def test_cosine_content(self): + rand = MDAnalysis.Universe(RANDOM_WALK) + pca_random = pca.PCA(rand.atoms).run() + dot = pca_random.transform(rand.atoms) + content = cosine_content(dot, 0) + assert_almost_equal(content, .99, 1) + + +def cosine_content(pc, i): + t = np.arange(len(pc)) + T = len(pc) + cos = np.cos(np.pi * t * (i + 1) / T) + return (2.0 / T) * (simps(cos*pc[:, i])) ** 2 / simps(pc[:, i] ** 2) diff --git a/testsuite/MDAnalysisTests/data/RANDOM_WALK_TOPO.pdb b/testsuite/MDAnalysisTests/data/RANDOM_WALK_TOPO.pdb new file mode 100644 index 00000000000..ace4d7dd241 --- /dev/null +++ b/testsuite/MDAnalysisTests/data/RANDOM_WALK_TOPO.pdb @@ -0,0 +1,105 @@ +TITLE MDANALYSIS FRAMES FROM 0, SKIP 1: Created by PDBWriter +CRYST1 0.000 0.000 0.000 0.00 0.00 0.00 P 1 1 +MODEL 1 +ATOM 1 X SYST 1 0.192 -0.797 -1.042 1.00 0.00 X +ATOM 2 X SYST 1 -0.150 0.308 0.957 1.00 0.00 X +ATOM 3 X SYST 1 -0.752 -0.125 -0.563 1.00 0.00 X +ATOM 4 X SYST 1 0.483 -1.342 -1.400 1.00 0.00 X +ATOM 5 X SYST 1 0.550 -0.672 -1.462 1.00 0.00 X +ATOM 6 X SYST 1 -0.644 -0.858 0.751 1.00 0.00 X +ATOM 7 X SYST 1 1.092 1.182 -1.172 1.00 0.00 X +ATOM 8 X SYST 1 -1.028 1.086 0.098 1.00 0.00 X +ATOM 9 X SYST 1 -2.109 0.854 -0.326 1.00 0.00 X +ATOM 10 X SYST 1 -0.430 3.140 0.412 1.00 0.00 X +ATOM 11 X SYST 1 -1.051 0.127 -0.654 1.00 0.00 X +ATOM 12 X SYST 1 1.934 1.939 -0.671 1.00 0.00 X +ATOM 13 X SYST 1 -1.622 2.113 -1.297 1.00 0.00 X +ATOM 14 X SYST 1 -0.964 -1.267 -1.626 1.00 0.00 X +ATOM 15 X SYST 1 0.079 -1.222 -1.151 1.00 0.00 X +ATOM 16 X SYST 1 1.328 1.413 -1.109 1.00 0.00 X +ATOM 17 X SYST 1 2.707 -3.728 2.036 1.00 0.00 X +ATOM 18 X SYST 1 -1.879 -1.751 -0.008 1.00 0.00 X +ATOM 19 X SYST 1 0.509 0.497 0.657 1.00 0.00 X +ATOM 20 X SYST 1 -1.859 -1.078 1.447 1.00 0.00 X +ATOM 21 X SYST 1 -1.484 0.366 1.113 1.00 0.00 X +ATOM 22 X SYST 1 1.004 -0.138 -1.316 1.00 0.00 X +ATOM 23 X SYST 1 -1.070 -1.564 -0.288 1.00 0.00 X +ATOM 24 X SYST 1 -1.634 -1.837 1.273 1.00 0.00 X +ATOM 25 X SYST 1 -0.578 -1.060 -1.521 1.00 0.00 X +ATOM 26 X SYST 1 0.855 1.954 -0.705 1.00 0.00 X +ATOM 27 X SYST 1 1.089 -2.214 0.386 1.00 0.00 X +ATOM 28 X SYST 1 0.654 -0.547 0.125 1.00 0.00 X +ATOM 29 X SYST 1 1.196 0.560 -1.884 1.00 0.00 X +ATOM 30 X SYST 1 1.537 -1.519 1.516 1.00 0.00 X +ATOM 31 X SYST 1 -0.333 0.382 -1.087 1.00 0.00 X +ATOM 32 X SYST 1 0.758 -0.508 0.637 1.00 0.00 X +ATOM 33 X SYST 1 0.912 -0.132 -0.114 1.00 0.00 X +ATOM 34 X SYST 1 0.514 -2.753 -1.568 1.00 0.00 X +ATOM 35 X SYST 1 -2.796 0.934 1.568 1.00 0.00 X +ATOM 36 X SYST 1 -0.316 -0.149 -1.460 1.00 0.00 X +ATOM 37 X SYST 1 1.589 -1.918 1.128 1.00 0.00 X +ATOM 38 X SYST 1 0.893 0.243 1.651 1.00 0.00 X +ATOM 39 X SYST 1 -1.718 -1.066 1.226 1.00 0.00 X +ATOM 40 X SYST 1 0.329 0.702 0.037 1.00 0.00 X +ATOM 41 X SYST 1 -0.093 -1.030 -0.573 1.00 0.00 X +ATOM 42 X SYST 1 -0.245 0.664 1.448 1.00 0.00 X +ATOM 43 X SYST 1 1.946 -2.568 0.803 1.00 0.00 X +ATOM 44 X SYST 1 -2.092 -1.243 -0.670 1.00 0.00 X +ATOM 45 X SYST 1 -0.315 0.386 -0.957 1.00 0.00 X +ATOM 46 X SYST 1 -2.788 -1.573 -0.351 1.00 0.00 X +ATOM 47 X SYST 1 -0.680 -1.665 0.872 1.00 0.00 X +ATOM 48 X SYST 1 1.844 1.731 0.719 1.00 0.00 X +ATOM 49 X SYST 1 -0.008 2.494 -1.046 1.00 0.00 X +ATOM 50 X SYST 1 0.594 0.757 -0.859 1.00 0.00 X +ATOM 51 X SYST 1 -1.617 -0.093 2.723 1.00 0.00 X +ATOM 52 X SYST 1 -0.160 1.120 2.035 1.00 0.00 X +ATOM 53 X SYST 1 -0.241 -1.276 1.334 1.00 0.00 X +ATOM 54 X SYST 1 -0.650 -0.517 0.589 1.00 0.00 X +ATOM 55 X SYST 1 0.140 -1.181 -2.026 1.00 0.00 X +ATOM 56 X SYST 1 0.263 0.025 -0.591 1.00 0.00 X +ATOM 57 X SYST 1 -2.602 0.590 -0.103 1.00 0.00 X +ATOM 58 X SYST 1 1.164 -1.947 -0.433 1.00 0.00 X +ATOM 59 X SYST 1 0.335 -2.007 -1.063 1.00 0.00 X +ATOM 60 X SYST 1 3.312 -0.641 0.481 1.00 0.00 X +ATOM 61 X SYST 1 0.318 -0.649 -1.398 1.00 0.00 X +ATOM 62 X SYST 1 0.846 1.393 -1.362 1.00 0.00 X +ATOM 63 X SYST 1 2.293 1.473 0.462 1.00 0.00 X +ATOM 64 X SYST 1 -0.713 -1.557 1.726 1.00 0.00 X +ATOM 65 X SYST 1 2.920 -1.216 -0.189 1.00 0.00 X +ATOM 66 X SYST 1 -1.824 1.101 0.334 1.00 0.00 X +ATOM 67 X SYST 1 0.151 -0.205 -0.040 1.00 0.00 X +ATOM 68 X SYST 1 -1.562 0.608 -0.604 1.00 0.00 X +ATOM 69 X SYST 1 2.523 1.642 0.167 1.00 0.00 X +ATOM 70 X SYST 1 -0.219 0.873 -0.553 1.00 0.00 X +ATOM 71 X SYST 1 1.307 -0.766 0.101 1.00 0.00 X +ATOM 72 X SYST 1 -0.034 -2.097 0.004 1.00 0.00 X +ATOM 73 X SYST 1 0.315 0.406 -0.990 1.00 0.00 X +ATOM 74 X SYST 1 -1.116 -0.050 -0.500 1.00 0.00 X +ATOM 75 X SYST 1 1.444 0.385 -0.008 1.00 0.00 X +ATOM 76 X SYST 1 1.172 1.610 -0.969 1.00 0.00 X +ATOM 77 X SYST 1 0.605 -0.021 -1.358 1.00 0.00 X +ATOM 78 X SYST 1 -2.316 0.837 0.168 1.00 0.00 X +ATOM 79 X SYST 1 0.430 -1.242 0.087 1.00 0.00 X +ATOM 80 X SYST 1 2.290 -0.587 -1.829 1.00 0.00 X +ATOM 81 X SYST 1 0.144 -0.724 -0.735 1.00 0.00 X +ATOM 82 X SYST 1 0.756 -0.873 -0.661 1.00 0.00 X +ATOM 83 X SYST 1 1.778 1.564 1.385 1.00 0.00 X +ATOM 84 X SYST 1 -1.414 1.198 0.396 1.00 0.00 X +ATOM 85 X SYST 1 -0.075 1.270 1.709 1.00 0.00 X +ATOM 86 X SYST 1 -3.111 -1.338 -1.163 1.00 0.00 X +ATOM 87 X SYST 1 0.021 3.037 -0.199 1.00 0.00 X +ATOM 88 X SYST 1 -1.968 -2.604 0.683 1.00 0.00 X +ATOM 89 X SYST 1 -0.082 1.604 -1.400 1.00 0.00 X +ATOM 90 X SYST 1 -0.772 1.781 -1.497 1.00 0.00 X +ATOM 91 X SYST 1 0.159 1.882 0.874 1.00 0.00 X +ATOM 92 X SYST 1 -2.455 -0.851 0.206 1.00 0.00 X +ATOM 93 X SYST 1 -0.461 1.583 3.403 1.00 0.00 X +ATOM 94 X SYST 1 1.318 -1.974 -0.323 1.00 0.00 X +ATOM 95 X SYST 1 0.742 -1.308 -1.363 1.00 0.00 X +ATOM 96 X SYST 1 1.760 0.902 1.184 1.00 0.00 X +ATOM 97 X SYST 1 1.740 -0.107 -1.524 1.00 0.00 X +ATOM 98 X SYST 1 0.084 -1.835 -0.533 1.00 0.00 X +ATOM 99 X SYST 1 2.022 2.084 0.112 1.00 0.00 X +ATOM 100 X SYST 1 -0.454 0.563 2.780 1.00 0.00 X +ENDMDL +END diff --git a/testsuite/MDAnalysisTests/data/xyz_random_walk.xtc b/testsuite/MDAnalysisTests/data/xyz_random_walk.xtc new file mode 100644 index 0000000000000000000000000000000000000000..39a68df7b7d0fc88e3627c2790baa82ecce5fa57 GIT binary patch literal 57380 zcmb5WWl$YKx3-Nt1b26Lf)m``o#0Mzx8Uv&Jh;1iaEIU;+}#P5@J@1`teV$!iW zeH^btZ}4oJx9pz%}R$4d{;{?sTaTbPw|e0X-dAQ`RI>pF(p>BMa;q9kIh!}^JvT`*Dd7!wp5%b3S3##N2!kyhc08mA{#A7wo7$Vx7&Ng%kY9F&Qg^nEhxj&+uRneOx$Fi30Yd`% zULXK6(1+iQz}R8##o1%--qIy<^{5I?2X~xvO7pQdPGoGjYum3dN%uyOnN2EgD}2^J zDt6seA+&P(L~rFX+Lzftdyh%OVb<8*H=x(zoU+`A$EOFGv;F2)Yew`EIl^aaaPaCT zp7mp}U|da;apu{T zC9+OxBU)C?1qWo#?ftbff1XT9yb0;c^q9X^%!FqA`^G^+u$3#Hql0g_7 zyF;jzAvv6Yv^2|b#PQU)Chif4&93OZQwyc@QW|CH@x)eqGzo~EpKo1DRDob2i;uF> z(#Rde;!uiK@E(cX?hm2NRAW`BJ*zF23YqYYmX7l+r$>cc`|*s(w>^3RY&$;OLHZfo1iI%&(5Z?Muk%`_3Ol5oc(wSO!Vww?7P*3-rYPp%Ma6!E(u8>~6!r}s;y;5I zU-dHXHS9~UU_Q^DuJ@V9Db@87TdK!|mz6CXBJE785-o%4sx;(lruD0AIrVQZ3=m%m*wn&XI%I8+n@Zq1G)bR2=I1r4?y97a)1~(K>jNNJ?_#(-Up$BM7VHw)gpf@ zf^ejso{-EZ@MZNuYbidPp_NJUD1?P}u-&KiCoPm2eu)UwWK}buVn_T)%j57U+l1l8 zAQqf3gD#PgQKZkyhNga~UN{i@vyj9sMVJP1kao9e98RX)y}yq*Pw_x7dMQ-dVEMcJ zg-r0f3GU7~i?8@M+aXm+p^%FS=C+7T2QMq)TTuy4)MM}HgHzHs2jegK!D{%n0s^;p zf}|Lx$vT?cu?obSmsQwamKdwVN1}LCm+s$~rN|SkbD%LaJWHK*u}el0olT)qo4Wt( z-4T|}GC#xB>X1-UGR1n=65QhknLGG0MehxT0Nvlv$=^J8e+|>?rJU6|B7Dn&WD$uB z#e4kWh_NQa*#V-OuG?5HmO(4yoE4J8CC=NHyCqI-4K>6 zAm(ezv~c2RTeu>YdQ#LTMQ%>1h4x*xycrA$)NB`08!>88UOw;vHCeUizJMbN)&JrE zQt(#>|E>Mu01aRu8t{I45D*9hfPs&If`Aa100;HB?6OJ1A_do^RvXH4pWR^Wi?7Z3 z1fA_RvLm0Wn))k+UXwVT9_moourcOItsk0F2@`2gl<`6AX4Y7lhpdBdUkGhLqk7PdMwF#<*;H2Xr&_$s_l0k(XWVR zaDF-In!61aMx+^6M;Gl-Ga^9Gz@2-s#EVJ-wrO$oJzYlLM;0KY7c;uNrvYs`#iGtR`e(Ig{NH4n^ z8#NUf3!zf2ZavvSLT#wbTgqdML3O4AE%zt3{`GUOB`X^{LZy~

Wx@H<0y3V)2e6 z_Ls2HErO(0e6uKF${@6P z6xonZK5~5}9KEn5Q_{rlpD#c}MOOal3|_7}2U(|bcTE4eRX0g66ZUfC-O!B99mh`{ zkrKS0KRNh8N>dOc-W%OXer-G$Q4u6wtR(1I< z45aSwS^!G0``=&yX8nhUXCZ(GOh9%ZAdnCM3&j8zaKPuk)&i}PeI(F^s0oHVm~F@Y zGNHjUQ{i>FoR`oo+!3_|wcZFUh8H!_)n8r}l*xrElFaC6)jy-~iug+B%b5#e2d$Aq zEtBe4NXU|pvF1K_ov0e?7tv@?L2ZHNXt%2W3|j6fA-cUJdCKMn0V(LQJhj$7W42o! zqOK)}V~b$u;72yovPe8It$U^&77iU#z%x1cB7DQDa2YY4se?E$^WK-X&LlRlgXfhB@vwPF%M}7Qwfd(`rodbse-L_VId#`Zw|0Y;MdHe^nYX(pgOx1U zl4i4P4-w$gbW7B8r^=-nfurOv*sEC!TFyT~o1QhVakozUtzUgnKAE~bRPYZFUF{x~ ziF*g@n8IRe(iYXl#fu6Rgs44IVen}5jWEcOhL^F!c_3$R9e)wodS~Q-&r{2nSe!wc zR{H5{q|OZFj_j+AiM?cJ&K}9$sWOd;^xPMMYGt-^G%F4^Lb-Bd#t&dmVj^vs?AI69 z=*c-qmmD4{NZ6uKoW#E|_{#;f;O6fe|8n`m#Ycb%6hP1b52XMX@PJ-`0Jy(=f|6?a zyGZ1KoqD+7{C7lT{y2y>22>yY>lN-)B9w-Cw>A=;9fX0-?}ADan;&k&8%#@m%NUcE zdrUjR^uD+^Grm+w?V3=1-ACC$AhTWUPs9^U?}NrM!#+j%S=|l-v0wqRTeY7(_Hk?@ zB2Dp~;9h1A{_c$USFwY4PRrV*1wn-w29%*E4)52YT;HzI;#o-EJkyb>?hm5kyggnm z5*>C^tC9pat~PAjG3^rm3 zuJwykm#U?6)lBgW;p~XVdcbhy_w?rIL!IC3^yzj=X@X?T;m1AmohT->&Vb^Ba?-!(wu5GJc-$qIF1T%P7GK6C?{zD=kTwVW&I)zZyZIpHZFb$@w)5&9Pr{%6jAcz9X_ScrOkeI5XO!Vq9#65yf+ znDT4=(BU&%tY`RvTbmx`BW_=Xg8ejhfjV9`R4G`w7)ivO*!Xcu&H9Ss1$O{YU(v`nKaRnf=g)@34$jTC)}?ebG&Tq%l7p^n}@qDj+r z^wQ1{`nZJr-KMuJ?Hw;ywHZwZsM_{6SXz@gCLaB0qGKX>M0g?(umEX-#2gUqSxl|F zf!`$~3hkq*#DTNzDmtkU03V=JSuEl8(zfxzAi1B8(Z+dxPTQx{%0`$_yn-7dH}2e`NgxbOj@-~n8;bK0%N7}}bhW0N|eVdEN2O=UNdBKp`Nm?j|-#xB?0vqVM) zWU%MuT;S4W=npwbuPLh~<4`b*oHw>3O@G`PU%G@Kj-N24JD!B(rlbUpTyeyUaszh_0|uJ?#yC?9Q-Yr=2N@Z=;qZ7 z6FNbagokbUN>KnAXB3jNF0}v#+L#+(5drZh3HpPp*o_6_n7-+}Jhd={#ScoqAKh z{m1}SQL$!6y@NCCyO4Ka9Lu~PiKQCnHjfyeiE{F!Mw5$Y5V#o5n~bqD{>|28?n8o< zyF1!^NFwA_AZ;e)-jKbNuW?ZF4QLFhrS$O(6_?7FiwTmoTj$Qn0*k1Sx4Pb^yss*) zEZ77fJ*S{w1JOo#Al3Y{MkqD&zT{P#IAf(a?x)qOEldtzySOd4wPUA13|u$jcOCnL z?s53sb0B>>Xd^!f?89(!n`nHDOy`tHMhwCFFwQopdRWjY1(M6If$LSmCJ&_WuV26k z5d*`2Yk#;v1(?tSWCU@@LybcHzAqBir zc8H*+h&umMnu02A14iw*>WGnr<V<$~XM zxyrsrf>K`6KQ5W#)~{?qOw%?g^oBF`*1)>thEaaU;#j*#MmdBg-O>ruJAU~v5^ zO9Ug4x9ro)8mME_VQSVfKUTjs30IPXwJNT71DT@FTM9`IvUg&b6F*L1Q#d!Z?j14V zX87KUe{ZKIYc2^Y3s2RXdw8q|gDqq9cmsDP)Bq+?jPWU{8%C26$@Y0AB#Pp5-3pan z-Rl97YTkl_0qDdIwqZ+2;yQPP_~(bmQrIli_B7kDTXN}{wejcgw|K#ggkn?b5s zwomP9z&6DHase;I{BKYIv;NePW`GGoK+ymXw}4l$fq*~(&osZe=!=I_&}Xi9R%L=| zX3rX!@~EtuVveaDQ~dUmbpj$t2FncApf zyPLXD5_;Rcq^j0!3|l+<`hZ+ONf~`3dsPKddc3-p!<+#be(Vg|zXRVWu@fdyD!@hj zu!r@e2p&{L@+odpIX;Jbc|!p|EFyR4yt9MV%AbFwLoeWa49a8XiL+4_XH&&2s~W_0 zZ$Y;6Bxr#RPrKsJaxEG3&m5zMdIl-#al}ufWSW>7Ohp7reW+h3OYOEYnh1|z&5G>l z^Q|CN*Jm@lb_%CL2=l2k=F#S{Yd}!j%MsHh&Yyuz4Vl+7`Hp@V8TUP@(9Q6)cmfo= zx+vk)&i#W_q(xbhGNPa%0a}PQzfODjGv#8!j*nFm{3gXZw-?igP!d*kFwpdv++1C;Y*39Wxk!>Du*Kcc%R%i-qL@tCVOoUFFpc} zQI_{@UJj#z^w=g8{C-df;v);#C~|Ca`maaHzrKMWBnS-ut^MI64PZkQP!hm}B*2F< z(2oKkem##2(jHnDFrS6Qfk1z6-%=*LsaS?)v|O{&T#?{QEzln_n-lJR4i}U0qO~^? zZg~|`qV$xwNm|fL>J=j?mclM`)n%GaleoKs)P*_RUifzNlF0jQ3SBwbfNB^-1sD3Me5vS^>J2k(fnFnp;vHg@TJ*>2<&WRs>}P90PAH73O#7{Pb~x z4|A>M%3w&IBEmTDjL*Z9Ysa6>r?=V5>AtE7FOu|u2#bk2i80Bxo9MVUT#Yl=TXKbP z8P!=MeC}nZ?6>l2J#ECCmk5Q3++^PnHo1^+g5AAM@npC@g@6cxplacd3`|Yy*;_Bz z|E{`^YfR)fPHPo93S$agI zq}wU$_7(EVVcfLkVOm>>)QfNn`>oA|%zo45c8w4DJ{5y^9|mEMflVCIh(LPKXs?SF zlfWW0Y>*Vb)`&tot8vi~VyjFuS5dVVgy@m*t$qqu!LvSE&&<(BKenf*q>jDatTSsZ zqycON0kys9AefikEgDov7(-k^f9i)|)dAaLT;$ch2P~W^O5@`vMwhE0-#A78UqJk2 z1W`!w*Z7~>e;9eZ1o)5!qz-t<6(HctkoCaj7b64uv7>Fw8RHW0pt}#?v_4!C4ndv{ zFd0-7Gi9H0Y;wjxJp4IQ&^ypRhR<<}M~g;5VI+uL4yKhRGaJ?!t=HdADY=fgc;{3L z5sV{DK+}`@XXwmO*+PqO)ekX`D$?Ct;fNJ}gETRi>;*LmnrwngnzWyjWs|%~NlYIT zmc95c=$vprpgercl+vOS?)bf|%v0ZCLefVjuV^f0WBEu+v@0@s&=me*lPjcJ7=Md) zKfjO~;rf6zKqYCDQD*ppQeS&>Q{kdIN5U!!SBW5x8FEDu`6T88-(e*^bsJg;TgLAE z@bTGuoRUjM5Rz6hCrc87tuPTq)LUIjuGfz5s};~TI2fXJA>CZxRckBNvU|RTa4*Fo z-TICVE5UQyve_(It_@i>v6Mm~3d9u0h_@GDb!c9;3YbFLl(0?DZQXo4!ZwOG4@WIG zJ-KAWvl`jM4gNHCZY(E}7hD?Dv&T*BJS(d!HV%uq=$xSg63}d3UVPP&Rl|cDV1J>n=j{FoNvzu zXX!P(6H*c7aNcF*LgVTy{Xfglzl+fDNEZL68IeEAaWR zb?G3E1x*_BHX|J}Oak6RFCP0v!{GqiZk+pej@N=2ix`cLcuSNtHDA=L!6nia+$gs)lH(A{eqK_{*rZM&~h?+R}&HZZ~)TngQq zL+5jxH_Dhf((K|Wj5Y|Ccyq!(4J_p{Uw)^gXVj%rZy?0Ign;}QwCWyHkm*g{@R3}* zq*I&g6An}WHpI1YJn`)MI)QXUXz839IW;?o&yMN10lgo`u)M;Onz#1(nEsr}r5uBw zE6In(MfwXN`8zlX*u@m#w5^dX=+ZVe1FNmcE_Tsqs1Hnu*=aA$5*y~QKOhkbec?@O zVKzSbm|y9Rb&zmm37Nc7Wyhfx?=oPey|=t0bKvPR17G4a9XwfADjbOTmF~Zcyb-eh zHU6jeA4VP?0X}8{0qYNlO<;|B4b&P|Ai}RbjtT!HMk2+V)`cjT2pe`h++rD&uu_U< zW?aRBa_@mik!kO3aA64@;m^fv$(W&T!4VNss}GNG2A{5yD*9jyG>K|lzExmX9xW*Q zzMCOm#>D3!?&0!Y^}F;zD%JyAJ3NGYUJzHDHEew2_JN=We%LKg3dfKvH^e-P=ldIY z61SCvZDc&~7`E?+*7WU1v}C@MvMC~A;JNVk{y4BPtXS>IuMcI{N!$*IWeMqFtIE>cXQMD_ z+u<^T+RR z;Tepc3A7gx8YolkAw$3Lmd6$~+7W`R;_=vLh0+<^0~!;u|bU(%^vH2|4#=#@duW!i{0v3QV!o~2gJ(v&SKjQw^H zes@a4)jlg-gwr**8ng*AMJ?Ag_F1=jsrxq;e_25m^8bC~UoLU0qAzAYpcd%M3S1|I?naRU(6HBDT6tf?lw>3z2_C)NcAMVhJ zK?@@yRPd)CMaiG{bju@Lxxb^j*C^l}EpM80LuSsEc?Z>g0()pRYy zWU@cnTgle_Km)dQnd3^jG`r8$KgW#k^nvo9i%{iWJ5XTnY_Jj~CWFAJ)zLX&dUJVJ z4;x|L>NbWQUl`*mD}F;nfJ993j)PRr86|GA0SiiIM0jv^&|#5e+v{>O0Ycr~BZ@&6%wB649gk7{o%p?jdH-M+iyQ!%C<6xnGvG z#hbAyKn?6Nn7xVN_5SC6SiU;#F7I^nB>RYfFyxtdyVgKBy$+pnMT{;SH?Lah!u>JF z$8rZ#X(?K0bz9yb-dOXJRgL&+t2S9%o^BZ5)Df;GlSBh5&wVRBD6`jv;@;31*&Rk| zfQi4Hpa_Ni3nIVg{P7eefD>9k=>Jo1Ab?Xdkig!^@3qNTIBgN8;?ryVm!LZ}H!?%H z20Y3VO8X@5S0@RQH|7I$?>?qr)o=zCp+eSDgEMasCBdn6XKMD&tyuCRexj=y_7{Or zvuG>0%jZD8F;$8RiAf8rcMu>~pbF(0PDxN>{mK&bVLVvWcQb^AI@hL#ah?qmV!h6x z4F}q{{}haUysRIJ0?NopF9>yx;afk3_<9v-L&ZW8o<#hl1i?GAqU~y<4kag<)^?8d zA2tYDjxb>)v~5nbp7wjGq(##)ogyS;D|OP0Mt!h9Nv?X_I-yVEd>p?NBsrYSy$y;b zYw)2fl+7x^h^TGsRs~bmE_>|7=_f=r_AIZuPhKt&vE9WBTGDXYr{Wp^2Ld7N; zdkC7rRj9LejtW-TkwzA*r4!_XNJly{((KHU-t23DW(%unP>#);<3dt^JNBlm->2`$ zj6{F#m4WZKZ&BX;*81*E2t^5VXIDDKru0)C0oqe7H_M@7&mfL`V`9PCR|3l7e|L^??P9woO40x31XU>R|IEqAFE0V)2>1voP=DS4 zy&@0=|5xo17j!vfzJnwIPp#t$q1H86S6{GVl)8KN!1(O&lozqih7{iq{#o?2wWyxT zuuAP{g314?Ib_ct&dOM|-1*1k@p%DS!O+3{=W`$rCpQjKL}h@`m#czWs$@)%^~JLbB8* zC3wjvRD%%`LIu;2qC_|Z|4i+l^glT&zOB}m(w4_yh! zyso#|xH4aqH`zyeb894hYuZ@3b{Tm9!f}LVx(wNxtiI)bjuoOKRBp!LZKFvfpHe#< zFBZ$rbsOyhyKv&;J4;v1XIyP0>!=lIwOyI8a^Pn8_SdtO0p2Ct?N*C5G9& zC#h5#UgtptcJk9a+Musada9}I^5&@b%3}!Pc6f7a+e>1K9UtapnCK9FkWR~gOyDz@ zG{U)HvmD2$Pd`3YqV9+etfb^HM4cG^pyY>|$?W7RsW)-Kyruj+-~FlqyZuTRN?*8{ z_6|5~0kprKf+m#u`^LXq{&-3TzzG7NNWe>gZ-Z|EtVjU?ezStez^|crvCeiI=Fc`y zWTN>2b%2pfSLECRYqoDp$>UQf`EX#gsR4R&%hSQ*avcp7`a_#-K=_%xu!nB*zOL}^!wDb>+SJaOy-8;@Lu6OG9@3%ig%lJngitl)?a5= zIs$BZjt$84vQRB?$}1U_IxvgybQVj> z4La{RDBM!AHTWBDXbQPIkvm;xsdvdKb}tUKU}!J#v*N+ZFVzTeRh$IxI>TaF3aiAx zSF%-ge;VTXeBNGpq9HGK#5*5SKK35sCE5r8(=|4rKIT+n5FYmm)-N+sRQK1a^D>uO z_v0mW7&Y`gzcnqU<0f%l^nQn^46$G2{dB_zlxSYq8{w7KewoZZA)Mv*Dlhqz%GZfz z2ffpVcCeB}zaV*JWKUB0N79zZamH%_2F)FutGu;Cl@UVL2PbqjItrUg4eNI0>-HLG z;yQKKQ$ss;`BEoTEcppL4{LD)w3v%xcmWjJu>QIuNX|a#Pp&I?&ko3XCe`$dsAenH zYK&CKO`I=X**R;*cl(i4nd`MZ>YIJlI-$0edZ|zUA6C$Xa{dLA-*f)3;{E#ikN{`_ zU}YQdm0$05WB$X+blT7si8o^z=%rDOwzo)x8ARE&XgV&Fi-0YAueGg@} zSV(EV#vHQdWP+C4>d2U6%FkJ%nti`#FcVy1Z0f@jyHT@`^`p0Hk=WxQchQmQuJikD zaIzV1mRTnh7Z2gqF^{s;TA8;%GyHQEZ<`feF_tTUD}C?k$~=}LNEKPnT{BHoUcpuY zU3l9LLMQ(SUv>3M9w%c=t~2^kEAvZoyza?RZu0G&jo?f4ddAX{-2yVId26wc-F=L7 z44zXOfnIKMCWYZ3`G%<@`t9j(<-JuTdq@8K)yd=vkb-c(WT_25|XnBCr#*~OFeDOV-LsQjR${j8&&Qvq~2vL`bxzV&~zl|JJ=mIHotru zT|<;@@?R7A(ta2HE~Oom$-@9ax${@$YeOf%wpAiaN? z!4#_bKQjZI3wR&`xVZ&bDFnO)I7b8x`iB`tTW1t%abKwjVFrA67q_sDVJ=uwZR(?P zDWR-NCRY(%v=j=MzM7rb7G4Cf_aB&8G{c<5x>g!Ar-|Q3sH^whlby3?Vp*PSYm=B< z+&NS*R5%EH74vHX7tF#^fllmNJTuq211 zOJMh$$+LRgiwV|Op(2T_D|5rUPV*^zLnkvk@Qre+&~@~LVT+Td@8Oj57HJMi&Pgip zN40@%WAe~x#K-L10akA!VvUkOaR?=3Hdh+$7nhmolR-S?rJ`NGNSx3?nS<}m8M)4| zrb2@2=);cb787h!Yfk8ZsjKAswqvLmIHZZ;M75LP5C>ch&$+w_t>~FH9EOdrcnm<9 zeO!Hif!Rm@W2-$bv=s*1Cv%zE(;xn;`G*m$KC&B21}^r|I7^VJoMcc&+$%qo)r{Wij~1D3Abd(| zm$j~qp@G*jSwO23>L-_}N74F27Ch{PKEIbxrov<@ z^AG?*Wur?MNv_gf2}xMS(h|RyV*C&{{5W?3qnZwBJc{HKa_u%PtO{dpju5dFs0%{; zEi>P$i=!+=K9k=3c_y_>u%@>wTXN$p^9IT<4Z?;+1-XBdZnQ>DcWn_wKEd+bK_9BY z&xzxJSu_%+KZ`bn-&m%%U}Cty>2e}gY}okmy@pA9!9ee^UoeFw-M_~F)c&kvYW}A-ea{A%fdXoiI|vBOuf4fn z`?oXMIM8$0P6WVnQ@f0jF`ZdkXhUSi8PQdUCUF-X`UESBb;kDo)%cH7;A}oxfeuBz z<5Nm2X-I;|q`UuT(##b-q!%(|PRq$O!fGv{0;r){bzB7(fj}gm?zVo>dVK>Ku`wg! zjzR3qauOsex=hA*Ljkp9!Tlq-XWyE+!hE!BALLofo%p>klpFXidk_nT;94YKTMNmt zTE%_i@LFT$Z5=Ej7=t#GVm_~C5zy;+Hv%^CMhAbGWzz#_} zl%qRXFz3wN{u?;i=48A)Xoj@`aOpq z+$9u}?!dxmi_sSf^ZyuN!BtzOntrMomsLhUw!9D8cs?E_eivTRzIr(8!HcCl$4IH~ z;~FAv`a$)8RTrKuys0PN`qFRmDdgdFW{-z+Ykshes+~LI2Az-`&)m&*H?1Uu?OEXh z)$(eRh*8hdnLI^&omG!dB`}l0T+{KaJm=91m0p+=#Y>r?{{xl)r$L)CVK` zTj5`9{N)8lX!xI4`8DYeFAG3j$^)2r0{kT$P$&=&*n9cC@3PE8W3h^@D@5QAdH}U1 z1ZBO0V~m_iP|9)(t4w30G1j&th|ahuMRbKGhDzU`6hT2B!;O?5dBeRmdFUrsq9H7p zrCDZo+^ILLo@UJA+|Lk6JNv2FdnryRpeAQ4^hg8LM}H-rNC!2Y3KOq_)qwiPGFg4h zWZ@8evnLa#T!dQ&_dd2wahkiWK)sA6+i@#Nql0@aB>G%tDwFErtpE zNmo=I*z?3NyM|45oNw^jA9y&U7qFthLhIbZ4@fGjgIg-}KsW-52;!j@s}@$d>Rqvy zQP4PR7wKJ;KZ_=!5#D=Snz+sz9+NlSNnCsY4-Ab6K0%n8Z{3%0zy8X}?7C@!SJ_+H z9T6+0*+=F=R9$LHfczx#kRz>IR4ilC=dxZX-!e5eQL!tcFu{6s72p#l_=|x5i3O*x%}Z~4`Aj65EH;m9N;s+`y%LRAi%G4Dl1RA)T7u+(RuzL=NhmYOhR(|PL>{D zHQ3lk8d37zu|~uz&$oueBc4Q)Pco4&D!9mcT@~byJ2uhGI2kL*Anyvzg7qt?CP=@=M@SeW9^<8#YW{g0L2bEi8=!=VBHPiAiKa z)1yn0jdB@PvSY=q*4q08yBuAIRwOW%T*_@}iX9vqRwNiBUpc}uW!zNt!T1NrP0l7x zA{yj!%u9=Yr)12@l6VXcshHpiwB3=c2)jqUeZi}s*+U|<tt^F(#cGIA z7QK=i;7Ml{dz z*8QSE_PtN2EyzKrq~e(sB7G9>@p5wc?wcY-db zKHIl*gJPm~{wg1Ho+d@|jn5@J`{8|h*fht+SGa?sEs+4k5p!70)a_kyDj|@HbIG5u z#A^_8L+Jt)b%|-LRElbN_$Cn_Xt09r!WJ#^n?RW!c$gVJ$P)iTo%IoVr@cS9aaX`R3B&baw0w%!Ctf#KmCbeWho&MkH)9CB9p zuy)+P=L`ZuK1060jW-Ls-xy#-K9m?p=Ey~XMXMCEK9X8PT@EqY59RMbf#$SZjYs=6 zsigvh(myRGsUdP<4B}XBoFdQeg=JjD)kD;@kz*j(X2fLmxW^*@&N75&*bhlZIJRmR*fBnf_S&D=yO#c!mD@sDYr^&$70)m z>vqve9E-86_L~!2Y}a>oTEkT=0#Vj7pMoxKA+zbFF`snX*!(04tLwM7C>O6)EWqm% zfrtn2W;vBC^S%aue!j8kP#N0YSiO{)v0rTbTw%`sVB%vtmSX5ZC`cU>2`T5?VckH0 zZWWs%h*+NSupkck=F>YTbb$qyRyC3o-tKT5v+`x+D$V)#B9^$i_}iY>$vMS5A+kmz z^2Hpi66ScN*~m)ovqj4zY4kXV4`*qkjy$uL+3# zV6BK9H__wL>%PobTH#cCWm78ZmVZvtL{TJ=_Ws_ZJ=T|ua$r|m{L3W)$6@m2r5Fbd zr)2f%2Zl=_3D-^g0Zt>Ao8(b3Y!tpnT82fom9(}-Yph~R`^a_U(r}0N@ac!aI`Xi)^49L1EiPguvI}}kZm7wJ!mN<(mW4K-y!4juy95PT0P7bKq8BK|4de zM=g>pr$+rO-Md!oP=jGySh)d6s<>sYIe*MZ(b*M__wMp)Jz%1$wW{4!k&;|`Ymhs+ z9Xa?h3dk{h^tP2~aPWdF)HJwdnJpD441{?Y*DBo-hRcR&WM>0`(eii4^$kUP@GW6& zui{tp8nb3oiOZUTbjtFytP;_*uWv{{%Dy|ux__QGhds2$-&>z<`fQ=(w~Y%Hx~u%T zQbe29>sxbm_zO81m+6gttFq6@-em}ztiknmloQ`RglHsce8okEDdA$X5l;EFSYmW@ zaSB}k=~Nr8bxi77$a zoQr)*y`r^m)uW#gHeRKol$n58&88%MIix_uuUCKF@V;Q(Ci4A(*juFfAkwl7e(;$L z3E1dnm%@?vj^@6+%u!Z1NMNy4NAjt$#3D_3+t5ykDfmzwyGrs(I94y-%_qmOx)Hc$ zgg!!o*)4-0g$l>^9Tt}-+1IcQRAT-M_`CD&qs$tV#)>z+-g7AM6R7W=61%l+Kg6@S ztSU2wgTL^t1^#sy9lvaJ;P(X1~%>c%nV4oMQTk?k!Q~oe?8%wkFjW8fIHK zsuBpi3h+Fmvp%t7aHsc0&|FeIsH=2c^ZC{QXLohwRR^y9ma>wPem<>dd7DdpAC7nW zTzhaZljVj&u{z>yfkDZD$1;8*aG2kWlSdEzxtRn%&D;4#D~4Xvc){vR`E3Oq^JqA{ z*n4q$^7{P~xG2Jh?yu6Iw&TrytL#w^)Ynj)V{1#UPWp9p+15k(VQFl&qapEf?sxXS zmZEQ(8SOr6=nVJ;4kj%_Eu3lI~-`pP#;LKM6VlArorWs7hLD1|4qkV zW{89y{soucbN+Y@6;PWR0cO$xRuTbbc!8b}nDJ{bcty=PhZp93_X4OeWvmrR#rq4f z#Y~aX8NyoR(a!O)0sIIg-FSP;opSs(F;WHVCDD^o8De{Z>Et6vS?x`~6*ejLph{Es zgsCSP!?z#P<6W7lsDgOJdaG)9`r>{tch?LBvOq$|S#COg0ei68=DOtqKVI0NOUns> z-2QngP)QAE5wR@bu_mt}H^iHyfPUE>+dZeFPLAOs-K7#FHZ~+D5GP!k^Q=%DvHFd$ z9X*}z8wC5#9oK9R+6dx?5f(jY__bN6D(epju{o$W=rbUVt}8qfKl&m(Uh*-#G84~M zw>8(^ZKZAVs@;KGPr7}uP8ucOfns<7qX+9v(>v-c^`l(b&ylQC%4*wf6wzmi=&hJ}jb>6Z>xf$Ft zMjHeFF`=J;{cVElmhg=2>ogM%2HN0awRnnm1etos7+Xcb4oXozDT)SYOH)x2AIZ*X z*>N4tn74Q!-Q{gGB0U8~HLf-RJZ>+?C%HWYRy)1-s)NQxpqPzT703_whd8mtOga-% z?N;qNWt!2dgf2xrht^b|yQ_qDWVH*c;nti}W@ig3X$h4ojrN;ZZod2fV&E?~#J?ZI zf8F~3Z-2Z76yU}K&?kVGcz_#5pf3Wx{%fCRp@+Sz3r5VK0a^620m2|>9gb8DUP|>u z#vHNUk<`!aPF2hYmEa0q=txGH>f!avzMj(oGQMQ`^Ov~BPdqmknS$dDgye|CO?(My zn1(pslOVJt7P-^~#-U5ka{aDzg2kvTrNOx9kOYjR4A-i=E(YPak-Ld~j{c^%nhBM{ z>;j&l2ngAQXyFIY7Ch2uE=zP2rl8m!MtRZ%sa^;rYzt)q>*i!ETc0U>v=0wvs31}3 zC4CG_7LiY4WfGj?IaHQ+=TdO+mmBvf@y>U8ux!z{xK0yukdcfVg36oNqGtSu%J*#O)624PoNzz1k7IAQp%^)~)j1`BrLpbPQBG1OGTTS8B8G!Y8_~8o<91#%Xs!hN!ibMp z3T8+%@RnzqAqQWN(OVhD)k^)fR8_Zkq|F-rVezx+XPOR{)rV5iT6ZbNmiKG0y=MM_ zz2}|XHIsN22={m@%#7_FV{%@8VmNUj7bX+$((GCUf5GvW84_WLU*msj|6%6)GEkq| z0A@mfI0*nV!1?Z9YvSMaX*tJJ#T?sinX_JKOt+#uVs}7qJ_1_fWt21?GWQ`YlUEWy zfqFA7={C7k&O~we94kJr5bT18s}Mzi88_H=3qzV_Sb;wJ+H>@rMY#l@v-i#MrMqyi z1VvIay0uMX_x!2TJ$lY)z&DlMdGtwK>IAek9yOUP4g#i}Jdo$}iUD_L57-*usM(St zLo%wUOSL67QsyOTE85|GKY95_*=C&Pe6qCrr`@&LHdp zxn$bx5DlpWwfM(uLw=7eMfWa>$>n>;u1`6&%5I1g^&dYQ!t7UR)B1u;3>4>nj=3v{ zZp@;!fVrPcy6hsb;&&?80p*mrvmum{pJB{36Pc_XyM_ zD_#o?yGWkT34c9LHjP0@jxuLDWFEE<()A0N`K;-!_!lw|{Q~x5SsFe6)x$ZopKA!L zAFbFdvOhues-No2`aKekR;x?KD~z`t9Dj8Y3F0M=**gQJ1kUHzqihnarM7h($uaj} zNe_|tuYpy=FBn!P7@2KJ*<1C%(ibcG|3>65H>AS9GX8JvkJrco-24OtyhniU2Sg5V z;|7EP>J;dP65+8Vwk^GCO)qWZ9X&nmM&u8AH47(^T2xFH#yK(Sn~L6AkMm$P^6E{u z?a^ycy_bejcnQU6&#Lgpq5|V1X^AM$Lz6RXx87s3ap*+#xGnd<+^@9N&{RugB0;t& zV1~KtQ^}t&6@wz!BrZN-AP+r@T*Mi~RkMGqa3;1IlpTaHvq(~He1Ehx0nbNUu7?mT zzds+>-HHWv%|CH0wtXX#_Q?A5RW?20Zhvh>z%diFcZR>D@tPouPV^|HpmxphJGy3^ z=bM&^ryRtVq=Y^-D;+}NcIx`hj5gg80uFtQD1Wi-x(3psJi3yYw!#|BLP7TSiIO{I z%X%Saqn9_cOd%zz_*wstvb&6`n`s*bOmQ!6r4;w#?rhxM-JRl2ad$25?oiy_p}3af zTA;YgNuTe1dcO1P^qsZt9kLVdKe^V-Br};Uq%9(nI^T`B0Nvn%hj;px!`Z8B?X>m` z9WCN?wzEa3@FP^`WKCF!tjS`1ACKadCXTCB))1eL3|cTBMMvbzPf}h*)+M!NXBWY& zy5ocyxLSK)YENcFpXYXsnuHUDpIb@7=kzwRT*pUc;$5m>Ho9q0SL63=JIA{$>ur6@ ziYJ91~7^b9It@O{KpOn2=VXW0B8MchZt~k z0@OLcI5&(W;0Bl<0sXtV`#pzdwdUlo2389n&ljw$#LT|v)UQihcdk#l;-^;h3ZneT#|WIrNq6r8NTPX_xGV z1?JFCr50(5My=tUxtr4kKDKQUyie8oVbq*RUS9RMaO^ zwR7D@5Qzg@9>W)`FP@Kh=Jt(eHtI;vj4W*$3SiqPrnB|x$~iezdN1uxIfd=8wc84xA zctzX|lw+0-kxHeNioAq(m0W5;|MGS6pz-)r(^18rqNfKYK{WUi=jlO}C>?z0HyLn} z8l;UKfh?lp7LV-_F0D>!-^oK0EjEcjsA60Qm9CIUpXLMMPo{9h&OyV~EeP{)d*x?% z+Hv>6=$WMM(M`o)Hc|?#y3kmUcgzXgK9Q!Xv%jY8Q#TQ*L~||%25UtkMG_M`xrnj6 z+!qd#R#QX4Crl?{{$Qdh#>VMyScpExtX527`+Ya{#|4lt>6{LQ&R8Q}@T_xe2YPF1d47yai>OAE6R9p-RC;EG%5uS7u^KLcbTg76Hby|IuzLoeZ*!iL0cYLD55qPVrL;7IUJB z`e3XT8jTM}d%wpAXMbQe?sZ8JztX6xazjyrh@X#lT!34;caIna#-1X^+0=)@Pw+_D ziGgeKqzn8pnCR+W%bi@2W70$Tz}qmmH)n6RFl^>`sGCx+6-_9AKR8!#ysoPPFlXyGI2t?}PnjSA;b9 z@QygP`S^2poH-tbM=*Ei4`~R0bnHT`$q#QJH3=++%wv9ZDu)J~z*t0m?W%68Yt`<7 zkI%}$tr+FzOMVmwy%dFqM{-X?PCET#9EaouVK>UfKux`Fp6>Cxb6JMLFfPO9IGXI(He98uxD+it|EHhWcKKhUg#&LG8LiWGI2ewv>>z~s zhLlm8{of;7OshrPE}iv*RV8&Ub}Q(j_0699j0Y@yGmv{I6RcQR-oY{EHzEHg33u}I*;Xh<_;bn(Dgsp2u6u}dGeZi7?UD=T z>O+b^@4bZmT2U_qvtve8I~7mr*a#-|*CO?x*W({_TKbglQx|@Pjv`Cv3PL)u@exD6 zQ(ck*fu|Z|1j0MuXvT^ph6c3j>dR#s26H7@zD|uw6Bk%8h@9{ z`DE0MyL%bUnDZ<7I&DxB{s*v-;G%23ad&*7Ua>B= z&s^38Si)&*ZW0*W_#<-v(=`3DLk`0EZ|&UT0Cq$HJHR*;Yzg2d5C|$D%kQ>oUGTB? z1(BZlS2tM1Y1( z8mO7!+6g{{SiiW+#iN*ZBnJh|(5>wcbCZy_zFSIC#!Y#_7XBx4PRtj!olPJ>qPjdERV&TX?bY4dnmrQ*_ z5#k%o!_;*uvR!HT!hRI3dy+S+YGADIf1emf<1SP3J)5kQt6`cc?ry=%XJww4Neo5w z^RV4>?$b{Nm-h=@Fjhgg(5c=jZ~--LkfPgSY1HOz1d4+beA>1;%=!ZLdrvEuuf`*- zRmX)S^>HWx$ePo=L2OjX5g|G4rigwXzf8tOn9bu*>bhm(MF9e?_+OTP+ba&h4hdlA39vH@xN!l3|37v%eUzE6FgW0+z`)Q~ zHd(qqfzvNIaiwF&6O&h2e9#s$mc!dazHHEpCv$sLb{KG+NkfT8ghrv{@YY#V#r8YV#CY~ zUW^~q&xjHn6N8~;xSm3pDl!J+BEtFFX!hC@pq(whYT^8n{())?22;-Vg2_o0({*Vj=~DOktoG5Xaut>~Ch8XLZDwcGq7zYH-{DHlj&O^t zOj;WjkHn90)o1%)qsZdJ%L)26xLw9)IBpwWajB-2Njd0Os(voyOXMcQHC!(v4wkXk ziC|S3_c&!Daq;azYaCzdyCv~h2-kl$g|%Wa$Uo$l2zB~FSJF386Z%4yS>a;TZEB$( zVN^zLs5`OoBfmak9wa!eKIzsT9&i+e#w`q`r{eoV*pQ6RKA!z+>4btRDrdUFL72Rz zWO?r@0`&m-b+c)E>3th%h!M0l&%|Rj?winuw6&#MgQTZ)<`5kC9><|{r0N1mXPe4d zI&H##u>7feC_(@8#XsNve?IQAXR&)$7NTp%{r`q*7l(d)E4nT|0^>8vn9gS7iR*gYXG=!zeTKi*hRkLGggMudJ!d%%&^ z16z(yc>vu@4_NrWA=&2rt`LvGze%LKPUHfShX2gQPr@h@)~(_J(yXw!kiMS1n!bIuIqDJ4NnIiPV_SJ#3Dj-!@7%s8MX_r;&>+Y?*P%EXsHd z=f6X#BT#5mn8hs>o6M4i8TBV{Z8Gn)f!=l1j3?!rRZ`U#>AFixGY+ZoVZyezgWNm< z+v(hihYoM9UZ8rHAsQ!ZN$AX;y{Rp_SO~jo&hw;ij$3X`cIo;uCTP1=k;UVO3!@P{ z(}>~bm_1)0{;R#8!F|j1ZhG2#JbN?8p7BgELP*Df3r}=TpZq0kW1T1jQ4^H!dQxgP zd7j9O5sY}k6d^oF6gJKY?_Cy4p)kmGlsXY=*MHywim?m{0-t$IjG zc=;yav*Wnsj-jEVYfz}Mj=YSb8QBopF-BT;9HIJG;#_;Qq?~a;#4p1D&WSO6FR+P{FmMyyUn@;Qj2k|-cjdj{XD8VxO>*03M zmJc7YA?aZbu#Cre+( zu}j2#wCgzx^0XK*s#@5Z?9}?DWl0=m#DPVx)S~Y-FStHiB~Xlc!)<2;hURv5qeE`W zm8h(c_x3fTl~g0-RrlkHhlwZzP9L;%8!sYSwI6>*%CY z5pO5SOOfWhmrM2!8!EnBnQouCUtDa>qCT{S}PTk!{{~;d7gdL^q>=RrI!Npw@h3XmjO%6}<<@)U{ z`{@q0FuOKU|55{cwp6?+SEM^(jM1jS)??A7Gt<1>3uEpW)Ia zU|Lxxu>!LmJ-&ZBFqjaC;JF1u=Zqr$o^wrihnXlE{gDF4hf}0zY5*C22$LBZbCr-R zLOXGz6hBOjT>&4KNl}U`iy^?KNG681!j3oU+b5wC8$q-m{U|C*z7KtzDMM92!}Z4w zEr{X2wR7_Yuww|+KlOl}ZXg%1pa?7;kmdJt<(*~-XCZug@flyR@Kj0XGAkV3=EQ+W zxGz)9>23MKndtdzK6SVu-<^WLyi~duMe3k$8Vh;EI8Y^p4`4Ps?|)W_aW0zkpD@lJ zXDUHNo27nAD-<`vT;EU`N;k7^G^nI}#ULA%v5E~6j> zNdZRt(MY^* z@~24WH5x)3JSSAOaPGXXI#)cHMUL8t{t8b&tx?MqgI=~dq5^&nKzL)%Oz7dga_WWA zt9ihI(JC7244sSL&fu9R6TUfWvexKu?Zhm2v}B!Xnl{&cL8fi(y(*=osS-0ZGc0{u zcgfQKWcfM3$)fAi{VgXq3!^D?+0$$a5BO7fw|uj+0dq=ry`g3|a#B!+lxOwu9v!MD z$9yH{aN-ty>``dU4*MBCnaMBbFh+^KKEDyzsH?EUKKplrs%#So$%%LNB7UdajWoHg zuT2mqs2ArWi>hBN73~bf8IbW6DU)gXz&RBR`s>zHY2=Fvx#eaB%w`pb1B-9v|CgVS zAeR5y&o|>0}G(wK4VzHtz!Dz}MqV*Z{wdj6Tfg?MkutnB7o4-vW}i z(IxIhSLCee(StPJg?#V-{(>wY#1In1M2w-|iS7oz?U}KK+bCFA9^O2&Av8Gs#xS** z)+8UqDb>i&3TaaIp&#H5-ikXF9s(o3DS9N0;e zfo(qXZk7G*Fb=BRCZkpUuD8jK0uv57BFH2|S(}+vzJfOCp4`sPhnc9L<{~l45P1U^ zah!0v7iWV&1tT}(T^xpQS>$M*{cTu{Cp$9owe9r}|4W(}M64BdgmPIp3cP+n`tyR^<4o zg;m(Llbi7Et!ZIS)iQFj_>D)(Invr2LG(8ii>*ID<9=Syd#6UITBq&pc>^K$yfh$m zs0EVMd%nFplo;X)@2OtI^Qt|xWKDstvikl>y(n~zvvB(8D3BTptTj@|gnwBHYwt@$ zI6rFbxQ~rWC7x)rP&l&EZ7$ub-nnzR-CYq*S9H+WN z$=UgFh;tr^f6jb1M}%qD($^TdPGPbmE;3&sfxTuBKfcFkqZpaLyQ69|4)5A}O+K$` z{ML=GG3?mt-YD{!+II#OHp}@9wlVwE>sV1KK3&BY`Qpz(*R|sN^Z(_C9>nwC`nk~u z>K)+s9IoyGHxWPpb1z_j_nUvu+u3Kq*jvV6+V&hK3V+LTp0!z~jesrOFc8zdRN_&7 zn2Tyx zK1l0T4yElQXcHxTKBBn9kzqxqsi*=|zOlAv(n4dNZO}GVh5#$Mfu=o$#2- z6>T5CFx455?9S<07sw0YqBS!5;}Z8Nd@+b8l!#7x6!951F&lxCjcdy{p+lB>=*ekl zXM%&=(|O1WB5FGJoX_gY8sz(|aW$r<;AAyI2yfn1FD0DcNaG|WGDiYu%QG}+c$P(& zC?Whz@70pF0O5|Ui#HUFS<#+vJR3tq<53C2B0UmQAkp8Q#PrDQIO4SXMx`(R-gn+{ z7MJ59gMp_Jb2B+m#n1|;Z4`P}2R`od07|IZ0$BY6){=Tj<;m?UFRTD?x&PEb44~h; ze^>s!Z$q` z7tUHf+LPDDLq;HWs{|F;nl=^{ql=;00ILKouo7(HyX}Lv%(v>q?JCU3`~HOvakRzK zxhey*qK&=6A@TOp5YAihW%0BQmx|-C{06lUhGUPkyF z6l{N0s^_A^I4yqrfl>-t$oULK*}@e=n`3Z3OR^!Pi%B>J0@kVh@nME`jCV~7zX z{C9wWv;H-t2()1vfT6d5#Y@0X3J@~D@!##({^vkQM+{ca`!0d-EGxTj6b;#XIKmsW zNmX~XDG*lq*kKwHdFl|*%8?N((bJ*l;7F1mL@DmRlq*Zkp(}G3h7^i9t@^)+zA+R} zG~?x#;9P!NJ40_ppo6xTsT3;94|W$ftTiW~?H(%A_;x@q%PbE4?Oj6l65 z8|A^B@W~SG%+06}Og>-mHyjrsgP1UopwNy`ehUfHfCLBH<WES(TNVRB%95Wz#ftl1}@eh$gCEZ(@k!9AlBtUH3T`dbnT!$r*7F}=F(|yb5wB8RGQGNzGdV3v+(*6BCsV>Fjhg?y zc=WjSf-jEc!XC3#RKkZQ!qIOk0Wk@S(iIWb(phDxGDRcm$Mrhj7(2RbwYG7}LwhKb za}LM_Y4lvBNt@ClYmwyFPt`E`==RJ9kLE%5TUos>QffP;3^!B+gFaA18td8`!O`rm z=u8-dPF=llLdl#DwFEcOJnFgC!t&NsPEI-L;mfNEjow(q*O%g96?tFpE z{l^d!NF3PxyYe4H*TB56Yv4WMD<2?_2jIsK2;g15|1l(rlUjhm_4fPYMtF-+@j6+& zaA9+QUMYeQY$0>^tpVEAYjUrl@<)>83Oqf4n z+Lv3Z5cQUQu==$zQ2Z#r(EL`Hml3k#tD(b+@Ean9HmO|OweiCdlRH(th@zB=;Fh9m z)g=YRPHQqpHbG}dJgK#S7qYTf8MmDlnjc`+di%IE)k7;Xxw+4`EW3Ik?eE0=s10Lo z#TMCk-pspl7KD2ysBBf(l!?oS{k<)zRRg;_u<%43u zdgCYiz01mcXI@msKa!c)|A$AQ^;dMe5@y9sl_WvNrztHTnp3DF;{#Jf3k;^|Csok1 z&+xL%YKlX0Ljp(!V|`LBR0#E_e<2(nH}O(+*LS}OSNWn}^3Rc@CT*K(%XpcjeMUv& zwJmk1l&y=eRm|MMk_m=~^NFdsjMGATVj<@%)`GZcKEa0(aBxrZ{GRv`C{9A zCgC{M2QK#?L(HJxBM*O9{xt;55xa^70tv{Y06a@v0StKoDSo$Qd&%)PJNQJg-)4Nn zFvkSGeLEeoP7~ov(Sj2%CQ|rihVGZ|W|+*(Nmv`os$l2MD7JQ2tx6yCjSq6f0pvwz z@%4NmI%{4raB}$D!PQsqgY7M6brkX)RB9atx!# zvczj;0<&@%_`aAiDumLQH}smOW{92UPM1rK$!y#eaU!8?J6kB^#d_iFC}9biDfTMbeLH zzzuZKyV@df!=$5qlGaEXx0UF~(PBa-_!JSlQKE)5Pu&&$=CD58qgnnd+NUP>o-rC~ zVe4yN_1W7-e{K&1zjknXF-0Cz{VS50r-Xsht9a=>IgKbhI`8#{-cibWgQPRWgp;8t z%@uUGM*UJ+)vy}NPe^4%at|l)SRa%!U`$$2xIYim*}qv#%l8@c=(yY;W;*#sIgr2cSSDEV_~Ii8j#BAbt0DXh%u@0 zO5fBrV`ph+bk1)ue+;pJfV=y@EdLq;=51Wzy}n)o`e1>7g#*SiV1JLH{%*?-O9t|# zkqPEKQNV~_!LKg7mjX<(DufkYo#4i-zik6A^TtmJi}+pbhlxzQpH-FCi&FS0H82EY zR41HQz&~c5*Ul`*p9MyXn!3kA^HPbJQoQeMqk@}6^vTJfL6eVJ%pFciDU=-?38WL- z6a5xzb4rvdR`Azuef%We&RPnMjSYk{4>RA0hb^fYQW`aC?=5;hjaivV__d>N_{n=G4-vkSlZQilzzk$CD(qsxSVP zdYl9qU8nQ>@D)~$gFPJ9&F#X#Fk_lfJovpTt#6yG#9Bf#T7wOhkcsu;zDfYZrC`7424c2bj|xF@rdi#18Ch zQ&GOKI!L35AF#8MaGX53#PRfWP;v*si|~ExwtQeCF)Y@r>jfe0~_-F6-`cC>t0GyD)n zN@c}18CLkf_pD)eYK~|%v|)H%?fdgM!&}Of&+B7g`wpM+gP#V*r#V{p18<}x+%+rm zDSEv4e}HknIHqMkv*E=u`TuaZ%*+rtH(d%kE85}wXja>twp7EpUCw@lCms(`98HU}UAYmq?HQP#Ee?eon8GCE<&=YxwT%8{3}9q|`19?IjO1%I`#jUX-6|43(EGh1pgQl`P^b$XW3D(V9d`Pt?dh z>hL@xnu$C*H0&Gi-JdCo*=ZJoTjkyeqLn8>$b`X~@2RxuRbiS4z>&A0+`1x3MP4ut zZA)%4sm`MoSu=}&{Kf7t6`mHPrW?Lv{(N`@b7A$E*`R;o-(;_^ZCEU1dqLIQaUY1% z1bH%A8k?gT@64d*=}=x|v@~UWHX!qbVUJ>{FbBBge;l!a)PdcOH9`b(|prwJRAe)dG?WcEi{{Qx9&w z|EdrQu42wolG1^GrW)qEXaF0AUOzmbV7#TQ0r$(ml-q*NV+2OrwzKnnNFl>et$3Zd z@9_7A^gBuhx_6hAN;vn;Vv-%#t=PIlWOAs%BCZcmh%U4qy;W5}}${@?8 z;?28}W~T#=ZECax+A1h@3M|WyFcry0u-)rZ#e9wmjHnU3`XStjOjc>Jr$&HlVRv==)274mWRAG6TZf8xwU$RdEb z-<{HPlQmCoKox}w8kZ>>;~OMPUFI{~&}MJHTQaN9-;kW0dfKNROejy3g-a&V4B`{_ zESxUG0{2D&l7j2btf#m{*iJ?Gxe1p(ic%Com}9vI-)EiN{zmi15jzM_@R#LZM}2@J zS|EV&RiNhv24?^S9FXG!;OLZ)&7}~Yz>9OF+wlFy`6n*{!?2YK^SFDBX2UyE?vB@gT%R|Q>6r_TR2ozVO*(fZd9W{KN(HI8ONZWsT z3C%wTmenePg8ZeS%loMXt?Yhw+pCLDuq)_sGw{UtTE8!Nu(AI{xRWLdW9 zG3$R6+iqUa-V;TkzPap@v4B}sR<($^)aIT~)g8cjT__chc$xM=DJYkMvHr!lL)V_< z!tJWOR>Ph0R?QiyJBGcQ<=f%AfTIj8S-zYnH#=nb`l5{O!?oc+c@9>!U~xD3<|^9( zO{EEArTGi@7i|{Z(n4Sfqf|TxjMbVBsn+@8a8Ll^+a5?^n7PYZZz9ch&TJP1aUMv$ z6_0{<`@-tjA5%ro465JhVs5#Wf#Mz(0$f2b2H4f4vcb)jVyHh?DZz(Ve!n#FzbtWp z^#A$dZ$AEQ+kmx@FRuYh$bh9Lz>zz!T?7*R?j!vS{4vM@kGuLFTRyLEVLHU!3^N;6l9no5Y3>3rgX9;aA^Fx~4oJmcv4LGQ2{2-ao}-yzj&y4uWI) zdYcHLC*mB&-pfEay3usklG2adduyjbaSBZf`J-3VxoI$Yu0i|c<5y*C!*x~181m-; zOsh>M=b{aYn&F~sZqQpZ{uUTT-`#n!xm2Gp@vL?W->9KVJ5=Of(y9d4Spj*F!FPOj zuod7`3|uHWz@|cDd&T=y88BLe?n;Q))N!_)k*nd(g~BjZ`t6q&^_cUom1}ZMJG+1@ z&J)JteZh01z=RyenEj9;hwo2M>9G@$yE>Q#zYW2#u=BWW%apLn1fe@LgCk@7~Sdm46-m?jMB$j;a7ZYJekVU<>^HPn~oUJ@z#lkzt9lY9p*O zo+3L`DbE>hAkPPuNs8;_SO#gIOq1)FA?;NFoW6%5xlY54RZyH> z{Iu2HPdv@tzJfdL!P#OeH-&9_pQkrfbJ5M2v}TpGapqPaQso={Tmf$9s>`O;zH9u+ zeL~>sy3FqP;Znuo#&Dj3H1W${NMeQN+Q^D#31YY?C;8-7Z?gQg!K?TV@==qRDQsLi zzjbFUS6(Pgu9YnRe0+vp;dM2w*8YY9pXe=#_GFFRl15lt7`gsAOC%ASm(d0O>8S6= zH@tRJ@i@e&8aoP;^;t3aTAOJW2m`_>-M8hCiEiT#54K+Bw`RuEs`rNX(M+?%lGYTT zBNREG?=2_{bJ)KWcU@#o32o(lf9Ji!c4~kSTy=3c(H-&>Jh37`;v^eRs?j8J!aGz~ ze&qm5@+pLp^6qJoMQi<=(yjx^`Cc7nKqSocWbO&F;|VQ(MV+Xt%s$o6xc6v{hxjZ? z#3d@vg(TwEe)Zjc>dQ??QBuovyIR3i)y~j(2Du4LGeI5R4AjQ`{83}H!bfLXttqu^ z3Z+$13zO==!g+;8Mp3Co@9cG7oz-4f!blxcOs=3uOLxkbtxPZQl+6Oa4IQ0ES`#OL0Ij0muL8Z@Uq< z1tPz1KqLes&d-cN2`pUlM68exjH4PCK8La;EYKj3X%!5p;QXcX{woc}7nF~Nrx6Oy zM)6jVl_H)C#6)TDDe5UNjmjt*fZtboFLBNEK~=cqD=s+O`oQ zf2Wb{p`6Y;weP#xZ>4N_9bFR6ZS3C-Ppa}!(!e>}foBO`#X!1_YoqvE|``XV`9CWP_Ir*J|cnCjxx{oO!F3}UuJ`i~5|Aa8M#nCz) zl$#lhy_`#U0LNOTFKxeql=xogx#{(Z&(?NFo@v8|T#LH987J0JsU-l4+TTjnPqTi3 z528xKIFsUBx;ow-&1fBOOiNhhPI`>nTAIa|%vPmru>3+@DOL`h&l>%6 zqs=PVdwHiebZJ?1GyQak5ZeQeLs`KPB@4&;sq+aY3YVY+G8}5 z>v#E(P6*_US`{Q`dwUi8LY>%RM4*_ z>EqLxW=b9yQ=t1NL{NS4j(-zg%ukmM;FAAw#0|0rcK@#YTPKABhJg1|FX8{k4;=9B z25cXY;`cnp^9~7tb$C1{k9si76t3Xv#HYYqWtQG|y!7be@d|m20uNJ>fo(4OTgoW9 zaiq!;saQ`kRd(QD0_bUA_R&_ti@%gP%ztut>B7{j(SI3`@BBf;rW#s8g5lI8ak`ec zb0xtZ$TB)061b)*G+i2fQf~w-gbltcMqRJw{2o^ZCcIWI{)-AYgo>@BNI{AC703?m0{1}Dqwh>wJQ(Y=5C;5b3`*hcG@txqOLpugKm?ARqY_;)`E67;_MX|rZ? zrwY-<$#`R61bM<++@bI7$5MhetNY?>f6^1<&x@R9BrOO+#D%UEg58-}oL3Xvemy3> zd*{=|{Ekig(U16E8o7cCVg9eIPxBa5U%-(GINni=JMD2vd~o+HAa(>>af|(CKonn1 zXrS^qkL@G=U7A-lpKFZ3vnWU<;)SSn`0xa?x|YA#xC+*y6rqL-FQcMGNv2_{=d?f| zu-*VH9o-kL2Wz`Vk~J&Ui#S65$I2Xic@$9y6FbTs40(`QQpIH@>qkz6ae*?8F@$}$ zoNkjb4rC^*bkJ5cet#rsmkGM>Q)J!>M6eHb_JRAes{AvF!18YAKX(3j;sM$H9VXzc ze?3{fzFxQh(E)hs0%8-0THyHazS4yzMG6@jZb#guZ=Qg+s9G3g@h7HwnOuHiqZn=SmUdk+fo?Qe?fN%xVY-8yaKBJHJj} z1Yg!AqnN1$uR|fva69DQX2=y^*#(8WT6wXo@Xjt~`8Ud4#|Pzs>@lC%AMmQ>#n-0a z_y#?nq(MDZ^WaGN&oPvEm>VD}iiI>|shSG+d>S_q81HsKe|CT*Sf$=Y&TFxlup1EyifK+6#tImB#u-@*R2?U#WOoWev2K2J*N?VU=npYoFV5+Cc@6yDCNBkZ{t-qXsoN@Q-pD#071X5;dqTL+(d z2wk>c3k*m;);i1)$S`FG60KN59Y+1EXM3Ycxxc6B^(JA;fD2<58gCgSp|a$8Q|%Q} z+(f7*vV_CSEA=;)S5?pyeEGve| zK|omj*{G1^?H9UptNJc()QJPtS^?c!YPH^0C(S^8Vb&3AiRNgi42O@H%(aykcBR)r z-Y!feXsdfj4*2+QJuwdN`Dp6i&jTP~}e!P$^B75X^XwX?Gzm z-^OQTgL|FmMlIwK>tIN5!qg@OjIqHQAqY<8^zC6^O$kK3=ge}t zA_7TqUP2R&Nswtlm7LT33vhn+B@K3H~aw*zxz$l z!13Sp)328 z?@$$w&V|m&FvE0l2`fuAN}2Vt2+6bzF?F2NezktM0rjNM@nhiYe8oCT*oIOqtk26% zjW*L{xT7Vl3m3pN|Iz_Zq{3Q`5a$6B#?foaN-q2wui6GzjQZ%sZ%r#Bx4IjBK82x~ z>QUfpg~k_^Jbjm+RpXCjxi=v?>ZWc+(tm^4I}{?Lvepy&SY6*nLJrO6yd4y5MzI#V{;Dwe=*>h3vGhuHcA*iss#@6MBYJUIEchOH4n~7 zI#_4)knWgF5gKo0mcK_k3D;slMTKnAO(jnzinlHGwum-^aY zOHGf4V;a$*8lWu1m{z-Qn4J^?t;CzXLqr%z_Gnjq;&xNI)5 zDCpd6Do-g(FJMD>g`)D_vP@v5#|FI&eJhXe>o8*^3TkQIz09K-n{vtOEd8MZxcW@} zFPt_0Auk7Ns+>^|6QivXp%i&jGpwJac2Uzwnl~tSJ3BY&icURx)v>gzrcw2fq(AK` zYUxbv1xWE1_-4?ReOQCkY4IPm2Dl5tXIF$KNQI@Z_e?l$?FERVjFEbq?gLY#?5dk~ zCZV=t>*uzJ#BcE!?qGr)wZ{EyPVAp$^j?L7KAT+=N>X{}Ai^<-yo6Z=(G(N(aqK&n z4&C8R!HLO|amPh?-j?R(Q~UXyJ@cAL6BcZO7t?AS3Np2Q;Gg~+Z!Tm4JX-v}js!ry ze}@P->t9C~fFocW_hKBVi^2g%zvo4c0>^(lIyGuo`GSWc!Mh)u*Gb2u8i(z`ar#~O z7|n}GM1nvkr(ZD))M7xTfbg;=Caf5Oi)p-o5irf=GdPG%PM3Q_^wFB@k@wz2}yK>4#mfMxrL-l2Zs`opb1`QDxqwOCaXZ@K)O!UK}5Vx zvu8|$!Lel*oU9hA)$%k!j_nxHg1tpTu_95_XA7U?}??M*g>j)7ZQa1(&kunq`*FI3m9zeI& z8tA9!J@Hz!k-`(+;1-ft6P=JRW)IMiIevdg?z})^IxEa&%+|ya)?L{YjLh&p_2EA? zGbR~5GWzIjJ47_ZUxX(56P+)>iLL`}xpHzTWqnGm09I$J_M7RzbTTLE5v6`bks=Sw zNSlScqe-p^SBVpZd-P^Z>W3h=1uOlVI&Rj28u_biO&wi&)t_ai;Kb-1Q1fArc7>bQ@kA-%=nhooNj&?K9^DEGxj z?-+eUBkfVj_#5?+GQvvGGf^@m!)mm z{9IAd(1S9_Cwuq;qg6f{3s&7^0X4vN${mhYppcgQ#)0V#dqOU9cAWY5c`KcRz@*g6 zv0CqU>E7b4O-A&xppPmTJFAaeg8kaR3UX88>oecfex7K2OtI<&v!MNT9iqn0S~@R? z%6UZ7R-jU;lHe+>1Q?P4`buNK7Wgf| ze=Hp@(8@J1;T1(f_-)m)Ut?BGonO&IilRZ5%n5%y*ODDqR!`^#WBDwQMrt`fN|cKs z^VZZ_!Dv`Ef{45hK8VxV|I-PPCM0I1tGry^w}Z`@Dh~4{R)pYc z4#;*D3--opWXySL^T3K>I4fmn6VecAZHV6O5=dfM?lR$rxKuhdzFb*ay?;)`0t;m! zWhB`ix@Ae6Guad+Js7OkKbSh?+Ki!mzgXzf+g8eHr+Exv&Hnbgu9cYrd?n#vDC* z^lW1@=86gSE@!K+hnV2s*5PdSC<}MF=+9Myx1ujTy}=IejvokG0=7zG+o|7&tVLpX%au1eq$GVqb zG6}{9@}c}9flv@IQq|;amD1rn=z_Uj!iJ2V{j>s^S#fMzIb-%8sbBq$pz5c48oww# zEJxKE9YtqtiWCslQ(HHJQ6Gij*oG!VT)YIslAp)Le(4OWK{;Ou5+eqB@YPNH=CvVJ zQLp&+=1@%BD=RpfkZ%WA<1sslEYfq^+~mfdYOsE&S>dtELa*N#`l-8GYD&x0JLGtj z3sfL|TA?|-Y)S55L-ouzT_c-wBumi<7yW?y^$N}T*J0uzxL9EfkNS-rhFLY2+*>ip z8Y|OcYzCc%l0H>Zh0M5R4lI-Wi^b9avgjFR>0?jFLU`wR<_JmH#*$#|#O}Jeu^Hd@ zO|4dfF2NMd*rYes{s{6N_Eu-N=S2KkjCvT;L_6-Qrvv_`R=q&_JgXo9VD`V z!NZWsN_47;8EqQ;y-;CRtE>J_U*KYF+;**h1%&2$jq3-)hF!2=pYNsX=2U{^F^R+p# zd3IZL=N~+O;*zLv=wE936PF%-=iFsLY#ISPWdNRlaysm9mH%AB;a0yHQsU`&D9;+F zL2my{tuYZ0-4LyNG)EErdV#OrDO?d}2twJEjkCpV^CpRf^0as4H4EQYsE^h(nh>6am6gk7IrPEfhSPA6O3r-wvlomFiS2(AFxY1f)$t2 zqE0v9UW?v92)r6+I#m zdcMeSmm3G>a_w{!iZ^vSIW4iDEs0Hcflr5X*fbJBB}~%U@dJ&Ys7Rajk6j{0jR{XB z9~xvvVsOW4yCmKCDAbpi+^}R+ZqBAmR8avDX}kA3wOUb)QfG@P??9WJ-eFc(Iv<*% zGp<)+|4Q`j!gkW-*N9t>XX7nyK%?q##E7rRb>yP=s*p7=$TL6 z*f^=-PR+VtJI7P6$l*c*l$xhX?PyL_RsLzWX1yTanOX*G&_zp6q)Hnjua zer2V1WC6R-sGYN$;Tl6d&$J=zysXQ(m0#d6ms^^4S+UrKhg3&3Z^hF#G=SoD+Us#R)6}vZ`JHB4y_$9 zXzANw@NNO@z-wL;wv)&VKpAq9<5D8VeVfpRYLyl&uDZv!dN)!`ajl z@j_?$rD#;6WZ)v?lSMH0rQV*wPEt~`*j#OwTIZ!LNozD)6m$Pz{ZJy`YXD4lobGo= z>yrGa+Cv#}0v@r38flwdi3D|0NN+ldJk~59;?gpgcb|Q-w<=SJzo2ldTImkH3}lKDtJ+2=4P^fcjqi= zux_;%W7N|-Vi)4pNC>`8GZ=*;Ug4bfsjeq2n4Cj-GS#;TaGBT`OQ)R+c!j7J?`t08 zCj?e(Uwzk3m6%q6h*HRku_CCktaI#_3uKuJ!;^P|$X}Zgm!*)aOK_n-=wT*vL3^c* zv%vb)n`sk8TZ4&osg$g~|3+4UJ|6h*8g(Lq`(Vv)f5A2(EI{_TPH(>x1 zh`;mhWMCFzz|m7qU~vj1X2!B@o-tJXF>W-e)$GyI9yCbmR6LWGE6H%aO-UsP7`aw% zA(>~8dwJ9{O*75&!gxE9s=`OPWF4}74Do|+-Jr}o?qoD1x4=eOhd%^wHKH;nU-g&vJ#8r!R3|4Yok7@8 zwListz!b#ZFL~~iNVG(f@ybq&6`SuXC>eTRn5ZPKvAn-k0ctz?CKM16Ue3uxO}jsN zqzZYfwV=dbA3>UVJ4!}Jc*W1Es=CswM}kRav0fWmf*}plq+mR-)V9@u!+Xb;o@s5i zHaiTu@OaKO@a(#HEc2wO$!v19`AfU8il4~|E`BVLUMQu0 z52m$oO<&prr>g=1W2(PPYT)(DEtQ?jmjs1E+hEqqlb-R7KR9J=s9BTN>io87UK3lR z<00LaWd>rHZBNg8i3j^U!FegEPvb;*YMH0 z60EJCy^QvvJIwq_`E#9zQvs=)A{nzqO??XR_7@@rA%l=*%x)w+Nm)tD zE2>uBZNv^PN0rM?$Wf|JKbpv?()z7k)|^`+bes=Q58TuT+r%xxXQF~n-^b@n^jKgS zA5Xo0e$J;YkWQH3+r3o{5 z&X&$O&2CSRwRQLTHA3N5W(i;aNc8)tZQlk77=!#s8U>f(C=|WKpctPv!(-q5yh39Y zPxy6pUqRZKLh{jf z_QthDLgo{&2%hrBx?I-{P9uuTt${60BK=>Rrb}$`hjMaBKkzYZv>==0uF!LE-awxl za;zhf5@yk$nJ2+t|CcFA;n@G#loa5J5pV_MO|UG0sgJ;U@js>>$uVsvD6ysr;QaEz zPppVMY5lGGT8obOeGBMW=z=ptyLg{KU|Kbre}w0iLaz0AAWQ|I2xya5Fv&y7V1y{7 z_&mE$3&zfJHJ0W_tR^q2+{}~ zHd$Gm=d9h$>BSOCUiCj8)(tvDUL23V={@U8#Y67d4{TNUV)RPaBxr1Wtq%2A4EiS~|IZ}lf@yai#knrUP3J~6_+Cv-3Mi_^T? zLg@r`jAFz6>OKbC9=cQ`o_U$AyP}T>zYYD#p`?U?U8?_b{F_5X0;Xtye7p*%5l{m; zR3C6w1LhzG*1m_Va@$x!Ea>=SMbZy6??OCDKQV@-TIz7`UsKCmi&sxn#)6QpQAi_v zCdbv;Ce8m@Rn<%N$=3#J8_w+U%bEWZWvy|;gd_2|9PeYT~DL>XTe6F3b@Z0+UMJk0g3k&yL4I%4)(Dl^o# zhK)EpS)_|t0LAdpKQrO%q|TaBMu>OxMoVf3*w*ntC#UOM^>ccO56RpyIKfBWVpCme z8=hL97_OFK`B7Ao08Ou$PkLI5F9b>6_=bFP=Ye?ZM^eXA^P;mSYI$@TDKzGJ(=L9) zk>RpviupS#FVXWCB++sER}-W!&BWGhMXphWnjWFg7H5I?tnxE2 z923ozzx%7=l-ihT9?ZfJhi7T=T%?t6C$27TN}_=~TLTZ4RP+oZk23@c%>ViMV@g^$ z@h?UF$)_G#0aIQ;oC3ZG2lUDz0Go%<}x1*)we&O*yiD%{DX zWg`Nue~o~^YSeRVcQ0w8I-%qHbxm1?nMsOLBJh8;fp>N$wSa4a38xscd_wVQIL5Zq zQ#;gHLa0*qwIR7tCT=cVm%j^7NTYC0Ct(~j<0%&JnI0TQ$4sr}gnGSb9~x+|$sshT z9_ffVcxK;y_uhYJNMO&WWAPJ7&+G#uUAdYMd5yexU!q1clftPFrKI{2xMOY>19@VM z9TcOm?xW(RGPxzPuLODJ0zHMEq8CBwg!Ya4OLKGN!qZUka)#U#q$;n@t#It%16&y) zLoaEXRN^NLn*OBe90IioG0hcxVafzsQGrHT8zb-L_v7g#bRf|~YCbGYatQg3WiyQz z_H?=-o{AE|a^N8LIf7@E3=8~zQ~9K~@Cu{nIX z!TsstNAM;)|M2IbUf(m;j71XT8&e8ls$NWO_ufG>dM`;*<2; z>`mk5>-KHS)SN9zi9GgFo9hqs?;}YCiFbE{#m0t>3I0;lMG@d=cEwu&Y zx64{?o7T~l2^KsHOCTaqQ+;hPId>j3stAta76TfbuJJ`=hjfLpeyAJ1BT6Cb!*G5e z-*19SrGy`SE3!Kr0kg52p^!JGC{k|z!8W5PJ)G#F1bXTl!~}*M;jjYV_lq-Me9(nm zc&UaFYSU+_Q&e}}G{+$F8`J8}9{L0RdxZQfpT$uJ?m)KJ)!}NwTnR$ciFl81i2~*> zjVF>1yrv6XH#JF2xl7bZSfG5B!FuKdlzs|Y-1YRHxw;&+gL=1k+9sSxBHX`RwzaYi z3CuCTaoUNLjd8lTPZMArU-nXF6Y!I6wO77nrZ0BuX@Te}xkXk5e{fMx*HcpGK+O+~ zct744`4ikVdj#BJ9QBI8gllE5-GV6!(?m6~0fh=_cI9yLY#|C-A4(dULOM1j!`joh z2Ycz~3bIa#8~11m3l_HLGi?n)?WIibKX(3jk`)G?;Qq_;A5V|iKy3N}9LPZZ>G$;> zu-?J{j&i^IFYn3gp&|+KeLYX}@<~RiX+HfT3k>g6o#LCl67#7`BA-}+uLObWcEUCp zCsrE++ap;VH7RO?EY;P7(pC&ijj}&aB8`Z;vG;)~!@66bkIV8;!z_Ix{goQ7fe2+h z9(e025?YgC+6DWP8mr9UN{X$7^sY%}Z{T&&lGnsl`>satOv-@Gp}I>uSjYc*^i8YQ*Twu zeDlB|EBNBvN+*ytCmSRVSa$P7LfcYzsH<|!cn(TzUl`jvoETe^Du`gXXRvd6+%l*n zRuikCD#gRcBsX z+Fy?Uc)Gv&zxD1$7Ra9fjR-)Tf(M?({HG3e8HN4zJLNlQL6Cs_mocruV9G#bSWZ$5 zzCw}D)Cb6(!TERkcR>As@mCqRWj8-n=0RZaIsGDH>Nj3hHi*rmVqxnqR8~{WTDQ2+(i?y6hk`7`vG3+17 zMC6rt-}M=nzrgaEsgk3spo&Xf{3*~&%Dh%t)X=G6uplfe+l;ANw~2(%8I-)Z zcd^?%gjv*H%@xPS-}o_e(B7r1bfC+R?212gO?9QX`8w#=)f*mh1Yl>GABfs z;%>LI2dEXf%HQh9eo%ix%=DUcAL{1K=?_Zj}X@h zw10}9?Uxu_S&)BoT4n85aXBk^F^}`E7g=K_zAs+G5z!?FcMfZ|%egU*k4Mwj#}reg z(ox90u*1LQ{IMi2{JX^YcjI46s(>ZnIl3F5XAYqfh)wT-^Y2cn-!;v9Or&RQJowKo z81Gt}F@{vpMC&-si0v8N?3me=Sg2kPeEfw+yaIsP(+HtKkoCwm=>Thto((vb=;ZIqgy) zm}@mLDV#BX7j}^+Rz+-+VDY}Q@Ou4#Lg|92B6XJG(bEv|SRxpyjGMO31`=&fZY|Ds z3P$TL_(;bnsdfqy4Hqp81OWUGsOWVJYH-Qpf zZNU8TGh2#Q)2B?Q=(Y!(_$(}3TO1gSNR!G{_Nnw#iIzTPm8uDSiDg?rw{o|$@U5qSZSa0-V3DkoF|=@_B>@> zyBNe_*kR=LQA4`D^nc99VUV{8kJS<9U&#G0PYS}nv!1^j|9Z*=JlO+l9Xa3#=qW%1 z_8UN$1E%{MeiS=%=_~{5WMXtsPZid-~7L$AwDzNJ8paq352VnyeO!ocN>j4`3CgI)~1^YX;g6XXojZi)+01y}Y(#MU4)j3CPTR@cj+ zZ*Y6a071ObWE}5g;(XRO0LySeGJ2eK!I_T=^cEt`rDVkstr4nKP-Q}0=$5sa{jdez z;qf8@^%f4rw3D*K&FyWEq>AH2e^A$WE!GfG`b-1Vy`!?^j!l?p4RIDSPvg|$4O{HV zr)McueYp6aPu}@(Z$j#@(rt@qyO66bf;S(WJ?9NF%v{vi=;Cvne_TJGBT;9xWPQUH z`PfJi)eh&>A4?lrvuoVMUZAZ zo9E-zwgj}?xL1a=jP1!m9W7}7+Zc+bSe`hhhg(bs(6|OlKb>eT4Q<%|Qz7^hqZEa6 z{|*%JseeuV-unawFjWp%0tF6d;L!SysXL?K;$lL)outBit%C7n^$Rb>XJ?35{ z?6$VnX_fe_o_XQ2tIQyC<1gKk?6RwiJAH=c8eL%C-+8n z^XBWzW+)x@A)Zd+Df}y>gajgI^ZPx;#|fVCnfjju9k9hchui7=f#+k*1d$tG*Di#l z>Mj>+djtcxnvufHpFftpz4iAh(%>VA45gWa6X%>jX2Motc%Z2cr8dJfftCg7v48tE zg^s&V>qv|>p9kB(-Ie;NEHLgS?i_NL0t)M+T=|L+M3r|zeHS79jh?bKW3~6*6~V(= zi~-mgV{ah-V6$&AZw%97RlZ?}Ly5?Y3k8!$w~;9BaC9wY_?n>twEo>7J&)#~#irWk zz}v3wB&2OBLOJmlqoT?6b`q>e8f0y}xMf$wvWznXiX9`^LydEmGm*pcf(9{<}CwxlN(?O=rus-1m^gikKdN9aS`L9 zc$syB?qSK&MhC&NQGwNB{DP!}l-L-FP{LORzNKfQ!x|KrLG@l&-3UKMgM72a7s>UN z%s36LLmu3;j&x2SA-^&8`DLX=7US3{;ZsV}AVr~I1BLX?7$0f0M4BS0+6Bm*uV!SH zfh-@Jg3vh+rXJ5}g7HBm%=wxDhCvV_iZH_%KMzWB%)riVT&pwbzJwYrCVDEpsvm}b zEnl&`{U;B&LGKos@cg?ux9I$yNKFoJjF}=>itYt+HB(ahQ0nSwCERjjY7XY4v80>z z$TgsNgCOrYYwz#{Q-(AdTi<t7lxq?N1UF?omSSO$bKk9_C^6ta zUsBwuFpnHre;bTObZQKcjFY6?YfdAp$wI0NaUj17lk`K6xV?V+=m4wtwQMXjRPXWr zoNGN{S2jk*&|3P(yjlf%Q{fa(lX6>VX}rh!FIum8hfhn1#HbQeKIYz~gUmvd-Tk)F za@Nf|OKvNzqgdXWyQ7UkjZ~0gxcu|V3BUFH@uVzV_?Mdg#HR-&z!NXvNe}R}23X1h z*1a)cmftzl6T-BPJ|a@ZHVSuWH|GJEQaBL|)f5bxr$4G_0ln>eC*o7PYY?EmVw))2 zU_B7Wf>rl?$|mSE6p9koN$6o|tCYgSX-%qL+F?WWP6<<*9g?*czJkrLQ@xyG*oo2^ z-dj4oNd`Xj?a}>FhC_hU#s){>7^`$SQlweh{=V6{cYNKfUa8`U&`oDp_X9zGskL$S zbIJX29>=#Vea6OFDFx=sBh^zQ%;DP&n9iN*j~(U*;6*03`FCy|biBv~HjSC5qgj+{ zP-5VLRrS59JabH^pjy69Bytwsw4+Uqg2CCnaqBz<5U*kC{9`tTt0Bb=AR(n?Q|EO? z^X&4y@3ESUJWom$S*^dgtLhb{d~EGPLQL2r{YrOcMA(&(a4?uInG;NbQ1wNC_0X?;>VcC1X5A!IFnZv)p3N^a%RW{=YQ0OQ zl8uB!^!A~YoIBTJGBHZFa;=^zAf|Dbb7~^KZlg9Go$S-+MAo^r9Wnf|VVH#!fzX7~ z7*cDbmTNI$Oc_vDII6k{X<9he?2k`*J(UIYQnNkrh z{+~_Z0j6jHQ_TNYD?pF~Tme0G|CdLpYM}Gs5+!$6l90Mw z{eYoyS+UE3=;{c`o=~g$kZ`QwXF7W!zk;QL@VKHQ|2audMWpp*=As_CYAC9Da4>!R zBE{=#J;{4-6DsG}+}L>UZtowr#03y}le)LuAcX72ncpL_I}ubGq?1OCQiD896=Zjp zi5;H&dcio?oWW;bkP*Q->+?4h++BeghZT0m2fbAHMHR?v7rGPAe+w8nI#m898`yh7 z%||oDLFh4Dsf%Crq_KFh%AiO)nMkb@A%mr$(z#?)?r;H4d)80HVU82-1#z5UG?loL zL7kg~)Bb;Dr$2E@Rk-wjHdP4ZP{8wcH-3OAW55$|Mga2je@xxier#RBBn%{Lb;zgM zsE6~PN4lLZ^LY_w>S76Ln+NaM;nY%gMx5J8qy9sFT5)?E;N1&XT+|!qM@P*JvKPs@Wz!?|!aT6a-9;d$&`D%2d z%TrhY*Lo?Zv{2^<)uUN=^7YQ--Xqttc_)E7|Jlhd$Rez50K($2VV1y>@?404f?#JU zZ!XSZBDNOk8&c`UHUDXd*e`2I@g3D}Nngcv2Gvvg^Mb|K?HSfTvKv6E@%p=om%-`sxsV zKhr`Ctb5Pr8hJ>VSg>{`EXIm^t#bh+yFC_ar0Io}q|vlKG3eioH8`e0z$i0zRGrxQ zD(hNivfi(FF2g{~%nlwZkLZx2?%wTiL0#&$FSE4BGG)srGHolI?0i)GwWsEHRD=VEXXw1+#`VJXl;-?A7>L8wmY~FfM3$%1T z2X54%>b!A@{W5UW9yX&!`)79qCKSBmX+Ie=iXDESd4Wxu4d70}n{u_RepkmF(A_HA zvvqR%6wBDLYe9f(xfDSzyc$JR5YR$VmA)x5WmPA4tjQb|ajJ-_hKK&6S0LxAS2*ol z7>Ttnvhj5W<|e$aZ?Zn=2K!nh&O+7`nH`t7Z+F8un6=N<=;H7A{}Sy>8DZ4p2tw^NahTK`fS`$g?MYE? z@aHDk$=Jzjo%-^uuYcu(;?Bc8@7?6VzqX92&W!ZU%wb(!P+GU4GwfNBu|4!3M1NeV z3jnK^iKlaG4J{ zM`lBvhq!3`A25XdZS1e)_We-YlsWC5n*X_ zmI{%U$BKEWHwh7(HMLRayL@T+Q7xGMN7%l&b1^9XE_T=vwq4GoQOh!KMh(xelVn}C zYw=~O2-dr%Q|B^HyTRfie9-64PM_;CL+#yi*oNLqV8Z_)_G6!Fg%E98}p?_9fY`Y z-jHQUB-!4^@R17Awxf7|rC_m<6(}Pn!_drgjD58fsa14#lwIxl$W7C)lt>n}{R4z( z5+sB)yzwYLYNA<=Ip$Mfy{6L2YucwlKd-UiT*x5>JskKhWahEafUOrqt93@_=@9%? z`1=^z!t;m%(4>X5>JAQeT}q#Eam5_zkFR=UXQp1vUd%HjA^8)#Z5iI<50~>s=@{@@ zFo#hEcSwb$J2!6&D&%4O(oQlAUxr8bbKn~sem1nS(v zgMAxjF*>d50*g_n-MJ@sCm9lCMIgOZ)-`Hzfaa^4g&guPDDSmX%{_XbHi>DuCDrjc z4dwPZ+U{2%Y3P!z9cH;y5T4*ai2k_J5dOUz%HNHD>j*i3E1*7g4fJoq(}92>aR9Dj zfhm5wdU3brIz~gLV(IiNP213oP#t0mPp4H=S&0u!r_lo4HdPYfvGD_m@Wi2c;Met> z$R;(GdAHuh=Gty$sDAi`VK4uTqcPHEqHU4fdF7Y54bQLL&q3Canr4`&9QVg~Ubd!M zqIGqe*Jw*ShZUx>CVjUh&X5xhP`E(~)+iiHc$f%y#kW?PqNcA@J0>M>^A**3K1~(W z52woS74NwuAM&br3J>Lmc+HuZtQW7t0=_jkd%{CmcmI3%l5l@4u{NKq6zI``?!A{Vk*C)%MJh|M#?9B+pjP7cuz zoXSjSaA9I?S+ytT#nm;T7)>%v3&|`;X9jN<{7afEw;&9EUQEljKv7l>K2Zyhi&H;+ zU@QScon?z})?pGaWs@y!lR3jForvS=*`!QWg3h7#JKfSZ6?XH#kQYw5<8||m=Pi#= z!7`A(Pncq{l`pK-+Iu4SuEsuKxiALo)T!ytViNwkVg*aW_i>Q}T_3dnwFY;Tsj;VI z2Vob~4jkc!Lu$?{VN4Cd7kMvTjAetL;k=%aPWJ2Suw#C1#Y<6baInZwpQ0!I>^nD> zG4-~#1HvoNlS4hM`ZNgc5A^FX$UU;jq*$Caiq zuqyuL_}A6%=M}U8SKNRrAk#+zd?D-sbNsF+{IJCl8Y6%ykc4pf)><-nrs|Kw=ogWn zR%C%r&ysS3Rt^2h4+K0%Cjm2Mmbe~z(47-~@J-qVR9@WiOZzkoehL(8op ztdaL*7;Ci3eHbAHW25AiG9^*|s)lx2__32#XV7Ojz~t#s5z}jEF>BI)a&C^W#t_a+ z`}P6ClyyR=2YSk52gEpK&;KPp5+r=YGBYtYfpq_+3$G`MmRPf0QWtdCL%j;oO+#ns zOLqoH@NyCc7kOd$9q=2+O%C;63#Ak=N70c=!HbyZ-n@t%`3cZb#JD_KT$ua{rOE7f zWD3fMa@1`VJ@ZFVRmELx-%oy)UXLo76mK4dzx%3ibd|jFT|xLJ{dQU+J@?u-RJ0jU z-bJf@dp3Q&+g_08D)Y4xly#?*jsY909mG1!nnO zS9mUV%0?uC(NW5j;?7esD8&kyc1Vx~udS?J7IVWgL8@lJo$-5z!0%oiOk1{?9XE?& z$%y@d{>zRHL0Da2hUSLK>B8hieqgrFxDb{G_~DPT6%(ByCeOAV_jiW^`SlF-B4GYO zBl_7FJhQ1o$(cf%bX!j}%7`Jmp|o2AQG9Q9Og|%dCNja6eXKN!KK~kI~|KNK?nQD^fjZ%2j4TZ*3dmm`;w;+mdT-o2=ftm zYVy15s7a%5o00axIp#L~8hRCexRbNBqB^k{yR1TK=0OJi993tVkQF=p8M#gE7&$~4 zjy*^#+2^HI5vC)@d0|1c+{QbBIDBmiKPKw}Xsa5;{sBVw^Y|0q@gTj+XE#R*Ejr{? zy^34@nPOe(c5@Ba;gHrot$~}oPZTlHQL{hptl=B^E~7 zyT|yva zN{mP!`X!|j1VTeDF5dW_Q7b6p2s}1(eDo~!Kp6pM`tz3amCW|&#+N~V@NJTa;Cpjr zx2I(}r7ReV-gW^T=(yxxkb^C4#cx<2P~t6*dT1ZNM>t1HmqKzkmxH$v8>}FmjV4#w z9#mO->_4FCKxG9nfTa+t$f;H*1rr6+s`U%CXie6|P+t)2v1tO^LV9&Ns|(~NsNwie zY0;}umCrQHF!K4JV^ze~d>C~f*c~lHsk9pZxeo(MF25r2-S<{J2KHt0VXQlVhk`_V zNeADc@`ytx%~}d>J?MVj$OUx*eXGdRsSll4iu`LtYqJ=8be109M@a|SQ%$O?z6qV= z{Kt+bpn7W)_`- z9$-mJsORM0f7`#@Ki#L&M)8FWVLJDwl2bi|x{Qd_(5Uh@BmdZv(kq zl>W-0bQ-UPBk1yr_ZHI%nA(BUhq}s&mpI$fJ4F3V|d zmR?}3pQhQG?L@|3-=+QkwseG>|5sbL%YduJU%#$^br8M=h*^NI-_ObXW9zk9T9E?> zImUXkO*8R}dU4v2ivpswK6$8JyxHtW9by#r2vGMY5C}>AdsMQZ)r3S@Uma2VCg!i_TvGgL(8@~^={(MhPkmY7J*KK zTDY%C;sYm|98*(`{_?&JM&cWv&WQRD3VA}u)os%qk>i^tovm{jj66~xv_N6VlH2FJ zV=d}h-ZH|AZ)Tcce$8lI{4Kb>(YtBKO$M+tv$RwCqX`+4P87AVjD?P?kVY$DSm^kh z*<9seUnAIQ)}798=GS~j)TKFNjb+arV zE?T+rS+uLM%!~;`R--~wlXK&3t?Z4^TeOE~p2$wc*kLG_Aq%*Swsl2Sb3u$G-bF@T zijyvD+uM&QWv-Dx^b6x7*DQ%}uedGpiHAiheT+z!yvYORtW+-QC3fAa4FYfjUX!)0 zbx;3qj}8@xKlyT`q~O_?#hzN6#ORTN=mJO8y7>I%31&#tThSHQb4 zuYTVHlm=oJ(36gU_8(V28FJD$u;9FUH$W|G_~cCR#LkvN6P2l6i;m1@Z4Trbn!_>A z`5_?G-YI`3bD2L&E7(jUf_#~Ij{3Y?Kj&kS6WqZQt^PpGvT}R?6_k7B%glrm#ubQ~s4UkZEaCVWRU}^auGj&W<=lYS z{5FVdY31&edj8DM6%j;=IB$Crt>-|nC)I5Zht9i|I20e-qX9@46URg$Rm4Wt>%=CE z*vGHV(b+?kjHVD1weGhZAfDp(IoX0Hhc1X`@jmZ{n3+aXzN;y=38)lDBLoz&Ij|#} zVe5r_$xhI8$se`rUE_-AJs#<6|>%@iR{z}D| zpq%BkX(B8Q*z(ECy9_jO4@*-f@B;GaSl(s14S3nw=(v_Po3fED6VyqmU0%XgPv}$y z_c386a&UQ|-m+ZA%I<7+P7arDTF2OPuo@3qp1eQ3s>frb^ONUu+R?2is{upLenID6JWX45xSdf%>&B7f}g^0V{-19nHL7J`AemfK&N$NM6A%~zf^>J|E7 zGA&#_1o{4gn&G+g?J$Z%J)LDnTjbpOUpC-W{6B@4w5EEqjy}>U$_~^NSv3^JmEa*w zjdjkd+`Jh6L+Q_Ys3+Y1KilF1Y%Kw{z5=FT09$}g1TSEgf9efIr9{JcaL5N`b(Yyf zdCqvE2~o|sTKa;^)BRL<-dC8jO0qveK=kj-B&x$_b3b4Qig`aR&lu-s3LEnjzkG+Y zWuY_Bd|@&i>dITg72JHpwT?D^Oo0TYWe~;)7H)xL!ITnn51Ap0Cmx8P#2w8|*GRMH zs(U94FAJT)1hGgiro!o6)xvxyYLYwmDEUmMEvnZ+?l7dpi=K{aOQHRZgI8GeHP2VW z1w(r^MDG{!t?gQ!PhB}rFYH$5Qv8o!2(YlNp98db5Mo2^4ULH?5fSo^?Tl5Q@h?3P zs_3I9Yx?C`vu0?<5TvWp{Io0bF-#vj3{X`zi@vA{?FNS{Sa}_z*Q`}=#13nKcXeRU zs%~;YGn;8|4yg}BuISvXoDD5o{>Vm<{VLrunJiz4`Sy$D>kl#ri(pvdy<{`{p*8}- zix;V19DHC`cEHwmxyq=+7xHuJXLDr%4UA-z2bxruaa6oa;K>y3ZKJCtgd)geiO-|G zOqz?Fx&^$;{nN_ozDR4_yE3?TTAwvww%pF)yD4}9_ykAVY zSHC-wF-or6usXB4aT@NTG$?JUYkAv!?oP3wm-J z=*n+_Ie@WeRj(IRi*3ZpS`>^^$=rWsQ#0GTi$%>%oKwi;f>fizSl*{jO;?|u0cC7Q zyrXM?cOoRB2ywc1*LiJ1gVr!f`L4*|;jsZJ0wR_0wApwp1f zZh{dKvm~Cw=@1gyM_JGe=7+ywCUpUE8)0#D#H1slK&ul};n?g8D zP~w0!VggxjI+4^#7Jd99+63k9nG_xGRJODxLeIN)^nsR*#c1Vy=R@;d%WKIo;}$vR z60!DQ2l9Iql>>0wE5&K?ymU}+h(5zf*7%&?1>u6U6p9@~MXgnQD8!4r!Bi~dTtd`1 zzRf_zJxDlKl0qIxgrATjNxE}+FXXq@jj|`RrxUC;Bwe(hMi&q5uUAV{M?~(dzQIA` zU)uDwn;$-KU&i+K*B38d$1<=f;JTC_`>lYxhkDU+iRxLPM+elWff+hiUx$-?$W*EE zWB8RC6T38P4{ae@CrdXjqsT0>W>ucJYb5|PV)wLwqQE;%6qar>lw>OZ{pZ?Bp+C&4 z)^+o{K5(@$FW8x>k$EYK>FvUp$NlhhH=;yqI`)X!oJ1|e@AL~LS6T89YPNwc^=LP5 zKYTq}Ip!D%mY`^yL{wZ`WBWN+$Md-Qs`cg9BE{4z1`oe;enz)jx@=N#OTit_+2H|5DYTdc%z)V9F9WssLAcfU9l5 zl^rn4KlRa;R|Pi`_{kV_9!oz%juU*bNSb08$Lqv#-95&^Ty=>$>Q5kkewe~1(2e1i z#f*`OB20dFEcsaI_TM?Q_^uES(q8qm+KrQmUE?h*o{RA8K6A&duAmxIQl~rp@?4Ex zh27WQJ2&+I;Iz^7I(?GtdGv^-@9IhIatfz#Y{YR8$okPUUU-<3I#DKF`zzjl`D#oj zIC4YkTZC%=A-35CoLe{g+5+AB$Ht15+ue30Vf*)k>;s@CFzSOMPIHaEq!<2hn~_d9 zlc9+{BA?Ur*Uc5>vDJpNKW*$u`F~(wWk4Y^y3nzJElLsPxt^u}nOF3m#q1a@pz z!vuvrRWN&VsO;w)FGQ=_qTadzYUF$67ZZ3#ngkvlis`aeXELvJP%ggk&w8abdOCvw52gR z3HICwfhvK~AcZ87+*h+4u4=YP%@d?Vtn^frxe#Re*0#sg?U1AJX!^t&@$5?c+tHu= z%1F5X_uqdt{*700fT>jA_c<8UCqbXJfD!9q7&1eid=Z_F3osMyF6TGPSafJfc+_Dh7FTlc4b}?Is3W2&s zU@Gop>^r{Z6%p6P9-^ogwFO?}B1hXf?Z)t=5%DP!bu=Z_pmxfu0FPwndsH${KFi1ji5{US=x$K%U^(s^@ z4%)yi&gow3Dlpt=d7`ThqAM(NCd4|n4rPNXNw^WVV118}u5=007+wUg_;3XMa6y+~ zwn5PpJ+k^%iKrf#JOOSnVb(?l1BTjd)$AK}gWiHQ8iBdB6j!9AIvCXy%(twuLwtsZ z)w({tj$87@@8h(?U97y`R#JCik^8^*Z5a+$WlAG z6e>c+&QLO1u-G^*DBE7FtDta$3{ipo;d0|2LVsKt3lIL4sox*^*Ht-SDjdkKfGWlB^$%GXINt#O z{~fP>K+dfjQ6axSjH_!+3#{SW912kz2@H(1+jKi3eU#e7W^r{5^Ydf1>s4>Xt-xr* zi};?6(;q0pc9soYME3>rYiJV%DLQrwom5cOUGeprH<+A85tfev2YUZfdd=bEy_@L9 zLO+=np7OZ_Zv{(`XzWm@z*l*sa&){b<27uObh5;Z9T6S5Y_7r|=Jv%Cj^ndK$DsSm zwVx8qt{cPMTqQ~$gD_lv=1;*!sQBS>KlLv;TV6VB)_H~UJlRlxecoMZiK?e&y6-g; zc=&&-(L8kgO0Mfe^JhgjXPRG($8$s5N` z?jP?K*WG<1;(zjwMp!}4F5R--^z2jGy4-oM%okq?IQqE3=gZ%g-cWIm zq{Gq098sA;T*p2O{Jrz=cTeXM$5ppZZRXr75iU{uEkyUiC7x_0jr4QX3JaPu8dj~6 z`7?i8_uRda8s*nnjtc=dBFEIvQ%{l-<|JbHpu46 zdojU9aGBh<1Fu0z1fH&vosNO{L;+a3Dg~x1V1Muf8!$cTF)#=;1L?Ow6F}+ekIMV| z`;|Gj_gr$k(a^5*SI#x4aMhkE7KsAKUcOA7(JrOPQYaw6crdMkUC^d#rMpb$yx1Rn zQzB}ayV{rfzngks#;%)-=VUAkf6jWbXQhG0igVYRC!Dt3bo~ayt%z;Ej&!HrE}9U| zlzygS_Y*y@Z4tcQ9oa=+;+X%&D03~3ow8x0G2=yD!DENA&DgzE=E^Zmf3-Cxe)mI; zSpkkkK8k@ZuQpG%ox(XmOYrT=Z_`#f2kqucd=s$rPHAIE^99BSnH8t^+f}W8{GK)O zO&OE!dsl1!irhv)Qyy*U&^7<6eJZYp*Q(9hziw9U%8YBL7AtLes2kd{Yo%StgG|XU zPpVfL%vheaJnrEnm4}J5Z*rVw^}i}0HL0aiH2YQAn%hZzf_oZzulY4BKNm9DdxCi$ zt4HPWa=rIK3pQ`xm%Zro<*O^2J)3XL?OA&y(Q&EweB&UFZFzsCkFrk*asIMY{_%$` zd-h(JxIahxGDAhA#7w<^FP1&w*5xsutt5IwDPz})L$fx`pVa$mVdzEc#l;Uyg+4mJ zHPsaK5S$Y8U+R4JzG5}CqjwyOCiVPu5KWcj*s$^z(}UBQx+^VNBd^>z{JKUvw&_Me zy8e}eE>4-VN^8=6p6M*E6SHLE-}BC)#DGUcS!!2yB7b>?Tk*5G*)jXBdfrOgZ|G$C Q=)8KpLzi0FR?Dyl08FkjIRF3v literal 0 HcmV?d00001 diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index df4432a62bf..ec40b57e00a 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -116,7 +116,9 @@ "COORDINATES_TOPOLOGY", "NUCLsel", "GRO_empty_atom", "GRO_missing_atomname", # for testing GROParser exception raise - "ENT" #for testing ENT file extension + "ENT", #for testing ENT file extension + "RANDOM_WALK", + "RANDOM_WALK_TOPO" # garbage topology to go along with XTC positions above ] from pkg_resources import resource_filename @@ -329,6 +331,8 @@ # Contains one of each residue in 'nucleic' selections NUCLsel = resource_filename(__name__, 'data/nucl_res.pdb') +RANDOM_WALK = resource_filename(__name__, 'data/xyz_random_walk.xtc') +RANDOM_WALK_TOPO = resource_filename(__name__, 'data/RANDOM_WALK_TOPO.pdb') # This should be the last line: clean up namespace del resource_filename From bc8d15a72cb7eb4004bfea2b901ae6151b358401 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 10 Aug 2016 11:56:44 -0700 Subject: [PATCH 31/39] Fixed tests for new random walk files. --- testsuite/MDAnalysisTests/analysis/test_pca.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 3da1bc8bfd7..46279c01a2e 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -24,7 +24,8 @@ assert_array_almost_equal, raises) -from MDAnalysisTests.datafiles import PDB, XTC, RANDOM_WALK, waterPSF, waterDCD +from MDAnalysisTests.datafiles import (PDB, XTC, RANDOM_WALK, RANDOM_WALK_TOPO, + waterPSF, waterDCD) class TestPCA(object): def setUp(self): @@ -76,7 +77,7 @@ def test_transform_universe(self): pca_test.transform(u2) def test_cosine_content(self): - rand = MDAnalysis.Universe(RANDOM_WALK) + rand = MDAnalysis.Universe(RANDOM_WALK_TOPO, RANDOM_WALK) pca_random = pca.PCA(rand.atoms).run() dot = pca_random.transform(rand.atoms) content = cosine_content(dot, 0) From 6257bd762194cf325305dee6952c252d8e5439b1 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 10 Aug 2016 12:36:26 -0700 Subject: [PATCH 32/39] Format change for mean progress meter. --- package/MDAnalysis/analysis/pca.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index f4e8e9257f0..1a965254eb3 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -197,8 +197,11 @@ def _prepare(self): if self._calc_mean: interval = int(self.n_frames // 100) interval = interval if interval > 0 else 1 + format = ("Mean Calculation Step" + "%(step)5d/%(numsteps)d [%(percentage)5.1f%%]\r") mean_pm = ProgressMeter(self.n_frames if self.n_frames else 1, - interval=interval, quiet=self._quiet) + interval=interval, quiet=self._quiet, + format= format) for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): if self.align: mobile_cog = self._atoms.center_of_geometry() From 159b9180b72df37a9f98d4f035ae20a7f2c9563d Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 10 Aug 2016 13:17:30 -0700 Subject: [PATCH 33/39] Moved cosine content, changed and added documentation. --- package/MDAnalysis/analysis/pca.py | 24 +++++++++++++++++++ .../MDAnalysisTests/analysis/test_pca.py | 12 ++-------- 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 1a965254eb3..9c4d78bb8fe 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -82,6 +82,7 @@ import numpy as np +from scipy.integrate import simps from MDAnalysis.core import AtomGroup from MDAnalysis import Universe from MDAnalysis.analysis.align import _fit_to @@ -299,3 +300,26 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, dot[i] = np.dot(xyz, self.p_components[:, :n_components]) return dot + + def cosine_content(self, pca_space, i): + """Measure the cosine content of the PCA projection. + + Cosine content is used as a measure of convergence for a protein + simulation. If this function is used in a publication, please cite: + + .. [BerkHess1] + Berk Hess. Convergence of sampling in protein simulations. Phys. Rev. E + 65, 031910 (2002). + + Parameters + ---------- + pca_space: array, shape (number of frames, number of components) + The PCA space to be analyzed. + i: int + The index of the pca_space to be analyzed for cosine content + """ + t = np.arange(len(pca_space)) + T = len(pca_space) + cos = np.cos(np.pi * t * (i + 1) / T) + return ((2.0 / T) * (simps(cos*pca_space[:, i])) ** 2 / + simps(pca_space[:, i] ** 2)) diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 46279c01a2e..1fd1dd49946 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -18,15 +18,14 @@ import MDAnalysis import MDAnalysis.analysis.pca as pca from MDAnalysis.analysis.align import _fit_to -from scipy.integrate import simps from numpy.testing import (assert_almost_equal, assert_equal, assert_array_almost_equal, raises) - from MDAnalysisTests.datafiles import (PDB, XTC, RANDOM_WALK, RANDOM_WALK_TOPO, waterPSF, waterDCD) + class TestPCA(object): def setUp(self): self.u = MDAnalysis.Universe(PDB, XTC) @@ -80,12 +79,5 @@ def test_cosine_content(self): rand = MDAnalysis.Universe(RANDOM_WALK_TOPO, RANDOM_WALK) pca_random = pca.PCA(rand.atoms).run() dot = pca_random.transform(rand.atoms) - content = cosine_content(dot, 0) + content = pca_random.cosine_content(dot, 0) assert_almost_equal(content, .99, 1) - - -def cosine_content(pc, i): - t = np.arange(len(pc)) - T = len(pc) - cos = np.cos(np.pi * t * (i + 1) / T) - return (2.0 / T) * (simps(cos*pc[:, i])) ** 2 / simps(pc[:, i] ** 2) From fc2f4d45a3c0696892fdd84180799e0519f6461d Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 10 Aug 2016 13:27:23 -0700 Subject: [PATCH 34/39] Fix up some docs --- package/MDAnalysis/analysis/pca.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 9c4d78bb8fe..0880c6fa150 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -305,11 +305,8 @@ def cosine_content(self, pca_space, i): """Measure the cosine content of the PCA projection. Cosine content is used as a measure of convergence for a protein - simulation. If this function is used in a publication, please cite: - - .. [BerkHess1] - Berk Hess. Convergence of sampling in protein simulations. Phys. Rev. E - 65, 031910 (2002). + simulation. If this function is used in a publication, please cite + [BerkHess1]_. Parameters ---------- @@ -317,6 +314,18 @@ def cosine_content(self, pca_space, i): The PCA space to be analyzed. i: int The index of the pca_space to be analyzed for cosine content + + Returns + ------- + A float reflecting the cosine content of the ith projection in the PCA + space. The output is bounded by 0 and 1, with 1 reflecting an agreement + with cosine while 0 reflects complete disagreement. + + References + ---------- + .. [BerkHess1] + Berk Hess. Convergence of sampling in protein simulations. Phys. Rev. E + 65, 031910 (2002). """ t = np.arange(len(pca_space)) T = len(pca_space) From 50e329fd4497a0f1ae15d2444fb99511e5e3d550 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 11 Aug 2016 13:54:38 -0700 Subject: [PATCH 35/39] Added autofunction autoclass docs, made cosine_content static method, updated tests accordingly. --- package/MDAnalysis/analysis/pca.py | 60 ++++++++++--------- .../MDAnalysisTests/analysis/test_pca.py | 2 +- 2 files changed, 34 insertions(+), 28 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 0880c6fa150..18b535fce09 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -75,6 +75,11 @@ From here, inspection of the pca_space and conclusions to be drawn from the data are left to the user. + +Functions +--------- +.. autoclass:: PCA +.. autofunction:: cosine_content """ from six.moves import range import logging @@ -301,34 +306,35 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, return dot - def cosine_content(self, pca_space, i): - """Measure the cosine content of the PCA projection. - Cosine content is used as a measure of convergence for a protein - simulation. If this function is used in a publication, please cite - [BerkHess1]_. +def cosine_content(pca_space, i): + """Measure the cosine content of the PCA projection. - Parameters - ---------- - pca_space: array, shape (number of frames, number of components) - The PCA space to be analyzed. - i: int - The index of the pca_space to be analyzed for cosine content + Cosine content is used as a measure of convergence for a protein + simulation. If this function is used in a publication, please cite + [BerkHess1]_. - Returns - ------- - A float reflecting the cosine content of the ith projection in the PCA - space. The output is bounded by 0 and 1, with 1 reflecting an agreement - with cosine while 0 reflects complete disagreement. + Parameters + ---------- + pca_space: array, shape (number of frames, number of components) + The PCA space to be analyzed. + i: int + The index of the pca_space to be analyzed for cosine content - References - ---------- - .. [BerkHess1] - Berk Hess. Convergence of sampling in protein simulations. Phys. Rev. E - 65, 031910 (2002). - """ - t = np.arange(len(pca_space)) - T = len(pca_space) - cos = np.cos(np.pi * t * (i + 1) / T) - return ((2.0 / T) * (simps(cos*pca_space[:, i])) ** 2 / - simps(pca_space[:, i] ** 2)) + Returns + ------- + A float reflecting the cosine content of the ith projection in the PCA + space. The output is bounded by 0 and 1, with 1 reflecting an agreement + with cosine while 0 reflects complete disagreement. + + References + ---------- + .. [BerkHess1] + Berk Hess. Convergence of sampling in protein simulations. Phys. Rev. E + 65, 031910 (2002). + """ + t = np.arange(len(pca_space)) + T = len(pca_space) + cos = np.cos(np.pi * t * (i + 1) / T) + return ((2.0 / T) * (simps(cos*pca_space[:, i])) ** 2 / + simps(pca_space[:, i] ** 2)) diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 1fd1dd49946..406b875f22a 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -79,5 +79,5 @@ def test_cosine_content(self): rand = MDAnalysis.Universe(RANDOM_WALK_TOPO, RANDOM_WALK) pca_random = pca.PCA(rand.atoms).run() dot = pca_random.transform(rand.atoms) - content = pca_random.cosine_content(dot, 0) + content = pca.cosine_content(dot, 0) assert_almost_equal(content, .99, 1) From a52d9b86990109f9c8a12430c0b1bc6c5c8f5159 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 12 Aug 2016 11:10:15 -0700 Subject: [PATCH 36/39] Doc fixes, test updates --- package/MDAnalysis/analysis/pca.py | 74 ++++++++++--------- .../MDAnalysisTests/analysis/test_pca.py | 11 ++- 2 files changed, 49 insertions(+), 36 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 18b535fce09..8be23708b40 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -22,20 +22,26 @@ :Year: 2016 :Copyright: GNU Public License v3 -This module contains the linear dimension reduction method Principal -Component Analysis. This module constructs a covariance matrix wherein each -element of the matrix is denoted by (i,j) row-column coordinates. The (i,j) -coordinate reflects the influence of the of the ith frame's coordinates on the -jth frame's coordinates of a given trajectory. The eponymous components are the +This module contains the linear dimensions reduction method Principal Component +Analysis (PCA). PCA sorts a simulation into 3N directions of descending +variance, with N being the number of atoms. These directions are called +the principal components. The dimensions to be analyzed are reduced by only +looking at a few projections of the first principal components. To run a +Principal Component Analysis, please refer to the :ref:`PCA-tutorial`. + +This module constructs a covariance matrix wherein each element of the matrix +is denoted by (i,j) row-column coordinates. The (i,j) coordinate is the +influence of the of the ith frame's coordinates on the jth frame's coordinates +of a given trajectory. The eponymous components are the eigenvectors of this matrix. -For each eigenvector, its eigenvalue reflects the variance that the eigenvector -explains. This value is made into a ratio. Stored in -:attribute:`cumulat_variance`, this ratio divides the accumulated variance -of the nth eigenvector and the n-1 eigenvectors preceding the eigenvector by -the total variance in the data. For most data, :attribute:`explained_variance` -will be approximately equal to one for some n that is significantly smaller -than the total number of components, these are the components of interest given +For each eigenvector, its eigenvalue is the variance that the eigenvector +explains. Stored in :attribute:`cumulated_variance`, a ratio for each number of +eigenvectors up to index i is provided to quickly find out how many principal +components are needed to explain the amount of variance reflected by those i +eigenvectors. For most data, :attribute:`cumulated_variance` will be +approximately equal to one for some n that is significantly smaller than the +total number of components, these are the components of interest given by Principal Component Analysis. From here, we can project a trajectory onto these principal components and @@ -68,8 +74,7 @@ `cumulated_variance`. The value at the ith index of `cumulated_variance` is the sum of the variances from 0 to i. - >>> n_pcs = next(x[0] for x in enumerate(PSF_pca.cumulated_variance) - >>> if x[1] > 0.95) + >>> n_pcs = np.where(PSF_pca.cumilated_var > 0.95)[0][0] >>> atomgroup = u.select_atoms('backbone') >>> pca_space = PSF_pca.transform(atomgroup, n_components=n_pcs) @@ -97,7 +102,12 @@ class PCA(AnalysisBase): - """Principal component analysis on an MD trajectory + """Principal component analysis on an MD trajectory. + + After initializing and calling method with a universe or an atom group, + principal components ordering the atom coordinate data by decreasing + variance will be available for analysis. Please refer to the + :ref:`PCA-tutorial` for more detailed instructions. Attributes ---------- @@ -118,14 +128,7 @@ class PCA(AnalysisBase): mean_atoms: MDAnalyis atomgroup After running :method:`PCA.run()`, the mean position of all the atoms used for the creation of the covariance matrix will exist here. - start: int - The index of the first frame to be used for the creation of the - covariance matrix. - stop: int - The index to stop before in the creation of the covariance matrix. - step: int - The amount of frames stepped between in the creation of the covariance - matrix. + Methods ------- transform(atomgroup, n_components=None) @@ -183,10 +186,12 @@ def __init__(self, atomgroup, select='all', align=False, mean=None, self.n_components = n_components self._n_atoms = self._atoms.n_atoms self._calculated = False + if mean is None: - logging.warn('In order to demean to generate the covariance matrix\n' - 'the frames have to be iterated over twice. To avoid\n' - 'this slowdown, provide an atomgroup for demeaning.') + warnings.warn('In order to demean to generate the covariance ' + 'matrix the frames have to be iterated over twice. ' + 'To avoid this slowdown, provide an atomgroup for ' + 'demeaning.') self.mean = np.zeros(self._n_atoms*3) self._calc_mean = True else: @@ -207,8 +212,9 @@ def _prepare(self): "%(step)5d/%(numsteps)d [%(percentage)5.1f%%]\r") mean_pm = ProgressMeter(self.n_frames if self.n_frames else 1, interval=interval, quiet=self._quiet, - format= format) - for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]): + format=format) + for i, ts in enumerate(self._u.trajectory[self.start:self.stop: + self.step]): if self.align: mobile_cog = self._atoms.center_of_geometry() mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, @@ -216,7 +222,6 @@ def _prepare(self): self._atoms, mobile_com=mobile_cog, ref_com=self._ref_cog) - self.mean += mobile_atoms.positions.ravel() else: self.mean += self._atoms.positions.ravel() mean_pm.echo(i) @@ -273,6 +278,7 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, step: int, optional Number of frames to skip over for PCA transform. Default: None becomes 1. + Returns ------- pca_space : array, shape (number of frames, number of components) @@ -310,16 +316,18 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None, def cosine_content(pca_space, i): """Measure the cosine content of the PCA projection. - Cosine content is used as a measure of convergence for a protein - simulation. If this function is used in a publication, please cite - [BerkHess1]_. + The cosine content of pca projections can be used as an indicator if a + simulation is converged. Values close to 1 are an indicator that the + simulation isn't converged. For values below 0.7 no statement can be made. + If you use this function please cite [BerkHess1]_. + Parameters ---------- pca_space: array, shape (number of frames, number of components) The PCA space to be analyzed. i: int - The index of the pca_space to be analyzed for cosine content + The index of the pca_component projectection to be analyzed. Returns ------- diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index 406b875f22a..ed430f3c91a 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -36,7 +36,7 @@ def setUp(self): def test_cov(self): atoms = self.u.select_atoms('backbone and name CA') - xyz = np.zeros((10, 642)) + xyz = np.zeros((self.pca.n_frames, self.pca._n_atoms*3)) for i, ts in enumerate(self.u.trajectory): xyz[i] = atoms.positions.ravel() @@ -45,6 +45,9 @@ def test_cov(self): def test_cum_var(self): assert_almost_equal(self.pca.cumulated_variance[-1], 1) + l = self.pca.cumulated_variance + l = np.sort(l) + assert_almost_equal(self.pca.cumulated_variance, l, 5) def test_pcs(self): assert_equal(self.pca.p_components.shape, @@ -69,13 +72,15 @@ def test_transform_mismatch(self): assert_equal(self.pca_space.shape, (self.u.trajectory.n_frames, 1)) - def test_transform_universe(self): + @staticmethod + def test_transform_universe(): u1 = MDAnalysis.Universe(waterPSF, waterDCD) u2 = MDAnalysis.Universe(waterPSF, waterDCD) pca_test = pca.PCA(u1).run() pca_test.transform(u2) - def test_cosine_content(self): + @staticmethod + def test_cosine_content(): rand = MDAnalysis.Universe(RANDOM_WALK_TOPO, RANDOM_WALK) pca_random = pca.PCA(rand.atoms).run() dot = pca_random.transform(rand.atoms) From e0617e07dbd5048be30a9033f99f73bcd0bec9fd Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 12 Aug 2016 13:43:06 -0700 Subject: [PATCH 37/39] More docs fixes --- package/MDAnalysis/analysis/pca.py | 11 ++++++++--- .../source/documentation_pages/analysis/pca.rst | 1 + .../source/documentation_pages/analysis_modules.rst | 1 + 3 files changed, 10 insertions(+), 3 deletions(-) create mode 100644 package/doc/sphinx/source/documentation_pages/analysis/pca.rst diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 8be23708b40..8bb49c24b62 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -31,7 +31,7 @@ This module constructs a covariance matrix wherein each element of the matrix is denoted by (i,j) row-column coordinates. The (i,j) coordinate is the -influence of the of the ith frame's coordinates on the jth frame's coordinates +influence of the of the ith coordinates on the jth coordinates of a given trajectory. The eponymous components are the eigenvectors of this matrix. @@ -106,8 +106,13 @@ class PCA(AnalysisBase): After initializing and calling method with a universe or an atom group, principal components ordering the atom coordinate data by decreasing - variance will be available for analysis. Please refer to the - :ref:`PCA-tutorial` for more detailed instructions. + variance will be available for analysis. As an example: + >>> pca = PCA(atomgroup, select='backbone').run() + >>> pca_space = pca.transform(atomgroup.select_atoms('backbone'), 3) + geneates the principal components of the backbone of the atomgroup and + then transforms those atomgroup coordinates by the direction of those + variances. Please refer to the :ref:`PCA-tutorial` for more detailed + instructions. Attributes ---------- diff --git a/package/doc/sphinx/source/documentation_pages/analysis/pca.rst b/package/doc/sphinx/source/documentation_pages/analysis/pca.rst new file mode 100644 index 00000000000..dc90112ad4b --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/analysis/pca.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.analysis.pca diff --git a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst index d4e241cd6a9..6c3e41a581c 100644 --- a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst @@ -113,3 +113,4 @@ Dimension Reduction :maxdepth: 1 analysis/diffusionmap + analysis/pca From 6099a79679cb3e9afd64c12cf3a529dc2823571b Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 12 Aug 2016 14:01:30 -0700 Subject: [PATCH 38/39] Doc fix [skip ci] --- package/MDAnalysis/analysis/pca.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 8bb49c24b62..9fb08210f39 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -74,7 +74,7 @@ `cumulated_variance`. The value at the ith index of `cumulated_variance` is the sum of the variances from 0 to i. - >>> n_pcs = np.where(PSF_pca.cumilated_var > 0.95)[0][0] + >>> n_pcs = np.where(PSF_pca.cumulated_var > 0.95)[0][0] >>> atomgroup = u.select_atoms('backbone') >>> pca_space = PSF_pca.transform(atomgroup, n_components=n_pcs) @@ -109,7 +109,7 @@ class PCA(AnalysisBase): variance will be available for analysis. As an example: >>> pca = PCA(atomgroup, select='backbone').run() >>> pca_space = pca.transform(atomgroup.select_atoms('backbone'), 3) - geneates the principal components of the backbone of the atomgroup and + generates the principal components of the backbone of the atomgroup and then transforms those atomgroup coordinates by the direction of those variances. Please refer to the :ref:`PCA-tutorial` for more detailed instructions. From 1cf340c71b376f7c204ca340d5fdf3c7514f6221 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sat, 13 Aug 2016 11:31:44 -0700 Subject: [PATCH 39/39] PCA docs updates --- package/MDAnalysis/analysis/pca.py | 20 +++++++++---------- .../MDAnalysisTests/analysis/test_pca.py | 1 + 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/analysis/pca.py b/package/MDAnalysis/analysis/pca.py index 9fb08210f39..2acce16c4b7 100644 --- a/package/MDAnalysis/analysis/pca.py +++ b/package/MDAnalysis/analysis/pca.py @@ -26,14 +26,13 @@ Analysis (PCA). PCA sorts a simulation into 3N directions of descending variance, with N being the number of atoms. These directions are called the principal components. The dimensions to be analyzed are reduced by only -looking at a few projections of the first principal components. To run a -Principal Component Analysis, please refer to the :ref:`PCA-tutorial`. +looking at a few projections of the first principal components. To learn how to +run a Principal Component Analysis, please refer to the :ref:`PCA-tutorial`. -This module constructs a covariance matrix wherein each element of the matrix -is denoted by (i,j) row-column coordinates. The (i,j) coordinate is the -influence of the of the ith coordinates on the jth coordinates -of a given trajectory. The eponymous components are the -eigenvectors of this matrix. +The PCA problem is solved by solving the eigenvalue problem of the covariance +matrix, a 3N x 3N matrix where the element (i, j) is the covariance between +coordinates i and j. The principal components are the eigenvectors of this +matrix. For each eigenvector, its eigenvalue is the variance that the eigenvector explains. Stored in :attribute:`cumulated_variance`, a ratio for each number of @@ -56,13 +55,14 @@ :data:`~MDAnalysis.tests.datafiles.DCD`). This tutorial shows how to use the PCA class. -First load all modules and test data :: +First load all modules and test data >>> import MDAnalysis as mda >>> import MDAnalysis.analysis.pca as pca >>> from MDAnalysis.tests.datafiles import PSF, DCD -Given a universe containing trajectory data we can perform PCA using -:class:`PCA`:: and retrieve the principal components. +Given a universe containing trajectory data we can perform Principal Component +Analyis by using the class :class:`PCA` and retrieving the principal +components. >>> u = mda.Universe(PSF,DCD) >>> PSF_pca = pca.PCA(u, select='backbone') >>> PSF_pca.run() diff --git a/testsuite/MDAnalysisTests/analysis/test_pca.py b/testsuite/MDAnalysisTests/analysis/test_pca.py index ed430f3c91a..29044fca5ac 100644 --- a/testsuite/MDAnalysisTests/analysis/test_pca.py +++ b/testsuite/MDAnalysisTests/analysis/test_pca.py @@ -27,6 +27,7 @@ class TestPCA(object): + """ Test the PCA class """ def setUp(self): self.u = MDAnalysis.Universe(PDB, XTC) self.pca = pca.PCA(self.u.atoms, select='backbone and name CA',