summaryrefslogtreecommitdiff
path: root/cad/src/geometry/NeighborhoodGenerator.py
blob: 14e6e181148a11ca92a4755da331d6f4ac3f7d37 (plain)
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
# Copyright 2004-2008 Nanorex, Inc.  See LICENSE file for details.
"""
NeighborhoodGenerator.py -- linear time way to find overlapping or nearby atoms.

@author: Will
@version: $Id$
@copyright: 2004-2008 Nanorex, Inc.  See LICENSE file for details.

History:

Will wrote this in bonds.py for use in bond inference.

Bruce 080411 moved it into its own module, classifying it
as geometry (i.e. moving it into that source package)
since it's essentially purely geometric (and could easily
be generalized to be entirely so) -- all it assumes about
"atoms" is that they have a few methods like .posn()
and .is_singlet(), and it needs no imports from model.
"""

import struct

from Numeric import floor

from geometry.VQT import vlen

# ==


# This is an order(N) operation that produces a function which gets a
# list of potential neighbors in order(1) time. This is handy for
# inferring bonds for PDB files that lack any bonding information.
class NeighborhoodGenerator:
    """
    Given a list of atoms and a radius, be able to quickly take a
    point and generate a neighborhood, which is a list of the atoms
    within that radius of the point.

    Reading in the list of atoms is done in O(n) time, where n is the
    number of atoms in the list. Generating a neighborhood around a
    point is done in O(1) time. So this is a pretty efficient method
    for finding neighborhoods, especially if the same generator can
    be used many times.
    """
    def __init__(self, atomlist, maxradius, include_singlets = False):
        self._buckets = { }
        self._oldkeys = { }
        self._maxradius = 1.0 * maxradius
        self.include_singlets = include_singlets
        for atom in atomlist:
            self.add(atom)

    def _quantize(self, vec):
        """
        Given a point in space, partition space into little cubes
        so that when the time comes, it will be quick to locate and
        search the cubes right around a point.
        """
        maxradius = self._maxradius
        return (int(floor(vec[0] / maxradius)),
                int(floor(vec[1] / maxradius)),
                int(floor(vec[2] / maxradius)))

    def add(self, atom, _pack = struct.pack):
        buckets = self._buckets
        if self.include_singlets or not atom.is_singlet():
            # The keys of the _buckets dictionary are 12-byte strings.
            # Comparing them when doing lookups should be quicker than
            # comparisons between tuples of integers, so dictionary
            # lookups will go faster.
            i, j, k = self._quantize(atom.posn())
            key = _pack('lll', i, j, k)
            buckets.setdefault(key, []).append(atom)
            self._oldkeys[atom.key] = key

    def atom_moved(self, atom):
        """
        If an atom has been added to a neighborhood generator and
        is later moved, this method must be called to refresh the
        generator's position information. This only needs to be done
        during the useful lifecycle of the generator.
        """
        oldkey = self._oldkeys[atom.key]
        self._buckets[oldkey].remove(atom)
        self.add(atom)

    def region(self, center, _pack = struct.pack):
        """
        Given a position in space, return the list of atoms that
        are within the neighborhood radius of that position.
        """
        buckets = self._buckets
        def closeEnough(atm, radius = self._maxradius):
            return vlen(atm.posn() - center) < radius
        lst = [ ]
        x0, y0, z0 = self._quantize(center)
        for x in range(x0 - 1, x0 + 2):
            for y in range(y0 - 1, y0 + 2):
                for z in range(z0 - 1, z0 + 2):
                    # keys are 12-byte strings, see rationale above
                    key = _pack('lll', x, y, z)
                    if buckets.has_key(key):
                        lst += filter(closeEnough, buckets[key])
        return lst

    def remove(self, atom):
        key = self._quantize(atom.posn())
        try:
            self._buckets[key].remove(atom)
        except:
            pass

    pass # end of class NeighborhoodGenerator

# end