1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
|
// $Id$
//
// Copyright (C) 2004-2008 Greg Landrum and Rational Discovery LLC
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#define PY_ARRAY_UNIQUE_SYMBOL rdalignment_array_API
#include <RDBoost/python.h>
#include <RDBoost/boost_numpy.h>
#include <RDBoost/import_array.h>
#include <RDBoost/PySequenceHolder.h>
#include <RDBoost/Wrap.h>
#include <Geometry/point.h>
#include <Numerics/Alignment/AlignPoints.h>
namespace python = boost::python;
namespace RDNumeric {
namespace Alignments {
void GetPointsFromPythonSequence(python::object &points,
RDGeom::Point3DConstPtrVect &pts) {
PyObject *pyObj = points.ptr();
unsigned int nrows, ncols;
double *data;
if (PyArray_Check(pyObj)) {
// get the dimensions of the array
PyArrayObject *ptsMat = reinterpret_cast<PyArrayObject *>(pyObj);
nrows = PyArray_DIM(ptsMat, 0);
ncols = PyArray_DIM(ptsMat, 1);
if (ncols != 3) throw_value_error("Wrong dimension for the points array");
data = reinterpret_cast<double *>(PyArray_DATA(ptsMat));
for (unsigned int i = 0; i < nrows; i++) {
auto *rpt =
new RDGeom::Point3D(data[i * 3], data[i * 3 + 1], data[i * 3 + 2]);
pts.push_back(rpt);
}
} else if (PySequence_Check(pyObj)) {
nrows = PySequence_Size(pyObj);
if (nrows <= 0) throw_value_error("Empty sequence passed in");
python::extract<RDGeom::Point3D> ptOk(points[0]);
if (!ptOk.check()) {
for (unsigned int i = 0; i < nrows; i++) {
PySequenceHolder<double> row(points[i]);
if (row.size() != 3)
throw_value_error("Wrong number of entries in the list of lists");
auto *rpt = new RDGeom::Point3D(row[0], row[1], row[2]);
pts.push_back(rpt);
}
} else {
for (unsigned int i = 0; i < nrows; i++) {
python::extract<RDGeom::Point3D> pt(points[i]);
if (pt.check()) {
auto *rpt = new RDGeom::Point3D(pt);
pts.push_back(rpt);
} else {
throw_value_error("non-Point3D found in sequence of points");
}
}
}
} else {
throw_value_error("non-sequence argument provided");
}
}
PyObject *AlignPointPairs(python::object refPoints, python::object probePoints,
python::object weights = python::list(),
bool reflect = false,
unsigned int maxIterations = 50) {
// The reference and probe points can be specifed in two formats
// 1. A Numeric array of dimensions (N,3) where N is the number of points
// 2. A list (or tuple) or lists (or tuples)
//
// The similar thing applies to weights
// 1. Can be a numeric vector of size N
// 2. A list of doubles of size N
// first deal with situation where we have Numerics arrays
RDGeom::Point3DConstPtrVect refPts, probePts;
GetPointsFromPythonSequence(refPoints, refPts);
GetPointsFromPythonSequence(probePoints, probePts);
unsigned int npt = refPts.size();
if (npt != probePts.size())
throw_value_error("Mis-match in number of points");
PyObject *weightsObj = weights.ptr();
RDNumeric::DoubleVector *wtsVec;
wtsVec = nullptr;
double *data;
if (PyArray_Check(weightsObj)) {
PyArrayObject *wtsMat = reinterpret_cast<PyArrayObject *>(weightsObj);
unsigned int nwts = PyArray_DIM(wtsMat, 0);
if (nwts != npt)
throw_value_error(
"Number of weights supplied do not match the number of points");
wtsVec = new RDNumeric::DoubleVector(nwts);
data = reinterpret_cast<double *>(PyArray_DATA(wtsMat));
for (unsigned int i = 0; i < nwts; i++) {
wtsVec->setVal(i, data[i]);
}
} else {
PySequenceHolder<double> wts(weights);
unsigned int nwts = wts.size();
if (nwts > 0) {
if (nwts != npt)
throw_value_error(
"Number of weights supplied do not match the number of points");
wtsVec = new RDNumeric::DoubleVector(nwts);
for (unsigned int i = 0; i < npt; i++) {
wtsVec->setVal(i, wts[i]);
}
}
}
RDGeom::Transform3D trans;
double ssd =
AlignPoints(refPts, probePts, trans, wtsVec, reflect, maxIterations);
npy_intp dims[2];
dims[0] = 4;
dims[1] = 4;
PyArrayObject *res = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_DOUBLE);
double *resData = reinterpret_cast<double *>(PyArray_DATA(res));
const double *tdata = trans.getData();
for (unsigned int i = 0; i < trans.numRows(); ++i) {
unsigned int itab = i * 4;
for (unsigned int j = 0; j < trans.numRows(); ++j) {
resData[itab + j] = tdata[itab + j];
}
}
delete wtsVec;
for (unsigned int i = 0; i < npt; ++i) {
delete probePts[i];
delete refPts[i];
}
PyObject *resTup = PyTuple_New(2);
PyObject *ssdItem = PyFloat_FromDouble(ssd);
PyTuple_SetItem(resTup, 0, ssdItem);
PyTuple_SetItem(resTup, 1, PyArray_Return(res));
return resTup;
}
}
}
BOOST_PYTHON_MODULE(rdAlignment) {
rdkit_import_array();
python::scope().attr("__doc__") =
"Module containing functions to align pairs of points in 3D";
std::string docString =
"Compute the optimal alignment (minimum RMSD) between two set of points \n\n\
\n\
ARGUMENTS:\n\n\
- refPoints : reference points sepcified as a N by 3 Numeric array or \n\
sequence of 3-sequences or sequence of Point3Ds \n\
- probePoints : probe points to align to reference points - same format \n\
restrictions as reference points apply here \n\
- weights : optional numeric vector or list of weights to associate to each pair of points\n\
- reflect : reflect the probe points before attempting alignment\n\
- maxIteration : maximum number of iterations to try to minimize RMSD \n\
\n\
RETURNS:\n\n\
a 2-tuple:\n\
- SSD value for the alignment\n\
- the 4x4 transform matrix, as a Numeric array\n\
\n";
python::def(
"GetAlignmentTransform", RDNumeric::Alignments::AlignPointPairs,
(python::arg("refPoints"), python::arg("probePoints"),
python::arg("weights") = python::list(), python::arg("reflect") = false,
python::arg("maxIterations") = 50),
docString.c_str());
}
|
