From aad7c4122c4367bcab2162d2e174088b4c46413b Mon Sep 17 00:00:00 2001 From: "Justin R. Porter" Date: Thu, 20 Dec 2018 09:54:23 -0600 Subject: [PATCH 1/2] wip rmsd in enspara --- enspara/geometry/libdist.pyx | 57 ++++++++++++++++++++++++++++++++++-- setup.py | 14 +++++++-- 2 files changed, 67 insertions(+), 4 deletions(-) diff --git a/enspara/geometry/libdist.pyx b/enspara/geometry/libdist.pyx index 9a8042cbc..a8acfbf41 100644 --- a/enspara/geometry/libdist.pyx +++ b/enspara/geometry/libdist.pyx @@ -29,18 +29,26 @@ cdef extern from "math.h" nogil: double fabs(double x) float fabs(float x) +cdef extern from "theobald_rmsd.h": + cdef extern float msd_atom_major( + int nrealatoms, int npaddedatoms, float* a, + float* b, float G_a, float G_b, int computeRot, float rot[9]) nogil + + def _check_is_2d(X): if len(X.shape) != 2: raise exception.DataInvalid( "Data array dimension must be two, got shape %s." % str(X.shape)) + def _check_is_1d(x): if len(x.shape) != 1: raise exception.DataInvalid( "Target point dimension must be one, got shape %s." % str(x.shape)) + def _prepare_for_2d_to_1d_distance(X, y, out): _check_is_2d(X) @@ -80,8 +88,10 @@ def _hamming(np.ndarray[INTEGRAL_TYPE_T, ndim=2] X, cdef long n_samples = len(out) cdef long n_features = len(y) - assert len(out) == X.shape[0], "Size of output array didn't match number of observations in X" - assert n_features == X.shape[1], "Number of features between X and y didn't match." + assert len(out) == X.shape[0], \ + "Size of output array didn't match number of observations in X" + assert n_features == X.shape[1], \ + "Number of features between X and y didn't match." cdef long i, j = 0 @@ -145,6 +155,48 @@ def _euclidean(np.ndarray[FLOAT_TYPE_T, ndim=2] X, return out.reshape(-1, 1) +cdef _rmsd(XA, XB): + cdef long i, j + cdef float[:, :, ::1] XA_xyz = XA.xyz + cdef float[:, :, ::1] XB_xyz = XB.xyz + cdef long XA_length = XA_xyz.shape[0] + cdef long XB_length = XB_xyz.shape[0] + cdef int n_atoms = XA_xyz.shape[1] + cdef float[::1] XA_trace, XB_trace + cdef double[:, ::1] out + + if XA._rmsd_traces is None: + XA.center_coordinates() + if XB._rmsd_traces is None: + XB.center_coordinates() + + if XA_xyz.shape[1] != XB_xyz.shape[1]: + raise ValueError('XA and XB must have the same number of atoms') + + XA_trace = XA._rmsd_traces + XB_trace = XB._rmsd_traces + + out = np.zeros((XA_length, XB_length), dtype=np.double) + + for i in range(XA_length): + for j in range(XB_length): + rmsd = sqrt(msd_atom_major(n_atoms, n_atoms, &XA_xyz[i,0,0], + &XB_xyz[j,0,0], XA_trace[i], XB_trace[j], 0, NULL)) + out[i,j] = rmsd + + return np.array(out, copy=False) + + +def rmsd(X, y, out=None): + + # out = _prepare_for_2d_to_1d_distance(XA, XB, out) + # _rmsd(X, y, out) + _rmsd(XA, AB, out) + + return out + + + def euclidean(X, y, out=None): """Compute the euclidean distance between a point, `y`, and a group of points `X`. Uses thread-parallelism with OpenMP. @@ -163,6 +215,7 @@ def euclidean(X, y, out=None): _euclidean(X, y, out) return out + def manhattan(X, y, out=None): """Compute the Manhattan distance between a point `y` and a group of points `X`. Thread-parallized using OpenMP. diff --git a/setup.py b/setup.py index 9aa51bc05..18d4d1631 100755 --- a/setup.py +++ b/setup.py @@ -37,6 +37,12 @@ 'or see http://docs.scipy.org/doc/numpy/user/install.html and' 'http://cython.org/ for more information.'])) +try: + import mdtraj + mdtraj_capi = mdtraj.capi() +except (ImportError, AttributeError): + print('MDTraj>=1.1 is required') + sys.exit(1) # this probably won't work for everyone. Works for me, though! # they'll need gcc 7 installed. Unfortunately, I don't have any idea how @@ -62,10 +68,14 @@ extra_compile_args=extra_compile_args, extra_link_args=extra_link_args, ), Extension( - "enspara.geometry.libdist", - ["enspara/geometry/libdist.pyx"], + name="enspara.geometry.libdist", + sources=["enspara/geometry/libdist.pyx"], + # msvc needs to be told "libtheobald", gcc wants just "theobald" + libraries=['theobald'], extra_compile_args=extra_compile_args, extra_link_args=extra_link_args, + include_dirs=[mdtraj_capi['include_dir'], np.get_include()], + library_dirs=[mdtraj_capi['lib_dir']], ), Extension( "enspara.msm.libmsm", ["enspara/msm/libmsm.pyx"], From 1387ed7e0fcc8c615d220d05b2415cb151f679a2 Mon Sep 17 00:00:00 2001 From: "Justin R. Porter" Date: Thu, 20 Dec 2018 13:41:58 -0600 Subject: [PATCH 2/2] fix setup.py's adminition about installing cython. --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 18d4d1631..6ddb75f20 100755 --- a/setup.py +++ b/setup.py @@ -30,7 +30,7 @@ except ImportError: sys.stderr.write('-' * 80) sys.stderr.write('\n'.join([ - 'Error: building mdtraj requires numpy and cython>=0.19', + 'Error: building enspara requires numpy and cython>=0.19', 'Try running the command ``pip install numpy cython`` or' '``conda install numpy cython``.',