diff --git a/package/MDAnalysis/analysis/align.py b/package/MDAnalysis/analysis/align.py index 490b3bb2e23..1593cd4f919 100644 --- a/package/MDAnalysis/analysis/align.py +++ b/package/MDAnalysis/analysis/align.py @@ -338,7 +338,6 @@ def _fit_to(mobile_coordinates, ref_coordinates, mobile_atoms, """ R, min_rmsd = rotation_matrix(mobile_coordinates, ref_coordinates, weights=weights) - mobile_atoms.translate(-mobile_com) mobile_atoms.rotate(R) mobile_atoms.translate(ref_com) @@ -629,9 +628,9 @@ def __init__(self, mobile, reference, select='all', filename=None, """ select = rms.process_selection(select) - self.ref_atoms = reference.select_atoms(*select['reference']) - self.mobile_atoms = mobile.select_atoms(*select['mobile']) - if in_memory or isinstance(mobile.trajectory, MemoryReader): + self.ref_atoms = reference.select_atoms(*select['reference']).atoms + self.mobile_atoms = mobile.select_atoms(*select['mobile']).atoms + if in_memory or not isinstance(mobile.trajectory, MemoryReader): mobile.transfer_to_memory() filename = None logger.info("Moved mobile trajectory to in-memory representation") diff --git a/package/MDAnalysis/analysis/clustering/cluster_methods.pyx b/package/MDAnalysis/analysis/clustering/cluster_methods.pyx index b95f37373e0..a5422815cd1 100644 --- a/package/MDAnalysis/analysis/clustering/cluster_methods.pyx +++ b/package/MDAnalysis/analysis/clustering/cluster_methods.pyx @@ -80,14 +80,12 @@ def affinity_propagation(similarity, preference, float lam, int max_iter, int co similarity[np.diag_indices(cn)] = preference indices = np.tril_indices(cn) - - # print(similarity) - # print(np.tril(similarity)) + cdef np.ndarray[np.float32_t, ndim=1] sim = np.ravel(similarity[indices]).astype(np.float32) # sim = np.ascontiguousarray(sim) - print(similarity[-2:]) + print(sim) cdef np.ndarray[long, ndim=1] clusters = np.zeros((cn), dtype=long) diff --git a/package/MDAnalysis/analysis/diffusionmap.py b/package/MDAnalysis/analysis/diffusionmap.py index 428ec8b404e..84fbf384822 100644 --- a/package/MDAnalysis/analysis/diffusionmap.py +++ b/package/MDAnalysis/analysis/diffusionmap.py @@ -231,13 +231,12 @@ def __init__(self, u, select='all', metric=rmsd, cutoff=1E0-5, Show detailed progress of the calculation if set to ``True``; the default is ``False``. """ - self._u = u + self.atoms = u.select_atoms(select).atoms + self._u = self.atoms.universe traj = self._u.trajectory # remember that this must be called before referencing self.n_frames super(DistanceMatrix, self).__init__(self._u.trajectory, **kwargs) - - self.atoms = self._u.select_atoms(select).atoms self._metric = metric self._cutoff = cutoff self._weights = weights diff --git a/package/MDAnalysis/analysis/encore/covariance.py b/package/MDAnalysis/analysis/encore/covariance.py index 349447e4f0d..2e6503ae57d 100644 --- a/package/MDAnalysis/analysis/encore/covariance.py +++ b/package/MDAnalysis/analysis/encore/covariance.py @@ -62,12 +62,14 @@ def ml_covariance_estimator(coordinates, reference_coordinates=None): else: # Normal covariance calculation: distance to the average coordinates_offset = coordinates - np.average(coordinates, axis=0) + # print(coordinates_offset[-10:]) # Calculate covariance manually coordinates_cov = np.zeros((coordinates.shape[1], coordinates.shape[1])) for frame in coordinates_offset: coordinates_cov += np.outer(frame, frame) + # print("cov", coordinates_cov[:10]) coordinates_cov /= coordinates.shape[0] return coordinates_cov @@ -121,6 +123,7 @@ def shrinkage_covariance_estimator( coordinates, # Call maximum likelihood estimator (note the additional column) sample = ml_covariance_estimator(np.hstack([x, xmkt[:, np.newaxis]]), 0) sample *= (t-1)/float(t) + # print("sample", sample[:10]) # Split covariance matrix into components @@ -131,6 +134,7 @@ def shrinkage_covariance_estimator( coordinates, # Prior prior = np.outer(covmkt, covmkt)/varmkt prior[np.ma.make_mask(np.eye(n))] = np.diag(sample) + # If shrinkage parameter is not set, estimate it if shrinkage_parameter is None: @@ -162,7 +166,10 @@ def shrinkage_covariance_estimator( coordinates, # Shrinkage constant k = (p-r)/c shrinkage_parameter = max(0, min(1, k/float(t))) - + # print(prior[0]) + # print("shrinkage_parameter", shrinkage_parameter) + # raise ValueError() + # calculate covariance matrix sigma = shrinkage_parameter*prior+(1-shrinkage_parameter)*sample @@ -222,6 +229,7 @@ def covariance_matrix(ensemble, reference_coordinates = reference_coordinates.flatten() sigma = estimator(coordinates, reference_coordinates) + # print(sigma[:5, :5]) # Optionally correct with weights if weights is not None: diff --git a/package/MDAnalysis/analysis/encore/similarity.py b/package/MDAnalysis/analysis/encore/similarity.py index 7d0db4331c9..7b9a330c7ce 100644 --- a/package/MDAnalysis/analysis/encore/similarity.py +++ b/package/MDAnalysis/analysis/encore/similarity.py @@ -286,8 +286,6 @@ def harmonic_ensemble_similarity(sigma1, # Difference between average vectors d_avg = x1 - x2 - - # Distance measure trace = np.trace(np.dot(sigma1, sigma2_inv) + np.dot(sigma2, sigma1_inv) @@ -296,16 +294,13 @@ def harmonic_ensemble_similarity(sigma1, d_hes = 0.25 * (np.dot(np.transpose(d_avg), np.dot(sigma1_inv + sigma2_inv, d_avg)) + trace) - print(d_avg.shape) - print((sigma1_inv + sigma2_inv).shape) - print((np.dot(sigma1_inv + sigma2_inv, - d_avg)).shape) - print(np.dot(sigma1, sigma2_inv).shape) - print((np.dot(np.transpose(d_avg), - np.dot(sigma1_inv + sigma2_inv, - d_avg)).shape)) + # print((sigma1_inv + sigma2_inv)[:10]) + print((np.dot(sigma1, sigma2_inv) + + np.dot(sigma2, sigma1_inv))[:10]) print(trace) - raise ValueError() + # print(sigma1[:5, :5]) + # print(sigma2[:5, :5]) + # print((sigma1_inv + sigma2_inv).dtype) return d_hes @@ -562,7 +557,12 @@ def dimred_ensemble_similarity(kde1, resamples1, kde2, resamples2, ln_P1P2_exp_P1 = np.average(np.log( 0.5 * (kde1.evaluate(resamples1) + kde2.evaluate(resamples1)))) ln_P1P2_exp_P2 = np.average(np.log( - 0.5 * (kde1.evaluate(resamples2) + kde2.evaluate(resamples2)))) + 0.5 * (kde1.evaluate(resamples2) + + kde2.evaluate(resamples2)))) + + print(ln_P1_exp_P1) + print(ln_P2_exp_P2) + print(ln_P1P2_exp_P1+ln_P1P2_exp_P2) return 0.5 * ( ln_P1_exp_P1 - ln_P1P2_exp_P1 + ln_P2_exp_P2 - ln_P1P2_exp_P2) @@ -952,6 +952,7 @@ def hes(ensembles, sigma2=sigmas[j]) values[i, j] = value values[j, i] = value + # Save details as required details = {} @@ -1495,6 +1496,7 @@ def dres(ensembles, if allow_collapsed_result and not hasattr(dimensionality_reduction_method, '__iter__'): values = values[0] + raise ValueError() return values, details diff --git a/package/MDAnalysis/analysis/reencore/similarity.py b/package/MDAnalysis/analysis/reencore/similarity.py index dbf422d0acc..a50681c55e9 100644 --- a/package/MDAnalysis/analysis/reencore/similarity.py +++ b/package/MDAnalysis/analysis/reencore/similarity.py @@ -4,7 +4,8 @@ from scipy.stats import gaussian_kde from ...core.ensemble import Ensemble -from .. import align +from ..align import AlignTraj +from ..pca import PCA from . import utils def clusters_to_indices(clusters, outlier_label=-1, outlier_index=-1): @@ -38,39 +39,50 @@ def ces(ensemble, clusters, outlier_label=-1, outlier_index=-1): outlier_i = n_cl # marker for new outlier columns for row, data in zip(frames_per_cl, ensemble.split_array(i_labels)): labels_, counts_ = np.unique(data, return_counts=True) + # treat outliers as individual clusters if labels_[0] == outlier_index: n_outlier_ = counts_[0] row[outlier_i:outlier_i + n_outlier_] = 1 outlier_i += n_outlier_ labels_ = labels_[1:] counts_ = counts_[1:] - # normalise over number of frames row[labels_] = counts_ - row /= len(data) - print(clusters[:20]) - print(clusters[20:]) - print(frames_per_cl) - - + row /= len(data) return utils.discrete_js_matrix(frames_per_cl) -def dres(ensemble, subspace, n_components=None, n_resample=1e3): +def dres(ensemble, subspace=PCA, n_components=3, n_resample=1000, + start=None, stop=None, step=None, **kwargs): + if isinstance(subspace, type): + subspace = subspace(ensemble, **kwargs) + if isinstance(subspace, PCA): + if not subspace._calculated: + subspace.run(start=start, stop=stop, step=step) + subspace = subspace.transform(subspace._atoms, + n_components=n_components, + start=start, stop=stop, + step=step) if n_components is not None: subspace = subspace[:, :n_components] - + n_u = ensemble.n_universes kdes = [gaussian_kde(x.T) for x in ensemble.split_array(subspace)] - resamples = np.concatenate([k.resample(size=n_resample) for k in kdes], - axis=1) # (n_dim, n_u*n_samples) - pdfs = np.array([k.evaluate(resamples) for k in kdes]) - pdfs = pdfs.reshape((n_u, n_resample)) - logpdfs = np.broadcast_to(np.log(pdfs).mean(axis=1), (n_u, n_u)) - - pdfs_ = np.broadcast_to(pdfs, (n_u, n_u, n_resample)) - sum_pdfs = pdfs_ + np.transpose(pdfs_, axes=(1, 0, 2)) - ln_pq_exp_pq = np.log(0.5 * sum_pdfs).mean(axis=-1) + resamples = [k.resample(size=n_resample) for k in kdes] + # resamples = np.concatenate([k.resample(size=n_resample) for k in kdes], + # axis=1) # (n_dim, n_u*n_samples) + # pdfs = np.array([k.evaluate(resamples) for k in kdes]) + # pdfs = np.array(np.split(pdfs, n_u, axis=1)) #pdfs.reshape((n_u, n_u, n_resample)) + pdfs = np.zeros((n_u, n_u, n_resample)) + for i, k in enumerate(kdes): + pdfs[i] = [k.evaluate(x) for x in resamples] + logpdfs = np.log(pdfs).mean(axis=-1) + pdfsT = pdfs.transpose((1, 0, 2)) + + sum_pdfs = pdfs + pdfsT + ln_pq_exp_pq = np.log(0.5 * (pdfs + pdfsT)).mean(axis=-1) + print(pdfs.shape) + print(logpdfs) return 0.5 * (logpdfs + logpdfs.T - ln_pq_exp_pq - ln_pq_exp_pq.T) @@ -116,7 +128,8 @@ def hes(ensemble, weights="mass", estimator="shrinkage", align=False, ensemble.transfer_to_memory() if align: - align.AlignTraj(ensemble, ensemble, weights=weights[0]).run() + AlignTraj(ensemble, ensemble, weights=weights[0], + in_memory=True).run() frames = [u.trajectory.timeseries(ag, order="fac") for ag, u in zip(ensemble._ags, ensemble.universes)] @@ -136,24 +149,41 @@ def hes(ensemble, weights="mass", estimator="shrinkage", align=False, for i, (coords, w) in enumerate(zip(frames[s], weights3)): avgs[s, i] = coords.mean(axis=0).flatten() cov = estimator(coords.reshape(len(coords), -1)) + # print(cov[:5, :5]) try: cov = np.dot(w, np.dot(cov, w)) except ValueError: raise ValueError("weights dimensions don't match selected atoms") + # print(cov[:5, :5]) covs[s, i] = cov inv_covs = np.zeros((n_s, n_u, n_u, n_a3, n_a3)) for i, sub in enumerate(covs): for j, arr in enumerate(sub): inv_covs[i, j] = np.linalg.pinv(arr[0]) + # if j == 0: + # print(arr[0]) + # print(inv_covs[i, j]) + # print(avgs.shape) diff = avgs - avgs.transpose((0, 2, 1, 3)) - cov_prod = covs @ inv_covs + + # print(avgs[0, 0, 0][:10]) + # print(avgs[0, 1, 0][:10]) + # print((avgs[0, 0, 0] - avgs[0, 1, 0])[:10]) + + # print(diff[0][0, 1][:20]) + + cov_prod = covs @ inv_covs.transpose((0, 2, 1, 3, 4)) cov_prod += cov_prod.transpose((0, 2, 1, 3, 4)) - trace = np.trace(cov_prod, axis1=-2, axis2=-1) - (2 * n_a3) + # print(cov_prod[0, 0, 1, : 10]) + # trace = np.trace(cov_prod - 2 * np.identity(n_a3), axis1=-1, axis2=-2) + trace = np.trace(cov_prod, axis1=-1, axis2=-2) - (2 * n_a3) + # print(trace) inv_cov_ = inv_covs + inv_covs.transpose((0, 2, 1, 3, 4)) - prod = np.einsum('ijklm,ijkl->ijkl', inv_covs, diff) + # print(inv_cov_.dtype) + prod = np.einsum('ijklm,ijkl->ijkl', inv_cov_, diff) similarity = np.einsum('ijkl,ijkl->ijk', diff, prod) similarity = 0.25 * (similarity + trace) if estimate_error: diff --git a/package/MDAnalysis/analysis/reencore/utils.py b/package/MDAnalysis/analysis/reencore/utils.py index c1aa143fdcd..ed64c4455e0 100644 --- a/package/MDAnalysis/analysis/reencore/utils.py +++ b/package/MDAnalysis/analysis/reencore/utils.py @@ -84,8 +84,11 @@ def discrete_js_matrix(matrix): def max_likelihood_covariance(coordinates, center=True): if center: coordinates -= coordinates.mean(axis=0) - - cov = np.einsum('ij,ik->jk', coordinates, coordinates) + _, y = coordinates.shape + cov = np.zeros((y, y)) + for frame in coordinates: + cov += np.outer(frame, frame) + # np.einsum('ij,ik->jk', coordinates, coordinates, out=cov) cov /= coordinates.shape[0] return cov @@ -110,7 +113,7 @@ def shrinkage_covariance(coordinates, shrinkage_parameter=None): if shrinkage_parameter is None: c = np.linalg.norm(sample - prior, ord="fro") ** 2 - inv_t = 1/len(offset) + inv_t = 1/t y = offset ** 2 p = inv_t * (y.T @ y).sum() - (sample ** 2).sum() rdiag = inv_t * (y ** 2).sum() - np.trace(sample ** 2) @@ -126,7 +129,10 @@ def shrinkage_covariance(coordinates, shrinkage_parameter=None): k = (p - r) / c shrinkage_parameter = max(0, min(1, k * inv_t)) - return shrinkage_parameter * prior + (1 - shrinkage_parameter) * sample + # print("shrinkage_parameter", shrinkage_parameter) + + cov = shrinkage_parameter * prior + (1 - shrinkage_parameter) * sample + return cov def get_bootstrap_frames(frames, n_samples=100): diff --git a/package/MDAnalysis/core/ensemble.py b/package/MDAnalysis/core/ensemble.py index 0f841879ce7..0c5e10de377 100644 --- a/package/MDAnalysis/core/ensemble.py +++ b/package/MDAnalysis/core/ensemble.py @@ -32,9 +32,27 @@ class AtomsWrapper: def __init__(self, ensemble): self.ensemble = ensemble + self.universe = ensemble + + def __len__(self): + return len(self.ensemble._atoms) def __getattr__(self, attr): - return getattr(self.ensemble._atoms, attr) + try: + return getattr(self.ensemble._atoms, attr) + except AttributeError: + return getattr(self.ensemble, attr) + + def select_atoms(self, *args, **kwargs): + return self.ensemble.select_atoms(*args, **kwargs).atoms + + @property + def positions(self): + return self.ensemble._atoms.positions + + @positions.setter + def positions(self, value): + self.ensemble._atoms.positions = value class Ensemble(object): @@ -42,7 +60,10 @@ class Ensemble(object): Should be able to plug this into any AnalysisBase class. """ - def __init__(self, universes, select=None, labels=None, frames=None): + def __init__(self, universes, select=None, labels=None, frames=None, + _ts_u=0, links=None): + self.ensemble = self + self.universe = self # set universe info universes = util.asiterable(universes) self.atoms = AtomsWrapper(self) @@ -68,6 +89,7 @@ def __init__(self, universes, select=None, labels=None, frames=None): raise ValueError(err) self._ag_indices = [ag.ix for ag in self._ags] + self._ags = tuple(self._ags) self.n_atoms = len(self._ags[0]) # label universes with filename or user-given name @@ -100,7 +122,18 @@ def __init__(self, universes, select=None, labels=None, frames=None): # pretend to be Universe self.trajectory = self - self._ts_u = 0 + self._ts_u = _ts_u + + if links is None: + links = set() + self.links = links + self.links.add(self) + + def __eq__(self, other): + return self._ags == other._ags + + def __hash__(self): + return hash(id(self)) def __getattr__(self, attr): try: @@ -113,17 +146,19 @@ def __getattr__(self, attr): raise AttributeError(f"{type(self).__name__} does not have " f"attribute {attr}") + def __repr__(self): + return (f"<{type(self).__name__} with {self.n_atoms} atoms, " + f"{len(self)} frames, frame={self.frame} at {hex(id(self))}>") + def __len__(self): return self.n_frames def __getitem__(self, i): if isinstance(i, (int, np.integer)): n_u, n_frame = self._get_relative_frame(i) - u = self.universes[n_u] - u.trajectory[n_frame] - self._ts_u = n_u + self._set_frame(ts_u=n_u, n_frame=n_frame) elif isinstance(i, str): - self._ts_u = self._universe_labels[i] + self._set_frame(ts_u=self._universe_labels[i]) elif util.iterable(i) or isinstance(i, slice): return self._slice_universe(i) else: @@ -131,11 +166,65 @@ def __getitem__(self, i): 'integers, slices, arrays, or string labels') return self.ts + def _set_frame(self, ts_u=None, n_frame=None): + if ts_u is not None: + for ens in self.links: + ens._ts_u = ts_u + + if n_frame is not None: + self._universe.trajectory[n_frame] + for ens in self.links: + ens._universe.trajectory[n_frame] + + @property + def frame(self): + ix = self._universe_frames + u = self.universes[self._ts_u] + u_ix = np.where(ix[:, 0] == self._ts_u)[0] + frame = u.trajectory.frame + fr_ix = np.where(ix[:, 1] == frame)[0] + try: + return u_ix[np.in1d(u_ix, fr_ix)][0] + except IndexError: + return None + + @property + def positions(self): + return self._atoms.positions + + @positions.setter + def positions(self, value): + self._atoms.positions = value + + def __iter__(self): for n_u, f in self._universe_frames: - self._ts_u = n_u - self.universes[n_u].trajectory[f] + self._set_frame(ts_u=n_u, n_frame=f) yield self.ts + self._set_frame(ts_u=0) + + def timeseries(self, atomgroup=None, start=None, stop=None, step=None, + frames=None, order="afc"): + n_atoms = self.n_atoms + frames = np.array(self._prepare_frames(start=start, stop=stop, + step=step, frames=frames)) + n_frames = len(frames) + fac = np.zeros((n_frames, n_atoms, 3)) + for i, frame in enumerate(frames): + self[frame] + np.copyto(fac[i], self.atoms.positions) + if order == "afc": + np.swapaxes(fac, 0, 1) + elif order == "caf": + np.swapaxes(fac, 2, 0) + elif order == "fca": + np.swapaxes(fac, 2, 1) + elif order == "cfa": + fac = fac.transpose((2, 0, 1)) + elif order == "acf": + fac = fac.transpose((1, 2, 0)) + return fac + @property def filename(self): @@ -150,13 +239,9 @@ def ts(self): return self.universes[self._ts_u].trajectory.ts @property - def universe(self): + def _universe(self): return self.universes[self._ts_u] - @property - def positions(self): - return self.atoms.positions - check_slice_indices = ProtoReader.check_slice_indices def _set_frames(self, frames=None): @@ -206,8 +291,7 @@ def iterate_over_atomgroups(self, start=None, stop=None, step=None, for i in frames: n_u, f = self._universe_frames[i] - self._ts_u = n_u - self.universe.trajectory[f] + self._set_frame(n_u, f) yield self._ags[n_u] def iterate_over_universes(self, start=None, stop=None, step=None, @@ -215,7 +299,8 @@ def iterate_over_universes(self, start=None, stop=None, step=None, frames_ = self._prepare_frames_by_universe(start=start, stop=stop, step=step, frames=frames) for ag, label, fr in zip(self._ags, self.labels, frames_): - yield type(self)([ag], labels=[label], frames=fr) + yield type(self)([ag], labels=[label], frames=fr, _ts_u=self._ts_u, + links=self.links) def _get_relative_frame(self, i): @@ -234,14 +319,16 @@ def _get_universe_from_frame(self, i): def _slice_universe(self, sl): frames = self.frames[sl] - return type(self)(self._ags, labels=self.labels, frames=frames) + return type(self)(self._ags, labels=self.labels, frames=frames, + _ts_u=self._ts_u, links=self.links) def select_atoms(self, *sel): sel = [s for s in sel if s is not None] ags = self._ags if sel: ags = [ag.select_atoms(*sel) for ag in ags] - return type(self)(ags, select=None, labels=self.labels) + return type(self)(ags, select=None, labels=self.labels, + _ts_u=self._ts_u, links=self.links) def split_array(self, arr): for i in range(self.n_universes): diff --git a/testsuite/MDAnalysisTests/analysis/test_encore.py b/testsuite/MDAnalysisTests/analysis/test_encore.py index 9285be4bc0b..4f22eec8098 100644 --- a/testsuite/MDAnalysisTests/analysis/test_encore.py +++ b/testsuite/MDAnalysisTests/analysis/test_encore.py @@ -250,23 +250,23 @@ def ens2(self, ens2_template): # result_value_custom = results_custom[0, 1] # assert_almost_equal(result_value, result_value_custom) -# def test_hes_align(self, ens1, ens2): -# # This test is massively sensitive! -# # Get 5260 when masses were float32? -# results, details = encore.hes([ens1, ens2], align=True) -# result_value = results[0,1] -# expected_value = 2047.05 -# assert_almost_equal(result_value, expected_value, decimal=-3, -# err_msg="Unexpected value for Harmonic Ensemble Similarity: {0:f}. Expected {1:f}.".format(result_value, expected_value)) - - def test_ces_to_self(self, ens1): - results, details = \ - encore.ces([ens1, ens1], - clustering_method=encore.AffinityPropagationNative(preference = -3.0)) - result_value = results[0,1] - expected_value = 0. - assert_almost_equal(result_value, expected_value, - err_msg="ClusteringEnsemble Similarity to itself not zero: {0:f}".format(result_value)) + # def test_hes_align(self, ens1, ens2): + # # This test is massively sensitive! + # # Get 5260 when masses were float32? + # results, details = encore.hes([ens1, ens2], align=True) + # result_value = results[0,1] + # expected_value = 2047.05 + # assert_almost_equal(result_value, expected_value, decimal=-3, + # err_msg="Unexpected value for Harmonic Ensemble Similarity: {0:f}. Expected {1:f}.".format(result_value, expected_value)) + + # def test_ces_to_self(self, ens1): + # results, details = \ + # encore.ces([ens1, ens1], + # clustering_method=encore.AffinityPropagationNative(preference = -3.0)) + # result_value = results[0,1] + # expected_value = 0. + # assert_almost_equal(result_value, expected_value, + # err_msg="ClusteringEnsemble Similarity to itself not zero: {0:f}".format(result_value)) # def test_ces(self, ens1, ens2): # results, details = encore.ces([ens1, ens2]) @@ -275,12 +275,12 @@ def test_ces_to_self(self, ens1): # assert_almost_equal(result_value, expected_value, decimal=2, # err_msg="Unexpected value for Cluster Ensemble Similarity: {0:f}. Expected {1:f}.".format(result_value, expected_value)) -# def test_dres_to_self(self, ens1): -# results, details = encore.dres([ens1, ens1]) -# result_value = results[0,1] -# expected_value = 0. -# assert_almost_equal(result_value, expected_value, decimal=2, -# err_msg="Dim. Reduction Ensemble Similarity to itself not zero: {0:f}".format(result_value)) + def test_dres_to_self(self, ens1): + results, details = encore.dres([ens1, ens1]) + result_value = results[0,1] + expected_value = 0. + assert_almost_equal(result_value, expected_value, decimal=2, + err_msg="Dim. Reduction Ensemble Similarity to itself not zero: {0:f}".format(result_value)) # def test_dres(self, ens1, ens2): # results, details = encore.dres([ens1, ens2], select="name CA and resnum 1-10") diff --git a/testsuite/MDAnalysisTests/analysis/test_reencore.py b/testsuite/MDAnalysisTests/analysis/test_reencore.py index 730598243bd..66a385f039e 100644 --- a/testsuite/MDAnalysisTests/analysis/test_reencore.py +++ b/testsuite/MDAnalysisTests/analysis/test_reencore.py @@ -21,7 +21,7 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # import MDAnalysis as mda -from MDAnalysis.analysis import encore, reencore, align +from MDAnalysis.analysis import encore, reencore, align, pca from MDAnalysis.analysis.diffusionmap import DistanceMatrix from MDAnalysis.analysis.clustering import Clusters, methods @@ -46,7 +46,7 @@ @pytest.fixture() def data(): - return np.arange(24, dtype=float).reshape((4, 6)) + return np.arange(24, dtype=np.float64).reshape((4, 6)) def test_max_likelihood_estimator(data): @@ -115,49 +115,58 @@ def ensemble1_aligned(self): a = align.AlignTraj(ens, ens, in_memory=True).run() return ens - @pytest.fixture() - def rmsd_mat(self, ensemble_aligned): - dist_mat = DistanceMatrix(ensemble_aligned, select="name CA").run() - return dist_mat.dist_matrix - - @pytest.fixture() - def rmsd_mat1(self, ensemble1_aligned): - dist_mat = DistanceMatrix(ensemble1_aligned).run() - return dist_mat.dist_matrix - def test_affinity_propagation(self, ens1): dist_mat = encore.get_distance_matrix(ens1).as_array() - # u = mda.Universe(PSF, DCD) - # u.transfer_to_memory(step=5) - # a = align.AlignTraj(u, u, select_atoms="name CA", in_memory=True).run() - # dist_mat = DistanceMatrix(u, select="name CA").run().dist_matrix clusters = Clusters(methods.AffinityPropagation) clusters.run(-dist_mat) assert len(clusters.cluster_indices) == 7 + def test_hes_to_self(self, ensemble1): + result_value = reencore.hes(ensemble1)[0, 1] + assert_almost_equal(result_value, 0, + err_msg=f"hes() to itself not zero: {result_value}") + + def test_hes(self, ensemble): + result = reencore.hes(ensemble, weights='mass')[0, 1] + min_bound = 1E5 + assert result > min_bound, "Unexpected value for Harmonic Ensemble Similarity" + + def test_hes_align(self, ens1, ens2): + ensemble = mda.Ensemble([ens1, ens2]).select_atoms("name CA") + result_value = reencore.hes(ensemble, align=True)[0, 1] + old, _ = encore.hes(ensemble.universes) + assert_almost_equal(result_value, old[0, 1], decimal=-2) - def test_ces_to_self(self, ensemble1, rmsd_mat1): + def test_ces_to_self(self, ensemble1_aligned): + dm = DistanceMatrix(ensemble1_aligned, select="name CA").run() + rmsd_mat1 = dm.dist_matrix + dm = DistanceMatrix(ensemble1_aligned).run() clusters = Clusters(methods.AffinityPropagation(preference=-3.0)) clusters.run(-rmsd_mat1) - result_value = reencore.ces(ensemble1, clusters)[0, 1] + result_value = reencore.ces(ensemble1_aligned, clusters)[0, 1] assert_almost_equal(result_value, 0, err_msg=f"ces() to itself not zero: {result_value}") - def test_ces(self, ensemble, rmsd_mat): - clusters = Clusters(methods.AffinityPropagation) - print(rmsd_mat[-2:]) + def test_ces_rmsd_enc(self, ensemble_aligned): + rmsd_mat_enc = encore.get_distance_matrix(ensemble_aligned).as_array() + clusters = Clusters(methods.AffinityPropagation()) + clusters.run(-rmsd_mat_enc) + result_value = reencore.ces(ensemble_aligned, clusters)[0, 1] + assert_almost_equal(result_value, 0.51, decimal=2, + err_msg=f"unexpected value") + + def test_ces(self, ensemble_aligned): + dm = DistanceMatrix(ensemble_aligned, select="name CA").run() + rmsd_mat = dm.dist_matrix + clusters = Clusters(methods.AffinityPropagation()) clusters.run(-rmsd_mat) - result_value = reencore.ces(ensemble, clusters)[0, 1] - assert_almost_equal(result_value, 0.51, + result_value = reencore.ces(ensemble_aligned, clusters)[0, 1] + assert_almost_equal(result_value, 0.69, decimal=2, err_msg=f"unexpected value") + def test_dres_to_self(self, ensemble1_aligned): + result = reencore.dres(ensemble1_aligned)[0, 1] + assert_almost_equal(result, 0) + + - def test_hes_to_self(self, ensemble1): - result_value = reencore.hes(ensemble1)[0, 1] - assert_almost_equal(result_value, 0, - err_msg=f"hes() to itself not zero: {result_value}") - - def test_hes(self, ensemble): - result = reencore.hes(ensemble, weights='mass')[0, 1] - min_bound = 1E5 - assert result > min_bound, "Unexpected value for Harmonic Ensemble Similarity" diff --git a/testsuite/MDAnalysisTests/core/test_ensemble.py b/testsuite/MDAnalysisTests/core/test_ensemble.py index 067863b52f6..74d35a3ea92 100644 --- a/testsuite/MDAnalysisTests/core/test_ensemble.py +++ b/testsuite/MDAnalysisTests/core/test_ensemble.py @@ -83,6 +83,7 @@ def test_analysis_multiple_in_memory(self, u1, u3): ens.transfer_to_memory() dist_mat = diffusionmap.DistanceMatrix(ens).run() dm = dist_mat.dist_matrix + # print(dm[-2:]) assert dm[-1, -2] != 0