root/trunk/src/space.py@778

Revision 778, 11.2 KB (checked in by apdavison, 3 years ago)

Starting to fill out the space module; minor tweaks from running doctest on highlevelapi.txt

Line
1# encoding: utf-8
2"""
3Tools for performing spatial/topographical calculations.
4
5"""
6
7# There must be some Python package out there that provides most of this stuff.
8
9import numpy
10import math
11from operator import and_
12from pyNN.random import NumpyRNG
13
14def distance(src, tgt, mask=None, scale_factor=1.0, offset=0.0,
15             periodic_boundaries=None): # may need to add an offset parameter
16    """
17    Return the Euclidian distance between two cells.
18    `mask` allows only certain dimensions to be considered, e.g.::
19      * to ignore the z-dimension, use `mask=array([0,1])`
20      * to ignore y, `mask=array([0,2])`
21      * to just consider z-distance, `mask=array([2])`
22    `scale_factor` allows for different units in the pre- and post- position
23    (the post-synaptic position is multipied by this quantity).
24    """
25    d = src.position - scale_factor*(tgt.position + offset)
26
27    if not periodic_boundaries == None:
28        d = numpy.minimum(abs(d), periodic_boundaries-abs(d))
29    if mask is not None:
31    return numpy.sqrt(numpy.dot(d, d))
32
33
34class Space(object):
35
36    AXES = {'x' : [0],    'y': [1],    'z': [2],
37            'xy': [0,1], 'yz': [1,2], 'xz': [0,2], 'xyz': range(3), None: range(3)}
38
39    def __init__(self, axes=None, scale_factor=1.0, offset=0.0,
40                 periodic_boundaries=None):
41        """
42        axes -- if not supplied, then the 3D distance is calculated. If supplied,
43            axes should be a string containing the axes to be used, e.g. 'x', or
44            'yz'. axes='xyz' is the same as axes=None.
45        scale_factor -- it may be that the pre and post populations use
46            different units for position, e.g. degrees and Âµm. In this case,
47            `scale_factor` can be specified, which is applied to the positions
48            in the post-synaptic population.
49        offset -- if the origins of the coordinate systems of the pre- and post-
50            synaptic populations are different, `offset` can be used to adjust
51            for this difference. The offset is applied before any scaling.
52        periodic_boundaries -- either `None`, or a tuple giving the boundaries
53            for each dimension, e.g. `((x_min, x_max), None, (z_min, z_max))`.
54        """
55        self.periodic_boundaries = periodic_boundaries
56        self.axes = numpy.array(Space.AXES[axes])
57        self.scale_factor = scale_factor
58        self.offset = offset
59
60    def distances(self, A, B, expand=False):
61        """
62        Calculate the distance matrix between two sets of coordinates, given
63        the topology of the current space.
64        From http://projects.scipy.org/pipermail/numpy-discussion/2007-April/027203.html
65        """
66        if len(A.shape) == 1:
67            A = A.reshape(3, 1)
68        if len(B.shape) == 1:
69            B = B.reshape(3, 1)
70        # I'm not sure the following line should be here. Operations may be redundant and not very
71        # transparent from the user point of view. I moved it into the DistanceDependentProbability Connector
72        #B = self.scale_factor*(B + self.offset)
73        d = numpy.zeros((len(self.axes), A.shape[1], B.shape[1]), dtype=A.dtype)
74        for axis in self.axes:
75            diff2 = A[axis,:,None] - B[axis, :]
76            if self.periodic_boundaries is not None:
77                boundaries = self.periodic_boundaries[axis]
78                if boundaries is not None:
79                    range = boundaries[1]-boundaries[0]
82            diff2 **= 2
83            d[axis] = diff2
84        if not expand:
85            d = numpy.sum(d, 0)
86        numpy.sqrt(d, d)
87        return d
88
89
90
91class PositionsGenerator(object):
92
93        def __init__(self, dimensions, axes=None):
94            """
95            dimensions -- either `None`, or a tuple/list giving the dimensions
96                          for each dimension, e.g. `((x_min, x_max), None, (z_min, z_max))`.
97            axes -- if not supplied, then the 3D distance is calculated. If supplied,
98                    axes should be a string containing the axes to be used, e.g. 'x', or
99                    'yz'. axes='xyz' is the same as axes=None.
100            """
101            self.dimensions = list(dimensions)
102            self.axes = numpy.array(Space.AXES[axes])
103            assert len(self.dimensions) == 3, "Dimensions should be of size 3, and axes should specify orientation!"
104            for item in self.dimensions:
105                if item is not None:
106                    assert len(item) == 2, "dimensions should be a list of tuples (min, max), not %s" %item
107                    assert item[0] <= item[0], "items elements should be (min, max), with min <= max"
108
109        def get(self, dims):
110            self.M         = numpy.prod(numpy.array(dims))
111            self.positions = numpy.zeros((3, self.M))
112            pass
113
114
115class RandomPositions(PositionsGenerator):
116
117        def __init__(self, dimensions, seed=None):
118            PositionsGenerator.__init__(self, dimensions)
119            self.seed = seed
120
121        def get(self, dims):
122            PositionsGenerator.get(self, dims)
123            numpy.random.seed(self.seed)
124            for axis in self.axes:
125                item = self.dimensions[axis]
126                if item is not None:
127                    bmin, bmax              = item
128                    self.positions[axis, :] = bmin + bmax * numpy.random.rand(self.M)
129            return self.positions
130
131
132class GridPositions(PositionsGenerator):
133
134        def get(self, dims):
135            PositionsGenerator.get(self, dims)
136            res = ""
137            for d in dims:
138                res += "0:%d," %d
139            grid = eval("numpy.mgrid[%s]" %res[0:-1])
140            for axis in self.axes:
141                item = self.dimensions[axis]
142                if item is not None:
143                    bmin, bmax = item
145                    data       = numpy.linspace(bmin+padding, bmax, dims[axis])
146                    self.positions[axis, :] = data[grid[axis].flatten()]
147            return self.positions
148
149
150class BaseStructure(object):
151    pass
152
153
154class Line(BaseStructure):
155
156    def __init__(self, dx=1.0, x0=0.0, y0=0.0, z0=0.0):
157        self.dx = dx
158        self.x0 = x0
159        self.y0 = y0
160        self.z0 = z0
161
162    def __eq__(self, other):
163        return reduce(and_, (getattr(self, attr) == getattr(other, attr)
164                             for attr in ("dx", "x0", "y0", "z0")))
165
166    def generate_positions(self, n):
167        x = self.dx*numpy.arange(n, dtype=float) + self.x0
168        y = numpy.zeros(n) + self.y0
169        z = numpy.zeros(n) + self.z0
170        return numpy.array((x,y,z))
171
172    def describe(self, n):
173        return "line with %d positions" % n
174
175
176class Grid2D(BaseStructure):
177
178    def __init__(self, aspect_ratio=1.0, dx=1.0, dy=1.0, x0=0.0, y0=0.0, z=0, fill_order="sequential"):
179        """
180        aspect_ratio - ratio of the number of grid points per side (not the ratio
181                       of the side lengths, unless dx == dy)
182        """
183        self.aspect_ratio = aspect_ratio
184        assert fill_order in ('sequential', 'random')
185        self.fill_order = fill_order
186        self.dx = dx; self.dy = dy; self.x0 = x0; self.y0 = y0; self.z = z
187
188    def __eq__(self, other):
189        return reduce(and_, (getattr(self, attr) == getattr(other, attr)
190                             for attr in ("aspect_ratio", "dx", "dy",
191                                          "x0", "y0", "fill_order")))
192
193    def calculate_size(self, n):
194        nx = math.sqrt(n*self.aspect_ratio)
195        ny = n/nx
196        return nx, ny
197
198    def generate_positions(self, n):
199        nx, ny = self.calculate_size(n)
200        x,y,z = numpy.indices((nx,ny,1), dtype=float)
201        x = self.x0 + self.dx*x.flatten()
202        y = self.y0 + self.dy*y.flatten()
203        z = self.z + z.flatten()
204        if self.fill_order == 'sequential':
205            return numpy.array((x,y,z)) # use column_stack, if we decide to switch from (3,n) to (n,3)
206        else:
207            raise NotImplementedError
208
209    def describe(self, n):
210        return "2D grid of size (%d, %d)" % self.calculate_size(n)
211
212
213class Grid3D(BaseStructure):
214
215    def __init__(self, aspect_ratioXY, aspect_ratioXZ, dx=1.0, dy=1.0, dz=1.0, x0=0.0, y0=0.0, z0=0,
216                 fill_order="sequential"):
217        """
218        If fill_order is 'sequential', the z-index will be filled first, then y then x, i.e.
219        the first cell will be at (0,0,0) (given default values for the other arguments),
220        the second at (0,0,1), etc.
221        """
222        self.aspect_ratios = (aspect_ratioXY, aspect_ratioXZ)
223        assert fill_order in ('sequential', 'random')
224        self.fill_order = fill_order
225        self.dx = dx; self.dy = dy; self.dz = dz
226        self.x0 = x0; self.y0 = y0; self.z0 = z0
227
228    def __eq__(self, other):
229        return reduce(and_, (getattr(self, attr) == getattr(other, attr)
230                             for attr in ("aspect_ratios", "dx", "dy", "dz",
231                                          "x0", "y0", "z0", "fill_order")))
232
233    def calculate_size(self, n):
234        a,b = self.aspect_ratios
235        nx = int(round(math.pow(n*a*b, 1/3.0)))
236        ny = int(round(nx/a))
237        nz = int(round(nx/b))
238        assert nx*ny*nz == n, str((nx, ny, nz, nx*ny*nz, n, a, b))
239        return nx, ny, nz
240
241    def generate_positions(self, n):
242        nx, ny, nz = self.calculate_size(n)
243        x,y,z = numpy.indices((nx,ny,nz), dtype=float)
244        x = self.x0 + self.dx*x.flatten()
245        y = self.y0 + self.dy*y.flatten()
246        z = self.z0 + self.dz*z.flatten()
247        if self.fill_order == 'sequential':
248            return numpy.array((x,y,z))
249        else:
250            raise NotImplementedError
251
252    def describe(self, n):
253        return "3D grid of size (%d, %d, %d)" % self.calculate_size(n)
254
255class Shape(object):
256    pass
257
258class Cuboid(Shape):
259
260    def __init__(self, height, width, depth):
261        """Not sure how h,w,d should be related to x,y,z."""
262        self.height = height
263        self.width = width
264        self.depth = depth
265
266    def sample(self, n, rng):
267        return rng.uniform(size=(n,3)) * (self.height, self.width, self.depth)
268
269
270class Sphere(Shape):
271
273        Shape.__init__(self)
275
276    def sample(self, n, rng):
277        raise NotImplementedError
278        #theta_phi = rng.uniform(size=(n,2))
279        #r = rng.???
280        # now transform from r,theta,phi to x,y,z
281
282
283class RandomStructure(BaseStructure):
284
285    def __init__(self, boundary, origin=(0.0,0.0,0.0), rotation=(0,0), rng=None):
286        """
287        `boundary` - a subclass of Shape
288        `origin` - a coordinate tuple (x,y,z)
289        `rotation` - (theta, phi) tuple giving the rotation to be applied (in degrees)
290
291        Note that the rotation is applied before the origin shift
292        """
293        self.boundary = boundary
294        self.origin = origin
295        self.rng = rng or NumpyRNG()
296
297    def generate_positions(self, n):
298        return numpy.array(self.origin) + rotate(self.boundary.sample(n, rng), *self.origin)
299
Note: See TracBrowser for help on using the browser.