From 4ad15a3db67bb1e0e578a8a835a5c9fd3b045ffd Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sat, 28 May 2016 11:31:08 +0200 Subject: [PATCH 01/25] [thermo] adding assertion messages to util --- pyemma/thermo/util/util.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index 427ebf9f5..43826f827 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -88,8 +88,8 @@ def get_averaged_bias_matrix(bias_sequences, dtrajs, nstates=None): def _ensure_umbrella_center(candidate, dimension): if isinstance(candidate, (_np.ndarray)): - assert candidate.ndim == 1 - assert candidate.shape[0] == dimension + assert candidate.ndim == 1, str(candidate) + " has ndim %d" % candidate.ndim + assert candidate.shape[0] == dimension, str(candidate) + "has shape[0] %d" % candidate.shape[0] return candidate.astype(_np.float64) elif types.is_int(candidate) or types.is_float(candidate): return candidate * _np.ones(shape=(dimension,), dtype=_np.float64) @@ -98,9 +98,9 @@ def _ensure_umbrella_center(candidate, dimension): def _ensure_force_constant(candidate, dimension): if isinstance(candidate, (_np.ndarray)): - assert candidate.shape[0] == dimension + assert candidate.shape[0] == dimension, str(candidate) + " has shape[0] %d" % candidate.shape[0] if candidate.ndim == 2: - assert candidate.shape[1] == dimension + assert candidate.shape[1] == dimension, str(candidate) + " has shape[1] %d" % candidate.shape[1] return candidate.astype(_np.float64) elif candidate.ndim == 1: return _np.diag(candidate).astype(dtype=_np.float64) @@ -155,8 +155,8 @@ def _get_umbrella_sampling_parameters( umbrella_centers = _np.array(umbrella_centers, dtype=_np.float64) force_constants = _np.array(force_constants, dtype=_np.float64) if kT is not None: - assert isinstance(kT, (int, float)) - assert kT > 0.0 + assert isinstance(kT, (int, float)), "kT has wrong type:" + str(type(kT)) + assert kT > 0.0, "non-positive kT: %f" % kT force_constants /= kT return ttrajs, umbrella_centers, force_constants, unbiased_state @@ -340,5 +340,5 @@ def assign_unbiased_state_label(memm_list, unbiased_state): if unbiased_state is None: return for memm in memm_list: - assert 0 <= unbiased_state < len(memm.models) + assert 0 <= unbiased_state < len(memm.models), "invalid state: " + str(unbiased_state) memm._unbiased_state = unbiased_state From 24a798be62501b18ca984c06544d2111d53c3e64 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sat, 28 May 2016 13:15:08 +0200 Subject: [PATCH 02/25] [thermo] adding unit tests for protected umbrella sampling util functions --- pyemma/thermo/tests/test_util.py | 189 +++++++++++++++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 pyemma/thermo/tests/test_util.py diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py new file mode 100644 index 000000000..ba369d806 --- /dev/null +++ b/pyemma/thermo/tests/test_util.py @@ -0,0 +1,189 @@ +# This file is part of PyEMMA. +# +# Copyright (c) 2016 Computational Molecular Biology Group, Freie Universitaet Berlin (GER) +# +# PyEMMA is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + +import unittest +import numpy as np +from pyemma.thermo.util.util import _ensure_umbrella_center +from pyemma.thermo.util.util import _ensure_force_constant + +# ================================================================================================== +# tests for protected umbrella sampling convenience functions +# ================================================================================================== + +class TestProtectedUmbrellaSamplingCenters(unittest.TestCase): + + def _assert_us_center(self, us_center, dimension): + self.assertIsInstance(us_center, np.ndarray) + self.assertTrue(us_center.dtype == np.float64) + self.assertTrue(us_center.ndim == 1) + self.assertTrue(us_center.shape[0] == dimension) + + def test_ensure_umbrella_center_from_scalar(self): + # dimension=1 + us_center = _ensure_umbrella_center(1.0, 1) + self._assert_us_center(us_center, 1) + np.testing.assert_array_equal(us_center, np.array([1.0], dtype=np.float64)) + # dimension=3 + us_center = _ensure_umbrella_center(1.0, 3) + self._assert_us_center(us_center, 3) + np.testing.assert_array_equal(us_center, np.array([1.0, 1.0, 1.0], dtype=np.float64)) + + def test_ensure_umbrella_center_from_tuple(self): + # dimension=1, type=tuple + us_center = _ensure_umbrella_center((1.0,), 1) + self._assert_us_center(us_center, 1) + np.testing.assert_array_equal(us_center, np.array([1.0], dtype=np.float64)) + # dimension=3, uniform + us_center = _ensure_umbrella_center((1.0, 1.0, 1.0), 3) + self._assert_us_center(us_center, 3) + np.testing.assert_array_equal(us_center, np.array([1.0, 1.0, 1.0], dtype=np.float64)) + # dimension=4, not uniform + us_center = _ensure_umbrella_center((1.0, 2.0, 3.0, 4.0), 4) + self._assert_us_center(us_center, 4) + np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) + # dimension=4x1, not uniform + us_center = _ensure_umbrella_center(((1.0, 2.0, 3.0, 4.0),), 4) + self._assert_us_center(us_center, 4) + np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) + + def test_ensure_umbrella_center_from_list(self): + # dimension=1 + us_center = _ensure_umbrella_center([1.0], 1) + self._assert_us_center(us_center, 1) + np.testing.assert_array_equal(us_center, np.array([1.0], dtype=np.float64)) + # dimension=3, uniform + us_center = _ensure_umbrella_center([1.0, 1.0, 1.0], 3) + self._assert_us_center(us_center, 3) + # dimension=4, not uniform + us_center = _ensure_umbrella_center([1.0, 2.0, 3.0, 4.0], 4) + self._assert_us_center(us_center, 4) + np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) + # dimension=4x1, not uniform + us_center = _ensure_umbrella_center([[1.0, 2.0, 3.0, 4.0],], 4) + self._assert_us_center(us_center, 4) + np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) + + def test_ensure_umbrella_center_from_ndarray(self): + # dimension=1 + us_center = _ensure_umbrella_center(np.array([1.0]), 1) + self._assert_us_center(us_center, 1) + np.testing.assert_array_equal(us_center, np.array([1.0], dtype=np.float64)) + # dimension=3, uniform + us_center = _ensure_umbrella_center(np.array([1.0, 1.0, 1.0]), 3) + self._assert_us_center(us_center, 3) + np.testing.assert_array_equal(us_center, np.array([1.0, 1.0, 1.0], dtype=np.float64)) + # dimension=4, not uniform + us_center = _ensure_umbrella_center(np.array([1.0, 2.0, 3.0, 4.0]), 4) + self._assert_us_center(us_center, 4) + np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) + # dimension=4x1, not uniform + us_center = _ensure_umbrella_center(np.array([[1.0, 2.0, 3.0, 4.0],]), 4) + self._assert_us_center(us_center, 4) + np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) + + +class TestProtectedUmbrellaSamplingForceMatrices(unittest.TestCase): + + def _assert_us_force_matrix(self, us_force_matrix, dimension): + self.assertIsInstance(us_force_matrix, np.ndarray) + self.assertTrue(us_force_matrix.dtype == np.float64) + self.assertTrue(us_force_matrix.ndim == 2) + self.assertTrue(us_force_matrix.shape[0] == dimension) + self.assertTrue(us_force_matrix.shape[1] == dimension) + + def test_ensure_umbrella_force_matrix_from_scalar(self): + # dimension=1 + us_force_matrix = _ensure_force_constant(1.0, 1) + self._assert_us_force_matrix(us_force_matrix, 1) + np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) + # dimension=2 + us_force_matrix = _ensure_force_constant(1.0, 2) + self._assert_us_force_matrix(us_force_matrix, 2) + np.testing.assert_array_equal( + us_force_matrix, np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.float64)) + + def test_ensure_umbrella_force_matrix_from_tuple(self): + # dimension=1 + us_force_matrix = _ensure_force_constant((1.0,), 1) + self._assert_us_force_matrix(us_force_matrix, 1) + np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) + # dimension=1x1 + us_force_matrix = _ensure_force_constant(((1.0,),), 1) + self._assert_us_force_matrix(us_force_matrix, 1) + np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) + # dimension=2, not uniform, diagonal + us_force_matrix = _ensure_force_constant((1.0, 2.0), 2) + self._assert_us_force_matrix(us_force_matrix, 2) + np.testing.assert_array_equal( + us_force_matrix, np.array([[1.0, 0.0], [0.0, 2.0]], dtype=np.float64)) + # dimension=2, not uniform, not diagonal + us_force_matrix = _ensure_force_constant(((1.0, 2.0), (3.0, 4.0)), 2) + self._assert_us_force_matrix(us_force_matrix, 2) + np.testing.assert_array_equal( + us_force_matrix, np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)) + + def test_ensure_umbrella_force_matrix_from_list(self): + # dimension=1 + us_force_matrix = _ensure_force_constant([1.0], 1) + self._assert_us_force_matrix(us_force_matrix, 1) + np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) + # dimension=1x1 + us_force_matrix = _ensure_force_constant([[1.0]], 1) + self._assert_us_force_matrix(us_force_matrix, 1) + np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) + # dimension=2, not uniform, diagonal + us_force_matrix = _ensure_force_constant([1.0, 2.0], 2) + self._assert_us_force_matrix(us_force_matrix, 2) + np.testing.assert_array_equal( + us_force_matrix, np.array([[1.0, 0.0], [0.0, 2.0]], dtype=np.float64)) + # dimension=2, not uniform, not diagonal + us_force_matrix = _ensure_force_constant([[1.0, 2.0], [3.0, 4.0]], 2) + self._assert_us_force_matrix(us_force_matrix, 2) + np.testing.assert_array_equal( + us_force_matrix, np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)) + + def test_ensure_umbrella_force_matrix_from_ndarray(self): + # dimension=1 + us_force_matrix = _ensure_force_constant(np.array([1.0]), 1) + self._assert_us_force_matrix(us_force_matrix, 1) + np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) + # dimension=1x1 + us_force_matrix = _ensure_force_constant(np.array([[1.0]]), 1) + self._assert_us_force_matrix(us_force_matrix, 1) + np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) + # dimension=2, not uniform, diagonal + us_force_matrix = _ensure_force_constant(np.array([1.0, 2.0]), 2) + self._assert_us_force_matrix(us_force_matrix, 2) + np.testing.assert_array_equal( + us_force_matrix, np.array([[1.0, 0.0], [0.0, 2.0]], dtype=np.float64)) + # dimension=2, not uniform, not diagonal + us_force_matrix = _ensure_force_constant(np.array([[1.0, 2.0], [3.0, 4.0]]), 2) + self._assert_us_force_matrix(us_force_matrix, 2) + np.testing.assert_array_equal( + us_force_matrix, np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)) + + + + + + + + + + + + From b22207d8bc4d043352b7511b5cb3a41d496a6349 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sat, 28 May 2016 13:35:14 +0200 Subject: [PATCH 03/25] [thermo] making util._ensure_umbrella_center() robust --- pyemma/thermo/util/util.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index 43826f827..724647843 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -87,14 +87,10 @@ def get_averaged_bias_matrix(bias_sequences, dtrajs, nstates=None): # ================================================================================================== def _ensure_umbrella_center(candidate, dimension): - if isinstance(candidate, (_np.ndarray)): - assert candidate.ndim == 1, str(candidate) + " has ndim %d" % candidate.ndim - assert candidate.shape[0] == dimension, str(candidate) + "has shape[0] %d" % candidate.shape[0] - return candidate.astype(_np.float64) - elif types.is_int(candidate) or types.is_float(candidate): - return candidate * _np.ones(shape=(dimension,), dtype=_np.float64) - else: - raise TypeError("unsupported type") + candidate = _np.asarray(candidate).astype(_np.float64).reshape((-1,)) + if candidate.shape[0] == 1 and dimension > 1: + return candidate[0] * _np.ones(shape=(dimension,), dtype=_np.float64) + return candidate def _ensure_force_constant(candidate, dimension): if isinstance(candidate, (_np.ndarray)): From 719c6ed429dddd17523d69dcd1895d77325e12ae Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sat, 28 May 2016 13:48:30 +0200 Subject: [PATCH 04/25] [thermo] making util._ensure_force_constant() robust --- pyemma/thermo/util/util.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index 724647843..7c08df5ac 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -93,19 +93,23 @@ def _ensure_umbrella_center(candidate, dimension): return candidate def _ensure_force_constant(candidate, dimension): - if isinstance(candidate, (_np.ndarray)): - assert candidate.shape[0] == dimension, str(candidate) + " has shape[0] %d" % candidate.shape[0] - if candidate.ndim == 2: - assert candidate.shape[1] == dimension, str(candidate) + " has shape[1] %d" % candidate.shape[1] - return candidate.astype(_np.float64) - elif candidate.ndim == 1: - return _np.diag(candidate).astype(dtype=_np.float64) - else: - raise TypeError("usupported shape") - elif isinstance(candidate, (int, float)): - return candidate * _np.ones(shape=(dimension, dimension), dtype=_np.float64) + candidate = _np.asarray(candidate).astype(_np.float64) + if candidate.ndim == 0: + candidate = candidate * _np.ones(shape=(dimension,), dtype=_np.float64) + if candidate.ndim == 1: + assert candidate.shape[0] == dimension, \ + str(candidate) + " has shape[0] %d" % candidate.shape[0] + _candidate = _np.zeros(shape=(dimension, dimension), dtype=_np.float64) + for i, x in enumerate(candidate): + _candidate[i, i] = x + candidate = _candidate else: - raise TypeError("unsupported type") + assert candidate.ndim == 2, str(candidate) + " has wrong ndim %d" % candidate.ndim + assert candidate.shape[0] == dimension, \ + str(candidate) + " has shape[0] %d" % candidate.shape[0] + assert candidate.shape[1] == dimension, \ + str(candidate) + " has shape[1] %d" % candidate.shape[1] + return candidate def _get_umbrella_sampling_parameters( us_trajs, us_centers, us_force_constants, md_trajs=None, kT=None): From 38f0cead5e7bcc8a2408e2d0ca0c9a54bbf3487c Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sun, 29 May 2016 13:13:54 +0200 Subject: [PATCH 05/25] [thermo] adding unit tests for util._get_umbrella_sampling_parameters() --- pyemma/thermo/tests/test_util.py | 80 ++++++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index ba369d806..cd6429690 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/pyemma/thermo/tests/test_util.py @@ -19,6 +19,7 @@ import numpy as np from pyemma.thermo.util.util import _ensure_umbrella_center from pyemma.thermo.util.util import _ensure_force_constant +from pyemma.thermo.util.util import _get_umbrella_sampling_parameters # ================================================================================================== # tests for protected umbrella sampling convenience functions @@ -177,6 +178,85 @@ def test_ensure_umbrella_force_matrix_from_ndarray(self): us_force_matrix, np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)) +class TestProtectedUmbrellaSamplingParameters(unittest.TestCase): + + def _assert_parameters(self, + ttrajs, umbrella_centers, force_constants, unbiased_state, + ref_ttrajs, ref_umbrella_centers, ref_force_constants, ref_unbiased_state): + for ttraj, ref_ttraj in zip(ttrajs, ref_ttrajs): + np.testing.assert_array_equal(ttraj, ref_ttraj) + for center, ref_center in zip(umbrella_centers, ref_umbrella_centers): + np.testing.assert_array_equal(center, ref_center) + for force_constant, ref_force_constant in zip(force_constants, ref_force_constants): + np.testing.assert_array_equal(force_constant, ref_force_constant) + self.assertTrue(unbiased_state == ref_unbiased_state) + + def test_umbrella_sampling_data_1x0(self): + ref_umbrella_centers = [0.0, 1.0] + ref_force_constants = [1.0, 1.0] + us_trajs = [np.array([0.0, 0.1, 0.2]), np.array([0.9, 1.0, 1.1])] + # no md data + ttrajs, umbrella_centers, force_constants, unbiased_state = \ + _get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants) + self._assert_parameters( + ttrajs, umbrella_centers, force_constants, unbiased_state, + [np.array([0, 0, 0]), np.array([1, 1, 1])], + ref_umbrella_centers, ref_force_constants, None) + # add md data + md_trajs = [np.array([0.0, 0.5, 1.0])] + ttrajs, umbrella_centers, force_constants, unbiased_state = \ + _get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs) + self._assert_parameters( + ttrajs, umbrella_centers, force_constants, unbiased_state, + [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], + ref_umbrella_centers, ref_force_constants, 2) + # with kT parameter + ttrajs, umbrella_centers, force_constants, unbiased_state = \ + _get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) + self._assert_parameters( + ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, + [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], + ref_umbrella_centers, ref_force_constants, 2) + + def test_umbrella_sampling_data_1x1(self): + ref_umbrella_centers = [0.0, 1.0] + ref_force_constants = [1.0, 1.0] + us_trajs = [np.array([[0.0], [0.1], [0.2]]), np.array([[0.9], [1.0], [1.1]])] + ttrajs, umbrella_centers, force_constants, unbiased_state = \ + _get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants) + self._assert_parameters( + ttrajs, umbrella_centers, force_constants, unbiased_state, + [np.array([0, 0, 0]), np.array([1, 1, 1])], + ref_umbrella_centers, ref_force_constants, None) + # add md data + md_trajs = [np.array([[0.0], [0.5], [1.0]])] + ttrajs, umbrella_centers, force_constants, unbiased_state = \ + _get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs) + self._assert_parameters( + ttrajs, umbrella_centers, force_constants, unbiased_state, + [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], + ref_umbrella_centers, ref_force_constants, 2) + # with kT parameter + ttrajs, umbrella_centers, force_constants, unbiased_state = \ + _get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) + self._assert_parameters( + ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, + [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], + ref_umbrella_centers, ref_force_constants, 2) + + + + + + + + From 40976a6c33b90e622f75766cc0eef968414f99a4 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sun, 29 May 2016 13:27:40 +0200 Subject: [PATCH 06/25] [thermo] making util._get_umbrella_sampling_parameters() robust --- pyemma/thermo/util/util.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index 7c08df5ac..bfa3c7399 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -118,12 +118,22 @@ def _get_umbrella_sampling_parameters( ttrajs = [] nthermo = 0 unbiased_state = None - for i in range(len(us_trajs)): + dimension = None + for i, traj in enumerate(us_trajs): state = None + try: + _dimension = traj.shape[1] + except IndexError: + _dimension = 1 + if dimension is None: + dimension = _dimension + else: + assert dimension == _dimension, "trajectory %i has unmatching dimension %d!=%d" % ( + i, _dimension, dimension) this_center = _ensure_umbrella_center( - us_centers[i], us_trajs[i].shape[1]) + us_centers[i], dimension) this_force_constant = _ensure_force_constant( - us_force_constants[i], us_trajs[i].shape[1]) + us_force_constants[i], dimension) if _np.all(this_force_constant == 0.0): this_center *= 0.0 for j in range(nthermo): @@ -134,10 +144,10 @@ def _get_umbrella_sampling_parameters( if state is None: umbrella_centers.append(this_center.copy()) force_constants.append(this_force_constant.copy()) - ttrajs.append(nthermo * _np.ones(shape=(us_trajs[i].shape[0],), dtype=_np.intc)) + ttrajs.append(nthermo * _np.ones(shape=(traj.shape[0],), dtype=_np.intc)) nthermo += 1 else: - ttrajs.append(state * _np.ones(shape=(us_trajs[i].shape[0],), dtype=_np.intc)) + ttrajs.append(state * _np.ones(shape=(traj.shape[0],), dtype=_np.intc)) if md_trajs is not None: this_center = umbrella_centers[-1] * 0.0 this_force_constant = force_constants[-1] * 0.0 @@ -150,8 +160,8 @@ def _get_umbrella_sampling_parameters( force_constants.append(this_force_constant.copy()) unbiased_state = nthermo nthermo += 1 - for md_traj in md_trajs: - ttrajs.append(unbiased_state * _np.ones(shape=(md_traj.shape[0],), dtype=_np.intc)) + for traj in md_trajs: + ttrajs.append(unbiased_state * _np.ones(shape=(traj.shape[0],), dtype=_np.intc)) umbrella_centers = _np.array(umbrella_centers, dtype=_np.float64) force_constants = _np.array(force_constants, dtype=_np.float64) if kT is not None: From 40a8c6ddbee18b5e863d671510b3016da04a3c87 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sun, 29 May 2016 13:51:29 +0200 Subject: [PATCH 07/25] [thermo] minor change --- pyemma/thermo/tests/test_util.py | 72 ++++++++++++++++---------------- 1 file changed, 35 insertions(+), 37 deletions(-) diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index cd6429690..4fbf2ee98 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/pyemma/thermo/tests/test_util.py @@ -17,9 +17,7 @@ import unittest import numpy as np -from pyemma.thermo.util.util import _ensure_umbrella_center -from pyemma.thermo.util.util import _ensure_force_constant -from pyemma.thermo.util.util import _get_umbrella_sampling_parameters +import pyemma.thermo.util.util as util # ================================================================================================== # tests for protected umbrella sampling convenience functions @@ -35,64 +33,64 @@ def _assert_us_center(self, us_center, dimension): def test_ensure_umbrella_center_from_scalar(self): # dimension=1 - us_center = _ensure_umbrella_center(1.0, 1) + us_center = util._ensure_umbrella_center(1.0, 1) self._assert_us_center(us_center, 1) np.testing.assert_array_equal(us_center, np.array([1.0], dtype=np.float64)) # dimension=3 - us_center = _ensure_umbrella_center(1.0, 3) + us_center = util._ensure_umbrella_center(1.0, 3) self._assert_us_center(us_center, 3) np.testing.assert_array_equal(us_center, np.array([1.0, 1.0, 1.0], dtype=np.float64)) def test_ensure_umbrella_center_from_tuple(self): # dimension=1, type=tuple - us_center = _ensure_umbrella_center((1.0,), 1) + us_center = util._ensure_umbrella_center((1.0,), 1) self._assert_us_center(us_center, 1) np.testing.assert_array_equal(us_center, np.array([1.0], dtype=np.float64)) # dimension=3, uniform - us_center = _ensure_umbrella_center((1.0, 1.0, 1.0), 3) + us_center = util._ensure_umbrella_center((1.0, 1.0, 1.0), 3) self._assert_us_center(us_center, 3) np.testing.assert_array_equal(us_center, np.array([1.0, 1.0, 1.0], dtype=np.float64)) # dimension=4, not uniform - us_center = _ensure_umbrella_center((1.0, 2.0, 3.0, 4.0), 4) + us_center = util._ensure_umbrella_center((1.0, 2.0, 3.0, 4.0), 4) self._assert_us_center(us_center, 4) np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) # dimension=4x1, not uniform - us_center = _ensure_umbrella_center(((1.0, 2.0, 3.0, 4.0),), 4) + us_center = util._ensure_umbrella_center(((1.0, 2.0, 3.0, 4.0),), 4) self._assert_us_center(us_center, 4) np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) def test_ensure_umbrella_center_from_list(self): # dimension=1 - us_center = _ensure_umbrella_center([1.0], 1) + us_center = util._ensure_umbrella_center([1.0], 1) self._assert_us_center(us_center, 1) np.testing.assert_array_equal(us_center, np.array([1.0], dtype=np.float64)) # dimension=3, uniform - us_center = _ensure_umbrella_center([1.0, 1.0, 1.0], 3) + us_center = util._ensure_umbrella_center([1.0, 1.0, 1.0], 3) self._assert_us_center(us_center, 3) # dimension=4, not uniform - us_center = _ensure_umbrella_center([1.0, 2.0, 3.0, 4.0], 4) + us_center = util._ensure_umbrella_center([1.0, 2.0, 3.0, 4.0], 4) self._assert_us_center(us_center, 4) np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) # dimension=4x1, not uniform - us_center = _ensure_umbrella_center([[1.0, 2.0, 3.0, 4.0],], 4) + us_center = util._ensure_umbrella_center([[1.0, 2.0, 3.0, 4.0],], 4) self._assert_us_center(us_center, 4) np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) def test_ensure_umbrella_center_from_ndarray(self): # dimension=1 - us_center = _ensure_umbrella_center(np.array([1.0]), 1) + us_center = util._ensure_umbrella_center(np.array([1.0]), 1) self._assert_us_center(us_center, 1) np.testing.assert_array_equal(us_center, np.array([1.0], dtype=np.float64)) # dimension=3, uniform - us_center = _ensure_umbrella_center(np.array([1.0, 1.0, 1.0]), 3) + us_center = util._ensure_umbrella_center(np.array([1.0, 1.0, 1.0]), 3) self._assert_us_center(us_center, 3) np.testing.assert_array_equal(us_center, np.array([1.0, 1.0, 1.0], dtype=np.float64)) # dimension=4, not uniform - us_center = _ensure_umbrella_center(np.array([1.0, 2.0, 3.0, 4.0]), 4) + us_center = util._ensure_umbrella_center(np.array([1.0, 2.0, 3.0, 4.0]), 4) self._assert_us_center(us_center, 4) np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) # dimension=4x1, not uniform - us_center = _ensure_umbrella_center(np.array([[1.0, 2.0, 3.0, 4.0],]), 4) + us_center = util._ensure_umbrella_center(np.array([[1.0, 2.0, 3.0, 4.0],]), 4) self._assert_us_center(us_center, 4) np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) @@ -108,71 +106,71 @@ def _assert_us_force_matrix(self, us_force_matrix, dimension): def test_ensure_umbrella_force_matrix_from_scalar(self): # dimension=1 - us_force_matrix = _ensure_force_constant(1.0, 1) + us_force_matrix = util._ensure_force_constant(1.0, 1) self._assert_us_force_matrix(us_force_matrix, 1) np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) # dimension=2 - us_force_matrix = _ensure_force_constant(1.0, 2) + us_force_matrix = util._ensure_force_constant(1.0, 2) self._assert_us_force_matrix(us_force_matrix, 2) np.testing.assert_array_equal( us_force_matrix, np.array([[1.0, 0.0], [0.0, 1.0]], dtype=np.float64)) def test_ensure_umbrella_force_matrix_from_tuple(self): # dimension=1 - us_force_matrix = _ensure_force_constant((1.0,), 1) + us_force_matrix = util._ensure_force_constant((1.0,), 1) self._assert_us_force_matrix(us_force_matrix, 1) np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) # dimension=1x1 - us_force_matrix = _ensure_force_constant(((1.0,),), 1) + us_force_matrix = util._ensure_force_constant(((1.0,),), 1) self._assert_us_force_matrix(us_force_matrix, 1) np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) # dimension=2, not uniform, diagonal - us_force_matrix = _ensure_force_constant((1.0, 2.0), 2) + us_force_matrix = util._ensure_force_constant((1.0, 2.0), 2) self._assert_us_force_matrix(us_force_matrix, 2) np.testing.assert_array_equal( us_force_matrix, np.array([[1.0, 0.0], [0.0, 2.0]], dtype=np.float64)) # dimension=2, not uniform, not diagonal - us_force_matrix = _ensure_force_constant(((1.0, 2.0), (3.0, 4.0)), 2) + us_force_matrix = util._ensure_force_constant(((1.0, 2.0), (3.0, 4.0)), 2) self._assert_us_force_matrix(us_force_matrix, 2) np.testing.assert_array_equal( us_force_matrix, np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)) def test_ensure_umbrella_force_matrix_from_list(self): # dimension=1 - us_force_matrix = _ensure_force_constant([1.0], 1) + us_force_matrix = util._ensure_force_constant([1.0], 1) self._assert_us_force_matrix(us_force_matrix, 1) np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) # dimension=1x1 - us_force_matrix = _ensure_force_constant([[1.0]], 1) + us_force_matrix = util._ensure_force_constant([[1.0]], 1) self._assert_us_force_matrix(us_force_matrix, 1) np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) # dimension=2, not uniform, diagonal - us_force_matrix = _ensure_force_constant([1.0, 2.0], 2) + us_force_matrix = util._ensure_force_constant([1.0, 2.0], 2) self._assert_us_force_matrix(us_force_matrix, 2) np.testing.assert_array_equal( us_force_matrix, np.array([[1.0, 0.0], [0.0, 2.0]], dtype=np.float64)) # dimension=2, not uniform, not diagonal - us_force_matrix = _ensure_force_constant([[1.0, 2.0], [3.0, 4.0]], 2) + us_force_matrix = util._ensure_force_constant([[1.0, 2.0], [3.0, 4.0]], 2) self._assert_us_force_matrix(us_force_matrix, 2) np.testing.assert_array_equal( us_force_matrix, np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)) def test_ensure_umbrella_force_matrix_from_ndarray(self): # dimension=1 - us_force_matrix = _ensure_force_constant(np.array([1.0]), 1) + us_force_matrix = util._ensure_force_constant(np.array([1.0]), 1) self._assert_us_force_matrix(us_force_matrix, 1) np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) # dimension=1x1 - us_force_matrix = _ensure_force_constant(np.array([[1.0]]), 1) + us_force_matrix = util._ensure_force_constant(np.array([[1.0]]), 1) self._assert_us_force_matrix(us_force_matrix, 1) np.testing.assert_array_equal(us_force_matrix, np.array([[1.0]], dtype=np.float64)) # dimension=2, not uniform, diagonal - us_force_matrix = _ensure_force_constant(np.array([1.0, 2.0]), 2) + us_force_matrix = util._ensure_force_constant(np.array([1.0, 2.0]), 2) self._assert_us_force_matrix(us_force_matrix, 2) np.testing.assert_array_equal( us_force_matrix, np.array([[1.0, 0.0], [0.0, 2.0]], dtype=np.float64)) # dimension=2, not uniform, not diagonal - us_force_matrix = _ensure_force_constant(np.array([[1.0, 2.0], [3.0, 4.0]]), 2) + us_force_matrix = util._ensure_force_constant(np.array([[1.0, 2.0], [3.0, 4.0]]), 2) self._assert_us_force_matrix(us_force_matrix, 2) np.testing.assert_array_equal( us_force_matrix, np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)) @@ -197,7 +195,7 @@ def test_umbrella_sampling_data_1x0(self): us_trajs = [np.array([0.0, 0.1, 0.2]), np.array([0.9, 1.0, 1.1])] # no md data ttrajs, umbrella_centers, force_constants, unbiased_state = \ - _get_umbrella_sampling_parameters( + util._get_umbrella_sampling_parameters( us_trajs, ref_umbrella_centers, ref_force_constants) self._assert_parameters( ttrajs, umbrella_centers, force_constants, unbiased_state, @@ -206,7 +204,7 @@ def test_umbrella_sampling_data_1x0(self): # add md data md_trajs = [np.array([0.0, 0.5, 1.0])] ttrajs, umbrella_centers, force_constants, unbiased_state = \ - _get_umbrella_sampling_parameters( + util._get_umbrella_sampling_parameters( us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs) self._assert_parameters( ttrajs, umbrella_centers, force_constants, unbiased_state, @@ -214,7 +212,7 @@ def test_umbrella_sampling_data_1x0(self): ref_umbrella_centers, ref_force_constants, 2) # with kT parameter ttrajs, umbrella_centers, force_constants, unbiased_state = \ - _get_umbrella_sampling_parameters( + util._get_umbrella_sampling_parameters( us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) self._assert_parameters( ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, @@ -226,7 +224,7 @@ def test_umbrella_sampling_data_1x1(self): ref_force_constants = [1.0, 1.0] us_trajs = [np.array([[0.0], [0.1], [0.2]]), np.array([[0.9], [1.0], [1.1]])] ttrajs, umbrella_centers, force_constants, unbiased_state = \ - _get_umbrella_sampling_parameters( + util._get_umbrella_sampling_parameters( us_trajs, ref_umbrella_centers, ref_force_constants) self._assert_parameters( ttrajs, umbrella_centers, force_constants, unbiased_state, @@ -235,7 +233,7 @@ def test_umbrella_sampling_data_1x1(self): # add md data md_trajs = [np.array([[0.0], [0.5], [1.0]])] ttrajs, umbrella_centers, force_constants, unbiased_state = \ - _get_umbrella_sampling_parameters( + util._get_umbrella_sampling_parameters( us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs) self._assert_parameters( ttrajs, umbrella_centers, force_constants, unbiased_state, @@ -243,7 +241,7 @@ def test_umbrella_sampling_data_1x1(self): ref_umbrella_centers, ref_force_constants, 2) # with kT parameter ttrajs, umbrella_centers, force_constants, unbiased_state = \ - _get_umbrella_sampling_parameters( + util._get_umbrella_sampling_parameters( us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) self._assert_parameters( ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, From 8f067a6b4d7f222f6e767748ed1f0c956601b94a Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sun, 29 May 2016 13:56:28 +0200 Subject: [PATCH 08/25] [thermo] minor change --- pyemma/thermo/tests/test_util.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index 4fbf2ee98..324dd0a93 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/pyemma/thermo/tests/test_util.py @@ -189,7 +189,7 @@ def _assert_parameters(self, np.testing.assert_array_equal(force_constant, ref_force_constant) self.assertTrue(unbiased_state == ref_unbiased_state) - def test_umbrella_sampling_data_1x0(self): + def test_umbrella_sampling_parameters_1x0(self): ref_umbrella_centers = [0.0, 1.0] ref_force_constants = [1.0, 1.0] us_trajs = [np.array([0.0, 0.1, 0.2]), np.array([0.9, 1.0, 1.1])] @@ -219,7 +219,7 @@ def test_umbrella_sampling_data_1x0(self): [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) - def test_umbrella_sampling_data_1x1(self): + def test_umbrella_sampling_parameters_1x1(self): ref_umbrella_centers = [0.0, 1.0] ref_force_constants = [1.0, 1.0] us_trajs = [np.array([[0.0], [0.1], [0.2]]), np.array([[0.9], [1.0], [1.1]])] @@ -262,6 +262,3 @@ def test_umbrella_sampling_data_1x1(self): - - - From 3350f7ed66b2b228349d3b846ff17dafbf9a2756 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sun, 29 May 2016 14:11:12 +0200 Subject: [PATCH 09/25] [thermo] adding unit tests for util._get_umbrella_bias_sequences() --- pyemma/thermo/tests/test_util.py | 33 ++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index 324dd0a93..65b8e8d7a 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/pyemma/thermo/tests/test_util.py @@ -248,6 +248,39 @@ def test_umbrella_sampling_parameters_1x1(self): [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) +class TestProtectedUmbrellaSamplingBiasSequence(unittest.TestCase): + + def _assert_bias_sequences(self, bias_sequences, ref_bias_sequences): + self.assertTrue(len(bias_sequences) == len(ref_bias_sequences)) + for bs, rbs in zip(bias_sequences, ref_bias_sequences): + np.testing.assert_array_equal(bs, rbs) + + def test_umbrella_sampling_bias_sequences_1x0(self): + trajs = [np.array([0.0, 0.5, 1.0])] + umbrella_centers = np.array([ + util._ensure_umbrella_center(0.0, 1), + util._ensure_umbrella_center(1.0, 1)], dtype=np.float64) + force_constants = np.array([ + util._ensure_force_constant(1.0, 1), + util._ensure_force_constant(2.0, 1)], dtype=np.float64) + self._assert_bias_sequences( + util._get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants), + [np.array([[0.0, 1.0], [0.125, 0.25], [0.5, 0.0]])]) + + def test_umbrella_sampling_bias_sequences_1x1(self): + trajs = [np.array([[0.0], [0.5], [1.0]])] + umbrella_centers = np.array([ + util._ensure_umbrella_center(0.0, 1), + util._ensure_umbrella_center(1.0, 1)], dtype=np.float64) + force_constants = np.array([ + util._ensure_force_constant(1.0, 1), + util._ensure_force_constant(2.0, 1)], dtype=np.float64) + self._assert_bias_sequences( + util._get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants), + [np.array([[0.0, 1.0], [0.125, 0.25], [0.5, 0.0]])]) + + +#bias_sequences = _get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants) From 8ff4e975417da6051bddb2b5eb00043fe8a1ef8b Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Sun, 29 May 2016 14:15:05 +0200 Subject: [PATCH 10/25] [thermo] making util._get_umbrella_bias_sequences() robust --- pyemma/thermo/util/util.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index bfa3c7399..fc604b825 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -174,6 +174,8 @@ def _get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants): from thermotools.util import get_umbrella_bias as _get_umbrella_bias bias_sequences = [] for traj in trajs: + if traj.ndim == 1: + traj = traj.reshape((-1, 1)) bias_sequences.append( _get_umbrella_bias(traj, umbrella_centers, force_constants)) return bias_sequences From d224e9d9b72235231bb4131292432d2656b4cfe6 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 6 Jun 2016 14:13:59 +0200 Subject: [PATCH 11/25] [thermo] adding negative tests for protected umbrella sampling helper functions --- pyemma/thermo/tests/test_util.py | 75 +++++++++++++++++++++++++++++++- 1 file changed, 73 insertions(+), 2 deletions(-) diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index 65b8e8d7a..dbc92f3f0 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/pyemma/thermo/tests/test_util.py @@ -94,6 +94,20 @@ def test_ensure_umbrella_center_from_ndarray(self): self._assert_us_center(us_center, 4) np.testing.assert_array_equal(us_center, np.array([1.0, 2.0, 3.0, 4.0], dtype=np.float64)) + def test_ensure_umbrella_center_catches_unmatching_dimension(self): + with self.assertRaises(ValueError): + util._ensure_umbrella_center([1.0, 1.0], 1) + with self.assertRaises(ValueError): + util._ensure_umbrella_center([1.0, 1.0, 1.0], 2) + with self.assertRaises(ValueError): + util._ensure_umbrella_center([1.0, 1.0], 3) + with self.assertRaises(ValueError): + util._ensure_umbrella_center([[1.0], [1.0]], 1) + with self.assertRaises(ValueError): + util._ensure_umbrella_center([[1.0], [1.0]], 3) + with self.assertRaises(ValueError): + util._ensure_umbrella_center([[1.0, 1.0], [1.0]], 3) + class TestProtectedUmbrellaSamplingForceMatrices(unittest.TestCase): @@ -175,6 +189,24 @@ def test_ensure_umbrella_force_matrix_from_ndarray(self): np.testing.assert_array_equal( us_force_matrix, np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float64)) + def test_ensure_umbrella_force_matrix_catches_unmatching_dimension(self): + with self.assertRaises(ValueError): + util._ensure_force_constant([1.0, 1.0], 1) + with self.assertRaises(ValueError): + util._ensure_force_constant([1.0, 1.0, 1.0], 2) + with self.assertRaises(ValueError): + util._ensure_force_constant([1.0, 1.0], 3) + with self.assertRaises(ValueError): + util._ensure_force_constant([[1.0], [1.0]], 1) + with self.assertRaises(ValueError): + util._ensure_force_constant([[1.0], [1.0]], 3) + with self.assertRaises(ValueError): + util._ensure_force_constant([[1.0, 1.0], 1.0], 3) + with self.assertRaises(ValueError): + util._ensure_force_constant([1.0, [1.0, 1.0]], 3) + with self.assertRaises(ValueError): + util._ensure_force_constant([[1.0, 1.0], [1.0, 1.0]], 3) + class TestProtectedUmbrellaSamplingParameters(unittest.TestCase): @@ -279,8 +311,47 @@ def test_umbrella_sampling_bias_sequences_1x1(self): util._get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants), [np.array([[0.0, 1.0], [0.125, 0.25], [0.5, 0.0]])]) - -#bias_sequences = _get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants) + def test_umbrella_sampling_bias_sequences_catches_unmatching_dimension(self): + # wrong centers + constants + with self.assertRaises(TypeError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + np.array([[1.0, 1.0]]), [[[1.0, 0.0], [1.0, 0.0]]]) + with self.assertRaises(TypeError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + [[1.0, 1.0]], np.array([[[1.0, 0.0], [1.0, 0.0]]])) + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + np.array([1.0, 1.0]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + np.array([[1.0, 1.0]]), np.array([[1.0, 0.0], [1.0, 0.0]])) + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + np.array([[[1.0, 1.0]]]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + np.array([[1.0, 1.0]]), np.array([[[[1.0, 0.0], [1.0, 0.0]]]])) + # conflicting centers + constants + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + np.array([[1.0, 1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + np.array([[1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences([np.array([[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]])], + np.array([[1.0, 1.0], [2.0, 2.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) + # traj does not match valid centers + constants + with self.assertRaises(TypeError): + util._get_umbrella_bias_sequences([[[0.0, 0.0], [0.5, 0.1], [1.0, 0.2]]], + np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences([np.array([0.0, 0.5, 1.0])], + np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) + with self.assertRaises(ValueError): + util._get_umbrella_bias_sequences( + [np.array([[0.0, 1.0, 2.0], [0.5, 1.0, 2.0], [1.0, 1.0, 2.0]])], + np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) From 2203da5c0711069e23814a161e83157f7b841770 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 6 Jun 2016 14:19:03 +0200 Subject: [PATCH 12/25] [thermo] refactoring protected umbrella sampling helper functions to satisfy new tests --- pyemma/thermo/util/util.py | 54 +++++++++++++++++++++++++++++++------- 1 file changed, 44 insertions(+), 10 deletions(-) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index fc604b825..255f1fad9 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -87,28 +87,38 @@ def get_averaged_bias_matrix(bias_sequences, dtrajs, nstates=None): # ================================================================================================== def _ensure_umbrella_center(candidate, dimension): - candidate = _np.asarray(candidate).astype(_np.float64).reshape((-1,)) + try: + candidate = _np.asarray(candidate).astype(_np.float64).reshape((-1,)) + except ValueError: + raise ValueError("Umbrella center " + str(candidate) + " cannot be cast as numpy.ndarray") if candidate.shape[0] == 1 and dimension > 1: return candidate[0] * _np.ones(shape=(dimension,), dtype=_np.float64) + elif candidate.shape[0] != dimension: + raise ValueError("Unmatching dimensions: umbrella center " + str(candidate) + \ + " is not compatible with dimension %d" % dimension) return candidate def _ensure_force_constant(candidate, dimension): - candidate = _np.asarray(candidate).astype(_np.float64) + try: + candidate = _np.asarray(candidate).astype(_np.float64) + except ValueError: + raise ValueError("Force constant " + str(candidate) + " cannot be cast as numpy.ndarray") if candidate.ndim == 0: candidate = candidate * _np.ones(shape=(dimension,), dtype=_np.float64) + if candidate.shape[0] != dimension: + raise ValueError("Force constant " + str(candidate) + \ + " has shape[0]=%d instead of %d" % (candidate.shape[0], dimension)) if candidate.ndim == 1: - assert candidate.shape[0] == dimension, \ - str(candidate) + " has shape[0] %d" % candidate.shape[0] _candidate = _np.zeros(shape=(dimension, dimension), dtype=_np.float64) for i, x in enumerate(candidate): _candidate[i, i] = x candidate = _candidate + elif candidate.ndim == 2: + if candidate.shape[0] != dimension: + raise ValueError("Force constant " + str(candidate) + \ + " has shape[1]=%d instead of %d" % (candidate.shape[1], dimension)) else: - assert candidate.ndim == 2, str(candidate) + " has wrong ndim %d" % candidate.ndim - assert candidate.shape[0] == dimension, \ - str(candidate) + " has shape[0] %d" % candidate.shape[0] - assert candidate.shape[1] == dimension, \ - str(candidate) + " has shape[1] %d" % candidate.shape[1] + raise ValueError("Force constant " + str(candidate) + " must be a 2d numpy.ndarray") return candidate def _get_umbrella_sampling_parameters( @@ -173,9 +183,33 @@ def _get_umbrella_sampling_parameters( def _get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants): from thermotools.util import get_umbrella_bias as _get_umbrella_bias bias_sequences = [] - for traj in trajs: + if not isinstance(umbrella_centers, _np.ndarray): + raise TypeError("umbrella_centers is not a numpy.ndarray: " + str(type(umbrella_centers))) + if not isinstance(force_constants, _np.ndarray): + raise TypeError("force_constants is not a numpy.ndarray: " + str(type(force_constants))) + if umbrella_centers.ndim != 2: + raise ValueError("umbrella_centers is not a 2d numpy.ndarray: " + \ + str(umbrella_centers.shape)) + if force_constants.ndim != 3: + raise ValueError("force_constants is not a 3d numpy.ndarray: " + \ + str(force_constants.shape)) + if umbrella_centers.shape[0] != force_constants.shape[0]: + raise ValueError("Unmatching number of umbrella centers and force constants: %d != %d" % ( + umbrella_centers.shape[0], force_constants.shape[0])) + dimension = umbrella_centers.shape[1] + if force_constants.shape[1] != dimension or force_constants.shape[2] != dimension: + raise ValueError("Dimension of force_constants does not match dimension of " + \ + "umbrella_centers: %d != %d,%d" % (dimension, + force_constants.shape[1], force_constants.shape[2])) + for i, traj in enumerate(trajs): + if not isinstance(traj, _np.ndarray): + raise TypeError("Trajectory %d is not a numpy.ndarray: " % i + str( + type(traj))) if traj.ndim == 1: traj = traj.reshape((-1, 1)) + if traj.shape[1] != dimension: + raise ValueError("Trajectory %d has unmatching dimension: %d!=%d" % ( + i, traj.shape[1], dimension)) bias_sequences.append( _get_umbrella_bias(traj, umbrella_centers, force_constants)) return bias_sequences From 3ca7f04234a4d7617088f3c39a25fe9307aaddad Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 6 Jun 2016 14:34:12 +0200 Subject: [PATCH 13/25] [thermo] replacing asserts is API and protected umbrella sampling helpers with if+raise constructs --- pyemma/thermo/api.py | 38 ++++++++++++++++++++++++++++---------- pyemma/thermo/util/util.py | 6 ++++-- 2 files changed, 32 insertions(+), 12 deletions(-) diff --git a/pyemma/thermo/api.py b/pyemma/thermo/api.py index 31bed65e3..eb19c7446 100644 --- a/pyemma/thermo/api.py +++ b/pyemma/thermo/api.py @@ -143,7 +143,8 @@ def estimate_umbrella_sampling( array([ 0.63..., 1.60..., 1.31...]) """ - assert estimator in ['wham', 'dtram', 'tram'], "unsupported estimator: %s" % estimator + if estimator not in ['wham', 'dtram', 'tram']: + raise ValueError("unsupported estimator: %s" % estimator) from .util import get_umbrella_sampling_data as _get_umbrella_sampling_data ttrajs, btrajs, umbrella_centers, force_constants, unbiased_state = _get_umbrella_sampling_data( us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT) @@ -284,7 +285,8 @@ def estimate_multi_temperature( can handle temperature changes as well. """ - assert estimator in ['wham', 'dtram', 'tram'], "unsupported estimator: %s" % estimator + if estimator not in ['wham', 'dtram', 'tram']: + ValueError("unsupported estimator: %s" % estimator) from .util import get_multi_temperature_data as _get_multi_temperature_data ttrajs, btrajs, temperatures, unbiased_state = _get_multi_temperature_data( energy_trajs, temp_trajs, energy_unit, temp_unit, @@ -470,11 +472,19 @@ def tram( # prepare trajectories ttrajs = _types.ensure_dtraj_list(ttrajs) dtrajs = _types.ensure_dtraj_list(dtrajs) - assert len(ttrajs) == len(dtrajs) - assert len(ttrajs) == len(bias) + if len(ttrajs) != len(dtrajs): + raise ValueError("Unmatching number of dtraj/ttraj elements: %d!=%d" % ( + len(dtrajs), len(ttrajs))) + if len(ttrajs) != len(bias): + raise ValueError("Unmatching number of ttraj/bias elements: %d!=%d" % ( + len(ttrajs), len(bias))) for ttraj, dtraj, btraj in zip(ttrajs, dtrajs, bias): - assert len(ttraj) == len(dtraj) - assert len(ttraj) == btraj.shape[0] + if len(ttraj) != len(dtraj): + raise ValueError("Unmatching number of data points in ttraj/dtraj: %d!=%d" % ( + len(ttraj), len(dtraj))) + if len(ttraj) != btraj.shape[0]: + raise ValueError("Unmatching number of data points in ttraj/bias trajectory: %d!=%d" % ( + len(ttraj), len(btraj))) # check lag time(s) lags = _np.asarray(lag, dtype=_np.intc).reshape((-1,)).tolist() # build TRAM and run estimation @@ -630,9 +640,13 @@ def dtram( # prepare trajectories ttrajs = _types.ensure_dtraj_list(ttrajs) dtrajs = _types.ensure_dtraj_list(dtrajs) - assert len(ttrajs) == len(dtrajs) + if len(ttrajs) != len(dtrajs): + raise ValueError("Unmatching number of dtraj/ttraj elements: %d!=%d" % ( + len(dtrajs), len(ttrajs)) ) for ttraj, dtraj in zip(ttrajs, dtrajs): - assert len(ttraj) == len(dtraj) + if len(ttraj) != len(dtraj): + raise ValueError("Unmatching number of data points in ttraj/dtraj: %d!=%d" % ( + len(ttraj), len(dtraj))) # check lag time(s) lags = _np.asarray(lag, dtype=_np.intc).reshape((-1,)).tolist() # build DTRAM and run estimation @@ -760,9 +774,13 @@ def wham( # check trajectories ttrajs = _types.ensure_dtraj_list(ttrajs) dtrajs = _types.ensure_dtraj_list(dtrajs) - assert len(ttrajs) == len(dtrajs) + if len(ttrajs) != len(dtrajs): + raise ValueError("Unmatching number of dtraj/ttraj elements: %d!=%d" % ( + len(dtrajs), len(ttrajs)) ) for ttraj, dtraj in zip(ttrajs, dtrajs): - assert len(ttrajs) == len(dtrajs) + if len(ttraj) != len(dtraj): + raise ValueError("Unmatching number of data points in ttraj/dtraj: %d!=%d" % ( + len(ttraj), len(dtraj))) # build WHAM from pyemma.thermo import WHAM wham_estimator = WHAM( diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index 255f1fad9..e6265da22 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -175,8 +175,10 @@ def _get_umbrella_sampling_parameters( umbrella_centers = _np.array(umbrella_centers, dtype=_np.float64) force_constants = _np.array(force_constants, dtype=_np.float64) if kT is not None: - assert isinstance(kT, (int, float)), "kT has wrong type:" + str(type(kT)) - assert kT > 0.0, "non-positive kT: %f" % kT + if not isinstance(kT, (int, float)): + raise ValueError("kT has wrong type:" + str(type(kT))) + if kT <= 0.0: + raise ValueError("non-positive kT: %f" % kT) force_constants /= kT return ttrajs, umbrella_centers, force_constants, unbiased_state From a0379a160d9d7ef8eaf5dc2d6302599e16b4a7c7 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 6 Jun 2016 17:04:59 +0200 Subject: [PATCH 14/25] [thermo] adding tests for the protected umbrella sampling helper functions for the case that only one trajectory is supplied --- pyemma/thermo/tests/test_util.py | 51 +++++++++++++++++++++++--------- 1 file changed, 37 insertions(+), 14 deletions(-) diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index dbc92f3f0..c48bf11db 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/pyemma/thermo/tests/test_util.py @@ -243,6 +243,12 @@ def test_umbrella_sampling_parameters_1x0(self): [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) # with kT parameter + with self.assertRaises(ValueError): + util._get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=0.0) + with self.assertRaises(ValueError): + util._get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT='kT') ttrajs, umbrella_centers, force_constants, unbiased_state = \ util._get_umbrella_sampling_parameters( us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) @@ -250,6 +256,18 @@ def test_umbrella_sampling_parameters_1x0(self): ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) + # single trajectory outside list + ref_umbrella_centers = 1.0 + ref_force_constants = 1.0 + us_trajs = np.array([0.0, 0.1, 0.2]) + md_trajs = np.array([0.0, 0.5, 1.0]) + ttrajs, umbrella_centers, force_constants, unbiased_state = \ + util._get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) + self._assert_parameters( + ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, + [np.array([0, 0, 0]), np.array([1, 1, 1])], + [ref_umbrella_centers], [ref_force_constants], 1) def test_umbrella_sampling_parameters_1x1(self): ref_umbrella_centers = [0.0, 1.0] @@ -272,6 +290,12 @@ def test_umbrella_sampling_parameters_1x1(self): [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) # with kT parameter + with self.assertRaises(ValueError): + util._get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=0.0) + with self.assertRaises(ValueError): + util._get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT='kT') ttrajs, umbrella_centers, force_constants, unbiased_state = \ util._get_umbrella_sampling_parameters( us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) @@ -279,6 +303,19 @@ def test_umbrella_sampling_parameters_1x1(self): ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) + # single trajectory outside list + ref_umbrella_centers = 1.0 + ref_force_constants = 1.0 + us_trajs = np.array([[0.0], [0.1], [0.2]]) + md_trajs = np.array([[0.0], [0.5], [1.0]]) + ttrajs, umbrella_centers, force_constants, unbiased_state = \ + util._get_umbrella_sampling_parameters( + us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) + self._assert_parameters( + ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, + [np.array([0, 0, 0]), np.array([1, 1, 1])], + [ref_umbrella_centers], [ref_force_constants], 1) + class TestProtectedUmbrellaSamplingBiasSequence(unittest.TestCase): @@ -352,17 +389,3 @@ def test_umbrella_sampling_bias_sequences_catches_unmatching_dimension(self): util._get_umbrella_bias_sequences( [np.array([[0.0, 1.0, 2.0], [0.5, 1.0, 2.0], [1.0, 1.0, 2.0]])], np.array([[1.0, 1.0]]), np.array([[[1.0, 0.0], [1.0, 0.0]]])) - - - - - - - - - - - - - - From be3f399e10863b24294b65ab4a6692cb0323f7af Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 6 Jun 2016 17:05:34 +0200 Subject: [PATCH 15/25] [thermo] amending the protected umbrella sampling helper functions for the case that only one trajectory is supplied --- pyemma/thermo/util/util.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index e6265da22..90e3d015d 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -129,6 +129,19 @@ def _get_umbrella_sampling_parameters( nthermo = 0 unbiased_state = None dimension = None + if not isinstance(us_trajs, (list, tuple)): + us_trajs = [us_trajs] + if not isinstance(us_centers, (list, tuple)): + us_centers = [us_centers] + if not isinstance(us_force_constants, (list, tuple)): + us_force_constants = [us_force_constants] + if len(us_trajs) != len(us_centers): + raise ValueError("Unmatching number of umbrella sampling trajectories and centers: %d!=%d" \ + % (len(us_trajs), len(us_centers))) + if len(us_trajs) != len(us_force_constants): + raise ValueError( + "Unmatching number of umbrella sampling trajectories and force constants: %d!=%d" \ + % (len(us_trajs), len(us_force_constants))) for i, traj in enumerate(us_trajs): state = None try: @@ -159,6 +172,8 @@ def _get_umbrella_sampling_parameters( else: ttrajs.append(state * _np.ones(shape=(traj.shape[0],), dtype=_np.intc)) if md_trajs is not None: + if not isinstance(md_trajs, (list, tuple)): + md_trajs = [md_trajs] this_center = umbrella_centers[-1] * 0.0 this_force_constant = force_constants[-1] * 0.0 for j in range(nthermo): From 5f6154674aa9ec5060fecab3e3d5e587572a655e Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 14:52:15 +0200 Subject: [PATCH 16/25] [thermo] reverting support for single trajectories in estimate_umbrella_sampling() --- pyemma/thermo/tests/test_util.py | 24 ------------------------ pyemma/thermo/util/util.py | 8 +++++--- 2 files changed, 5 insertions(+), 27 deletions(-) diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index c48bf11db..70551d34c 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/pyemma/thermo/tests/test_util.py @@ -256,18 +256,6 @@ def test_umbrella_sampling_parameters_1x0(self): ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) - # single trajectory outside list - ref_umbrella_centers = 1.0 - ref_force_constants = 1.0 - us_trajs = np.array([0.0, 0.1, 0.2]) - md_trajs = np.array([0.0, 0.5, 1.0]) - ttrajs, umbrella_centers, force_constants, unbiased_state = \ - util._get_umbrella_sampling_parameters( - us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) - self._assert_parameters( - ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, - [np.array([0, 0, 0]), np.array([1, 1, 1])], - [ref_umbrella_centers], [ref_force_constants], 1) def test_umbrella_sampling_parameters_1x1(self): ref_umbrella_centers = [0.0, 1.0] @@ -303,18 +291,6 @@ def test_umbrella_sampling_parameters_1x1(self): ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) - # single trajectory outside list - ref_umbrella_centers = 1.0 - ref_force_constants = 1.0 - us_trajs = np.array([[0.0], [0.1], [0.2]]) - md_trajs = np.array([[0.0], [0.5], [1.0]]) - ttrajs, umbrella_centers, force_constants, unbiased_state = \ - util._get_umbrella_sampling_parameters( - us_trajs, ref_umbrella_centers, ref_force_constants, md_trajs=md_trajs, kT=2.0) - self._assert_parameters( - ttrajs, umbrella_centers, force_constants * 2.0, unbiased_state, - [np.array([0, 0, 0]), np.array([1, 1, 1])], - [ref_umbrella_centers], [ref_force_constants], 1) class TestProtectedUmbrellaSamplingBiasSequence(unittest.TestCase): diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index 90e3d015d..d0a0e3c17 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -130,11 +130,13 @@ def _get_umbrella_sampling_parameters( unbiased_state = None dimension = None if not isinstance(us_trajs, (list, tuple)): - us_trajs = [us_trajs] + raise ValueError("The parameter us_trajs must be a list of numpy.ndarray objects") if not isinstance(us_centers, (list, tuple)): - us_centers = [us_centers] + raise ValueError( + "The parameter us_centers must be a list of floats or numpy.ndarray objects") if not isinstance(us_force_constants, (list, tuple)): - us_force_constants = [us_force_constants] + raise ValueError( + "The parameter us_force_constants must be a list of floats or numpy.ndarray objects") if len(us_trajs) != len(us_centers): raise ValueError("Unmatching number of umbrella sampling trajectories and centers: %d!=%d" \ % (len(us_trajs), len(us_centers))) From 8da1247caebc21c8f916c8dc4a86ea9452b8311c Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 14:53:36 +0200 Subject: [PATCH 17/25] [thermo] making docstring for estimate_umbrella_sampling() explicit on how to pass trajecttories, us_centers, and us_force_constants --- pyemma/thermo/api.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pyemma/thermo/api.py b/pyemma/thermo/api.py index eb19c7446..1f2a94122 100644 --- a/pyemma/thermo/api.py +++ b/pyemma/thermo/api.py @@ -57,10 +57,10 @@ def estimate_umbrella_sampling( us_dtrajs : list of N int arrays, each of shape (T_i,) The integers are indexes in 0,...,n-1 enumerating the n discrete states or the bins the umbrella sampling trajectory is in at any time. - us_centers : array-like of size N + us_centers : list of N floats or d-dimensional arrays of floats List or array of N center positions. Each position must be a d-dimensional vector. For 1d umbrella sampling, one can simply pass a list of centers, e.g. [-5.0, -4.0, -3.0, ... ]. - us_force_constants : float or array-like of float + us_force_constants : list of N floats or d- or dxd-dimensional arrays of floats The force constants used in the umbrellas, unit-less (e.g. kT per squared length unit). For multidimensional umbrella sampling, the force matrix must be used. md_trajs : list of M arrays, each of shape (T_i, d), optional, default=None @@ -148,6 +148,9 @@ def estimate_umbrella_sampling( from .util import get_umbrella_sampling_data as _get_umbrella_sampling_data ttrajs, btrajs, umbrella_centers, force_constants, unbiased_state = _get_umbrella_sampling_data( us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT) + for traj, dtraj in zip(us_trajs, us_dtrajs): + if traj.shape[0] != dtraj.shape[0]: + raise ValueError("unmatching trajectory lengths") if md_dtrajs is None: md_dtrajs = [] estimator_obj = None From cff81bcfdbb6dbc2e78f9ca7a5d860df5ea26fd5 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 15:13:39 +0200 Subject: [PATCH 18/25] [thermo] adding unit test to check for exceptions for trajectory mismatches in estimate_umbrella_sampling() --- pyemma/thermo/tests/test_api.py | 78 +++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 pyemma/thermo/tests/test_api.py diff --git a/pyemma/thermo/tests/test_api.py b/pyemma/thermo/tests/test_api.py new file mode 100644 index 000000000..ba224f1a1 --- /dev/null +++ b/pyemma/thermo/tests/test_api.py @@ -0,0 +1,78 @@ +# This file is part of PyEMMA. +# +# Copyright (c) 2016 Computational Molecular Biology Group, Freie Universitaet Berlin (GER) +# +# PyEMMA is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + +import unittest +import numpy as np +from pyemma.thermo import estimate_umbrella_sampling +from pyemma.coordinates import cluster_regspace + +# ================================================================================================== +# tests for the umbrella sampling API +# ================================================================================================== + +class TestProtectedUmbrellaSamplingCenters(unittest.TestCase): + + def test_exceptions(self): + us_centers = [1.1, 1.3] + us_force_constants = [1.0, 1.0] + us_trajs = [ + np.array([1.0, 1.1, 1.2, 1.1, 1.0, 1.1]), + np.array([1.3, 1.2, 1.3, 1.4, 1.4, 1.3])] + md_trajs = [ + np.array([0.9, 1.0, 1.1, 1.2, 1.3, 1.4]), + np.array([1.5, 1.4, 1.3, 1.4, 1.4, 1.5])] + cluster = cluster_regspace(data=us_trajs+md_trajs, max_centers=10, dmin=0.15) + us_dtrajs = cluster.dtrajs[:2] + md_dtrajs = cluster.dtrajs[2:] + # unmatching number of us trajectories / us parameters + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs[:-1], us_dtrajs, us_centers, us_force_constants) + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs[:-1], us_centers, us_force_constants) + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs, us_centers[:-1], us_force_constants) + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs, us_centers, us_force_constants[:-1]) + # unmatching number of md trajectories + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs, us_centers, us_force_constants, + md_trajs=md_trajs[:-1], md_dtrajs=md_dtrajs) + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs, us_centers, us_force_constants, + md_trajs=md_trajs, md_dtrajs=md_dtrajs[:-1]) + # unmatchig trajectory lengths + us_trajs_x = [ + np.array([1.0, 1.1, 1.2, 1.1, 1.0]), + np.array([1.3, 1.2, 1.3, 1.4, 1.4])] + md_trajs_x = [ + np.array([0.9, 1.0, 1.1, 1.2, 1.3]), + np.array([1.5, 1.4, 1.3, 1.4, 1.4])] + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs_x, us_dtrajs, us_centers, us_force_constants) + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs, us_centers, us_force_constants, + md_trajs=md_trajs_x, md_dtrajs=md_dtrajs) + + From 499193f73abb4d26872629a4a322d2c92157ab91 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 15:20:16 +0200 Subject: [PATCH 19/25] [thermo] checking for trajectory inconsistencies in estimate_umbrella_sampling() --- pyemma/thermo/api.py | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/pyemma/thermo/api.py b/pyemma/thermo/api.py index 1f2a94122..9b4f54033 100644 --- a/pyemma/thermo/api.py +++ b/pyemma/thermo/api.py @@ -148,11 +148,31 @@ def estimate_umbrella_sampling( from .util import get_umbrella_sampling_data as _get_umbrella_sampling_data ttrajs, btrajs, umbrella_centers, force_constants, unbiased_state = _get_umbrella_sampling_data( us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT) + if len(us_trajs) != len(us_dtrajs): + raise ValueError( + "Number of continuous and discrete umbrella sampling trajectories does not " + \ + "match: %d!=%d" % (len(us_trajs), len(us_dtrajs))) + i = 0 for traj, dtraj in zip(us_trajs, us_dtrajs): if traj.shape[0] != dtraj.shape[0]: - raise ValueError("unmatching trajectory lengths") + raise ValueError( + "Lengths of continuous and discrete umbrella sampling trajectories with " + \ + "index %d does not match: %d!=%d" % (i, len(us_trajs), len(us_dtrajs))) + i += 1 if md_dtrajs is None: md_dtrajs = [] + else: + if len(md_trajs) != len(md_dtrajs): + raise ValueError( + "Number of continuous and discrete unbiased trajectories does not " + \ + "match: %d!=%d" % (len(md_trajs), len(md_dtrajs))) + i = 0 + for traj, dtraj in zip(md_trajs, md_dtrajs): + if traj.shape[0] != dtraj.shape[0]: + raise ValueError( + "Lengths of continuous and discrete unbiased trajectories with " + \ + "index %d does not match: %d!=%d" % (i, len(md_trajs), len(md_dtrajs))) + i += 1 estimator_obj = None if estimator == 'wham': estimator_obj = wham( From b68fabcfcc31802860fccdc20bcfd46f14775008 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 16:42:36 +0200 Subject: [PATCH 20/25] [thermo] adding type check for md_trajs --- pyemma/thermo/util/util.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index d0a0e3c17..ec7143078 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -144,6 +144,9 @@ def _get_umbrella_sampling_parameters( raise ValueError( "Unmatching number of umbrella sampling trajectories and force constants: %d!=%d" \ % (len(us_trajs), len(us_force_constants))) + if md_trajs is not None: + if not isinstance(md_trajs, (list, tuple)): + raise ValueError("The parameter md_trajs must be a list of numpy.ndarray objects") for i, traj in enumerate(us_trajs): state = None try: From 090b80bc581d7a88d07cddcc475a4c0031d75caf Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 16:54:25 +0200 Subject: [PATCH 21/25] [thermo] adding exception tests for single trajectories and md_trajs/md_dtrajs None cases --- pyemma/thermo/tests/test_api.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/pyemma/thermo/tests/test_api.py b/pyemma/thermo/tests/test_api.py index ba224f1a1..d7ca7d4a5 100644 --- a/pyemma/thermo/tests/test_api.py +++ b/pyemma/thermo/tests/test_api.py @@ -74,5 +74,21 @@ def test_exceptions(self): estimate_umbrella_sampling( us_trajs, us_dtrajs, us_centers, us_force_constants, md_trajs=md_trajs_x, md_dtrajs=md_dtrajs) - + # unmatching md_trajs/md_dtrajs cases + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs, us_centers, us_force_constants, + md_trajs=None, md_dtrajs=md_dtrajs) + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs, us_centers, us_force_constants, + md_trajs=md_trajs, md_dtrajs=None) + # single trajectory cases + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs[0], us_dtrajs[0], us_centers[0], us_force_constants[0]) + with self.assertRaises(ValueError): + estimate_umbrella_sampling( + us_trajs, us_dtrajs, us_centers, us_force_constants, + md_trajs=md_trajs[0], md_dtrajs=md_dtrajs[0]) From 7b1359dbf3304df26d35226dcdcfb97adb079428 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 16:55:06 +0200 Subject: [PATCH 22/25] [thermo] moving type checks to estimate_umbrella_sampling() and minor refactoring --- pyemma/thermo/api.py | 32 +++++++++++++++++++++++++++++--- pyemma/thermo/util/util.py | 18 ------------------ 2 files changed, 29 insertions(+), 21 deletions(-) diff --git a/pyemma/thermo/api.py b/pyemma/thermo/api.py index 9b4f54033..ca0a1de2e 100644 --- a/pyemma/thermo/api.py +++ b/pyemma/thermo/api.py @@ -143,11 +143,25 @@ def estimate_umbrella_sampling( array([ 0.63..., 1.60..., 1.31...]) """ + from .util import get_umbrella_sampling_data as _get_umbrella_sampling_data + # sanity checks if estimator not in ['wham', 'dtram', 'tram']: raise ValueError("unsupported estimator: %s" % estimator) - from .util import get_umbrella_sampling_data as _get_umbrella_sampling_data - ttrajs, btrajs, umbrella_centers, force_constants, unbiased_state = _get_umbrella_sampling_data( - us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT) + if not isinstance(us_trajs, (list, tuple)): + raise ValueError("The parameter us_trajs must be a list of numpy.ndarray objects") + if not isinstance(us_centers, (list, tuple)): + raise ValueError( + "The parameter us_centers must be a list of floats or numpy.ndarray objects") + if not isinstance(us_force_constants, (list, tuple)): + raise ValueError( + "The parameter us_force_constants must be a list of floats or numpy.ndarray objects") + if len(us_trajs) != len(us_centers): + raise ValueError("Unmatching number of umbrella sampling trajectories and centers: %d!=%d" \ + % (len(us_trajs), len(us_centers))) + if len(us_trajs) != len(us_force_constants): + raise ValueError( + "Unmatching number of umbrella sampling trajectories and force constants: %d!=%d" \ + % (len(us_trajs), len(us_force_constants))) if len(us_trajs) != len(us_dtrajs): raise ValueError( "Number of continuous and discrete umbrella sampling trajectories does not " + \ @@ -159,9 +173,16 @@ def estimate_umbrella_sampling( "Lengths of continuous and discrete umbrella sampling trajectories with " + \ "index %d does not match: %d!=%d" % (i, len(us_trajs), len(us_dtrajs))) i += 1 + if md_trajs is not None: + if not isinstance(md_trajs, (list, tuple)): + raise ValueError("The parameter md_trajs must be a list of numpy.ndarray objects") + if md_dtrajs is None: + raise ValueError("You have provided md_trajs, but md_dtrajs is None") if md_dtrajs is None: md_dtrajs = [] else: + if md_trajs is None: + raise ValueError("You have provided md_dtrajs, but md_trajs is None") if len(md_trajs) != len(md_dtrajs): raise ValueError( "Number of continuous and discrete unbiased trajectories does not " + \ @@ -173,7 +194,11 @@ def estimate_umbrella_sampling( "Lengths of continuous and discrete unbiased trajectories with " + \ "index %d does not match: %d!=%d" % (i, len(md_trajs), len(md_dtrajs))) i += 1 + # data preparation + ttrajs, btrajs, umbrella_centers, force_constants, unbiased_state = _get_umbrella_sampling_data( + us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT) estimator_obj = None + # estimation if estimator == 'wham': estimator_obj = wham( ttrajs, us_dtrajs + md_dtrajs, @@ -200,6 +225,7 @@ def estimate_umbrella_sampling( maxiter=maxiter, maxerr=maxerr, save_convergence_info=save_convergence_info, dt_traj=dt_traj, init=init, init_maxiter=init_maxiter, init_maxerr=init_maxerr, **parsed_kwargs) + # adding thermodynamic state information and return results try: estimator_obj.umbrella_centers = umbrella_centers estimator_obj.force_constants = force_constants diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index ec7143078..db281ec57 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -129,24 +129,6 @@ def _get_umbrella_sampling_parameters( nthermo = 0 unbiased_state = None dimension = None - if not isinstance(us_trajs, (list, tuple)): - raise ValueError("The parameter us_trajs must be a list of numpy.ndarray objects") - if not isinstance(us_centers, (list, tuple)): - raise ValueError( - "The parameter us_centers must be a list of floats or numpy.ndarray objects") - if not isinstance(us_force_constants, (list, tuple)): - raise ValueError( - "The parameter us_force_constants must be a list of floats or numpy.ndarray objects") - if len(us_trajs) != len(us_centers): - raise ValueError("Unmatching number of umbrella sampling trajectories and centers: %d!=%d" \ - % (len(us_trajs), len(us_centers))) - if len(us_trajs) != len(us_force_constants): - raise ValueError( - "Unmatching number of umbrella sampling trajectories and force constants: %d!=%d" \ - % (len(us_trajs), len(us_force_constants))) - if md_trajs is not None: - if not isinstance(md_trajs, (list, tuple)): - raise ValueError("The parameter md_trajs must be a list of numpy.ndarray objects") for i, traj in enumerate(us_trajs): state = None try: From 3806b7d14e68e5f6020d2e191936d9e6cd2f875d Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 17:05:54 +0200 Subject: [PATCH 23/25] [thermo] adding test for unmatching trajectory dimensions in util._get_umbrella_sampling_parameters() --- pyemma/thermo/tests/test_util.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/pyemma/thermo/tests/test_util.py b/pyemma/thermo/tests/test_util.py index 70551d34c..105e4306c 100644 --- a/pyemma/thermo/tests/test_util.py +++ b/pyemma/thermo/tests/test_util.py @@ -292,6 +292,15 @@ def test_umbrella_sampling_parameters_1x1(self): [np.array([0, 0, 0]), np.array([1, 1, 1]), np.array([2, 2, 2])], ref_umbrella_centers, ref_force_constants, 2) + def test_umbrella_sampling_parameters_unmatching_dimensions(self): + ref_umbrella_centers = [0.0, 1.0] + ref_force_constants = [1.0, 1.0] + us_trajs_x = [ + np.array([[0.0], [0.1], [0.2]]), np.array([[0.9, 0.0], [1.0, 0.0], [1.1, 0.0]])] + with self.assertRaises(ValueError): + util._get_umbrella_sampling_parameters( + us_trajs_x, ref_umbrella_centers, ref_force_constants) + class TestProtectedUmbrellaSamplingBiasSequence(unittest.TestCase): From 2d798bb4710264cd8847fa6e8725f4e9eaa6a35d Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 17:06:46 +0200 Subject: [PATCH 24/25] [thermo] replacing an assert by if, raise ValueError --- pyemma/thermo/util/util.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pyemma/thermo/util/util.py b/pyemma/thermo/util/util.py index db281ec57..cb5be54c7 100644 --- a/pyemma/thermo/util/util.py +++ b/pyemma/thermo/util/util.py @@ -138,8 +138,9 @@ def _get_umbrella_sampling_parameters( if dimension is None: dimension = _dimension else: - assert dimension == _dimension, "trajectory %i has unmatching dimension %d!=%d" % ( - i, _dimension, dimension) + if dimension != _dimension: + raise ValueError( + "Trajectory %i has unmatching dimension %d!=%d" % (i, _dimension, dimension)) this_center = _ensure_umbrella_center( us_centers[i], dimension) this_force_constant = _ensure_force_constant( @@ -207,8 +208,7 @@ def _get_umbrella_bias_sequences(trajs, umbrella_centers, force_constants): force_constants.shape[1], force_constants.shape[2])) for i, traj in enumerate(trajs): if not isinstance(traj, _np.ndarray): - raise TypeError("Trajectory %d is not a numpy.ndarray: " % i + str( - type(traj))) + raise TypeError("Trajectory %d is not a numpy.ndarray: " % i + str(type(traj))) if traj.ndim == 1: traj = traj.reshape((-1, 1)) if traj.shape[1] != dimension: From c169834d39fc6c228f70380cce69104127dee323 Mon Sep 17 00:00:00 2001 From: Christoph Wehmeyer Date: Mon, 13 Jun 2016 17:17:16 +0200 Subject: [PATCH 25/25] [doc] added changes to CHANGELOG --- doc/source/CHANGELOG.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/source/CHANGELOG.rst b/doc/source/CHANGELOG.rst index 14df10980..94713ded0 100644 --- a/doc/source/CHANGELOG.rst +++ b/doc/source/CHANGELOG.rst @@ -10,6 +10,7 @@ Changelog **Fixes**: - clustering: fixed KMeans with minRMSD metric. #814 +- thermo: made estimate_umbrella_sampling more robust w.r.t. input and fixed doumentation. #812 #827 2.2 (5-17-16) -------------