topo.base.sheetcoords

1 """ 2 Provides SheetCoordinateSystem, allowing conversion between the two 3 coordinate systems used in Topographica. 4 5 'Sheet coordinates' allow simulation parameters to be specified in 6 units that are density-independent, whereas 'matrix coordinates' 7 provide a means of realizing the continuous sheets. 8 9 Hence we can have a pair of 'sheet coordinates' (x,y); floating-point 10 Cartesian coordinates indicating an arbitrary point on the sheet's plane. 11 We can also have a pair of 'matrix coordinates' (r,c), which are used 12 to address an underlying matrix. These matrix coordinates are also 13 floating-point numbers to allow precise conversion between the two 14 schemes. Where it is necessary to address a specific element of the 15 matrix (as is often the case in calculations), we also have the usual 16 matrix index coordinates (r_idx, c_idx). We refer to these as 17 matrixidx coordinates. SheetCoordinateSystem provies methods for converting 18 between sheet and matrix coordinates, as well as sheet and matrixidx 19 coordinates. 20 21 Everyone should use these facilities for conversions between the two 22 coordinate systems to guarantee consistency. 23 24 25 26 27 Example of how the matrix stores the representation of the Sheet 28 ================================================================ 29 30 For the purposes of this example, assume the goal is a Sheet with 31 density=3 that has a 1 at (-1/2,-1/2), a 5 at (0.0,0.0), and a 9 at 32 (1/2,1/2). More precisely, for this Sheet, 33 34 the continuous area from -1/2,-1/2 to -1/6,-1/6 has value 1, 35 the continuous area from -1/6,-1/6 to 1/6,1/6 has value 5, and 36 the continuous area from 1/6,1/6 to 1/2,1/2 has value 9. 37 38 With the rest of the elements filled in, the Sheet would look like:: 39 40 (-1/2,1/2) -+-----+-----+-----+- (1/2,1/2) 41 | | | | 42 | 7 | 8 | 9 | 43 | | | | 44 (-1/2,1/6) -+-----+-----+-----+- (1/2,1/6) 45 | | | | 46 | 4 | 5 | 6 | 47 | | | | 48 (-1/2,-1/6) -+-----+-----+-----+- (1/2,-1/6) 49 | | | | 50 | 1 | 2 | 3 | 51 | | | | 52 (-1/2,-1/2) -+-----+-----+-----+- (1/2,-1/2) 53 54 where element 5 is centered on 0.0,0.0. A matrix that would match 55 these Sheet coordinates is:: 56 57 [[7 8 9] 58 [4 5 6] 59 [1 2 3]] 60 61 If we have such a matrix, we can access it in one of two ways: Sheet 62 or matrix/matrixidx coordinates. In matrixidx coordinates, the matrix is indexed 63 by rows and columns, and it is possible to ask for the element at 64 location [0,2] (which returns 9 as in any normal row-major matrix). 65 But the values can additionally be accessed in Sheet coordinates, 66 where the matrix is indexed by a point in the Cartesian plane. In 67 Sheet coordinates, it is possible to ask for the element at location 68 (0.3,0.02), which returns floating-point matrix coordinates that can 69 be cropped to give the nearest matrix element, namely the one with 70 value 6. 71 72 Of course, it would be an error to try to pass matrix coordinates like 73 [0,2] to the sheet2matrix calls; the result would be a value far 74 outside of the actual matrix. 75 76 $Id: sheetcoords.py 10613 2009-10-03 11:42:49Z ceball $ 77 """ 78 79 __version__ = '$Revision: 10613 $' 80 81 82 from numpy import array,floor,ceil,round_,arange 83 84 from boundingregion import BoundingBox

85 86 87 88 89 90 # Note about the 'bounds-master' approach we have adopted 91 # ======================================================= 92 # 93 # Our current approach is a "bounds-master" approach, where we trust 94 # the user's specified x width, and choose the nearest density and y 95 # height to make that possible. The advantage of this is that when we 96 # change the density (which is often), each such simulation is the 97 # best possible approximation to the given area. Generally, the area 98 # is much more meaningful than the density, so this approach makes 99 # sense. Plus, the y height is usually the same as the x width, so 100 # it's not usually a problem that the y height is not respected. The 101 # disadvantage is that the user's area can only be trusted in one 102 # dimension, because of wanting to avoid the complication of separate 103 # xdensity and ydensity, which makes this approach very difficult to 104 # explain to the user. 105 # 106 # The other approach is density-master: trust the user's specified 107 # density as-is, and simply choose the nearest area that fits that 108 # density. The advantages are that (a) it's very simple to describe 109 # and reason about, and (b) if a user asks for a different area, they 110 # get a true subsection of the same simulation they would have gotten 111 # at the larger area. The disadvantage is that the simulation isn't 112 # the best approximation of the given area that it could be -- e.g. at 113 # low densities, the sheet area could be quite significantly different 114 # than the one the user requested. Plus, if we took this approach 115 # seriously, then we'd let the density specify the matrix coordinate 116 # system entirely, including the origin, which would mean that the 117 # actual area would often be offset from the intended one, which is 118 # even worse. Differences between the area and the offset could cause 119 # severe problems in the alignment of projections between sheets with 120 # different densities, which would make low-density versions of 121 # hierarchical models behave very strangely. 122 123 124 125 126 -class SheetCoordinateSystem(object):

127 """ 128 Provides methods to allow conversion between sheet and matrix 129 coordinates. 130 """

131 - def __get_xdensity(self):

132 return self.__xdensity

133 - def __get_ydensity(self):

134 return self.__ydensity

135 - def __get_shape(self):

136 return self.__shape

137 138 #### These are all properties so that they can't be set #### 139 # CB: unnecessary? if we're going to keep this, what about bounds and lbrt? 140 xdensity = property(__get_xdensity, 141 doc="""The spacing between elements in an underlying 142 matrix representation, in the x direction.""") 143 144 ydensity = property(__get_ydensity, 145 doc="""The spacing between elements in an underlying 146 matrix representation, in the y direction.""") 147 148 shape = property(__get_shape) 149 150

151 - def __init__(self,bounds,xdensity,ydensity=None):

152 """ 153 Store the bounds (as l,b,r,t in an array), xdensity, and 154 ydensity. 155 156 If ydensity is not specified, it is assumed that the specified 157 xdensity is nominal and that the true xdensity should be 158 calculated. The top and bottom bounds are adjusted so that 159 the ydensity is equal to the xdensity. 160 161 If both xdensity and ydensity are specified, these and the bounds 162 are taken to be exact and are not adjusted. 163 """ 164 if not ydensity: 165 bounds,xdensity = self.__equalize_densities(bounds,xdensity) 166 167 self.bounds = bounds 168 self.__set_xdensity(xdensity) 169 self.__set_ydensity(ydensity or xdensity) 170 171 self.lbrt = array(bounds.lbrt()) 172 173 r1,r2,c1,c2 = Slice._boundsspec2slicespec(self.lbrt,self) 174 self.__shape = (r2-r1,c2-c1)

175 176 177 ### we use xstep and ystep so that the repeatedly performed 178 ### calculations in matrix2sheet() use multiplications rather than 179 ### divisions, for speed

180 - def __set_xdensity(self,density):

181 self.__xdensity=density 182 self.__xstep = 1.0/density

183

184 - def __set_ydensity(self,density):

185 self.__ydensity=density 186 self.__ystep = 1.0/density

187 188

189 - def __equalize_densities(self,nominal_bounds,nominal_density):

190 """ 191 Calculate the true density along x, and adjust the top and 192 bottom bounds so that the density along y will be equal. 193 194 Returns (adjusted_bounds,true_density) 195 """ 196 left,bottom,right,top = nominal_bounds.lbrt() 197 width = right-left; height = top-bottom 198 center_y = bottom + height/2.0 199 # The true density is not equal to the nominal_density 200 # when nominal_density*(right-left) is not an integer. 201 true_density = int(nominal_density*(width))/float(width) 202 203 n_cells = round(height*true_density,0) 204 adjusted_half_height = n_cells/true_density/2.0 205 # (The above might be clearer as (step*n_units)/2.0, where 206 # step=1.0/density.) 207 208 return (BoundingBox(points=((left, center_y-adjusted_half_height), 209 (right, center_y+adjusted_half_height))), 210 true_density)

211 212

213 - def sheet2matrix(self,x,y):

214 """ 215 Convert a point (x,y) in Sheet coordinates to continuous matrix 216 coordinates. 217 218 Returns (float_row,float_col), where float_row corresponds to y, 219 and float_col to x. 220 221 Valid for scalar or array x and y. 222 223 Note about Bounds 224 For a Sheet with BoundingBox(points=((-0.5,-0.5),(0.5,0.5))) 225 and density=3, x=-0.5 corresponds to float_col=0.0 and x=0.5 226 corresponds to float_col=3.0. float_col=3.0 is not inside the 227 matrix representing this Sheet, which has the three columns 228 (0,1,2). That is, x=-0.5 is inside the BoundingBox but x=0.5 229 is outside. Similarly, y=0.5 is inside (at row 0) but y=-0.5 230 is outside (at row 3) (it's the other way round for y because 231 the matrix row index increases as y decreases). 232 """ 233 # First translate to (left,top), which is [0,0] in the matrix, 234 # then scale to the size of the matrix. The y coordinate needs 235 # to be flipped, because the points are moving down in the 236 # sheet as the y index increases in the matrix. 237 float_col = (x-self.lbrt[0]) * self.__xdensity 238 float_row = (self.lbrt[3]-y) * self.__ydensity 239 return float_row, float_col

240 241

242 - def sheet2matrixidx(self,x,y):

243 """ 244 Convert a point (x,y) in sheet coordinates to the integer row and 245 column index of the matrix cell in which that point falls, given a 246 bounds and density. Returns (row,column). 247 248 Note that if coordinates along the right or bottom boundary are 249 passed into this function, the returned matrix coordinate of the 250 boundary will be just outside the matrix, because the right and 251 bottom boundaries are exclusive. 252 253 Valid for scalar or array x and y. 254 """ 255 r,c = self.sheet2matrix(x,y) 256 r = floor(r) 257 c = floor(c) 258 259 # CB: was it better to have two different methods? 260 if hasattr(r,'astype'): 261 return r.astype(int), c.astype(int) 262 else: 263 return int(r),int(c)

264 265

266 - def matrix2sheet(self,float_row,float_col):

267 """ 268 Convert a floating-point location (float_row,float_col) in matrix 269 coordinates to its corresponding location (x,y) in sheet 270 coordinates. 271 272 Valid for scalar or array float_row and float_col. 273 274 Inverse of sheet2matrix(). 275 """ 276 x = float_col*self.__xstep + self.lbrt[0] 277 y = self.lbrt[3] - float_row*self.__ystep 278 return x, y

279 280

281 - def matrixidx2sheet(self,row,col):

282 """ 283 Return (x,y) where x and y are the floating point coordinates of 284 the *center* of the given matrix cell (row,col). If the matrix cell 285 represents a 0.2 by 0.2 region, then the center location returned 286 would be 0.1,0.1. 287 288 NOTE: This is NOT the strict mathematical inverse of 289 sheet2matrixidx(), because sheet2matrixidx() discards all but the integer 290 portion of the continuous matrix coordinate. 291 292 Valid only for scalar or array row and col. 293 """ 294 x,y = self.matrix2sheet((row+0.5), (col+0.5)) 295 296 # Rounding is useful for comparing the result with a floating point number 297 # that we specify by typing the number out (e.g. fp = 0.5). 298 # Round eliminates any precision errors that have been compounded 299 # via floating point operations so that the rounded number will better 300 # match the floating number that we type in. 301 return round_(x,10),round_(y,10)

302 303

304 - def closest_cell_center(self,x,y):

305 """ 306 Given arbitary sheet coordinates, return the sheet coordinates 307 of the center of the closest unit. 308 """ 309 return self.matrixidx2sheet(*self.sheet2matrixidx(x,y))

310 311

312 - def sheetcoordinates_of_matrixidx(self):

313 """ 314 Return x,y where x is a vector of sheet coordinates 315 representing the x-center of each matrix cell, and y 316 represents the corresponding y-center of the cell. 317 """ 318 rows,cols = self.shape 319 return self.matrixidx2sheet(arange(rows),arange(cols))

320 321 322 323 324 325 326 327 # Needs cleanup/rename: 328 # 329 # since it's different from slice. It's our special slice that's an 330 # array specifying row_start,row_stop,col_start,col_stop for a 331 # Sheet (2d array). 332 # 333 # In python, a[slice(0,2)] (where a is a list/array/similar) is 334 # equivalent to a[0:2]. 335 # 336 # So, not sure what to call this class. (Will need to rename some 337 # methods, too.) SCSSlice? SheetSlice? 338 339 from numpy import int32,ndarray

340 341 -class Slice(ndarray):

342 """ 343 Represents a slice of a SheetCoordinateSystem; i.e., an array 344 specifying the row and column start and end points for a submatrix 345 of the SheetCoordinateSystem. 346 347 The slice is created from the supplied bounds by calculating the 348 slice that corresponds most closely to the specified bounds. 349 Therefore, the slice does not necessarily correspond exactly to 350 the specified bounds. The bounds that do exactly correspond to the 351 slice are available via the 'bounds' attribute. 352 353 Note that the slice does not respect the bounds of the 354 SheetCoordinateSystem, and that actions such as translate() also 355 do not respect the bounds. To ensure that the slice is within the 356 SheetCoordinateSystem's bounds, use crop_to_sheet(). 357 """ 358

359 - def compute_bounds(self,scs):

360 spec = self._slicespec2boundsspec(self,scs) 361 return BoundingBox(points=spec)

362 363 __slots__ = [] 364

365 - def __new__(cls, bounds, sheet_coordinate_system, force_odd=False, 366 min_matrix_radius=1): # CEBALERT: min_matrix_radius only used for odd slice.

367 """ 368 Create a slice of the given sheet_coordinate_system from the 369 specified bounds. 370 """ 371 # I couldn't find documentation on subclassing array; I used 372 # the following as reference: 373 # http://scipy.org/scipy/numpy/browser/branches/maskedarray/ 374 # numpy/ma/tests/test_subclassing.py?rev=4577 375 if force_odd: 376 slicespec=Slice._createoddslicespec(bounds,sheet_coordinate_system, 377 min_matrix_radius) 378 else: 379 slicespec=Slice._boundsspec2slicespec(bounds.lbrt(),sheet_coordinate_system) 380 # numpy.int32 is specified explicitly in Slice to avoid having 381 # it default to numpy.int. int32 saves memory (and is expected 382 # by optimized C functions). 383 a = array(slicespec, dtype=int32, copy=False).view(cls) 384 return a

385 386

387 - def submatrix(self,matrix):

388 """ 389 Return the submatrix of the given matrix specified by this 390 slice. 391 392 Equivalent to computing the intersection between the 393 SheetCoordinateSystem's bounds and the bounds, and 394 returning the corresponding submatrix of the given matrix. 395 396 The submatrix is just a view into the sheet_matrix; it is not 397 an independent copy. 398 """ 399 return matrix[self[0]:self[1],self[2]:self[3]]

400 401 # CB: not sure if this is a good idea or not. In some ways, I 402 # think matrix[slice_()] would be clearer than 403 # slice.submatrix(matrix), but I'm not sure. cf.py is where 404 # to look to see this in action. 405 ## def __call__(self): 406 ## return slice(self[0],self[1]),slice(self[2],self[3]) 407 408 ### CLEANUP ### 409 410 # CEBALERT: unnecessary? use translate and crop and back. or is 411 # that more steps? rename

412 - def positionlesscrop(self,x,y,sheet_coord_system):

413 """ 414 Return the correct slice for a weights/mask matrix at this 415 ConnectionField's location on the sheet (i.e. for getting 416 the correct submatrix of the weights or mask in case the 417 unit is near the edge of the sheet). 418 """ 419 sheet_rows,sheet_cols = sheet_coord_system.shape 420 421 # get size of weights matrix 422 n_rows,n_cols = self.shape_on_sheet() 423 424 # get slice for the submatrix 425 center_row,center_col = sheet_coord_system.sheet2matrixidx(x,y) 426 427 c1 = -min(0, center_col-n_cols/2) # assume odd weight matrix so can use n_cols/2 428 r1 = -min(0, center_row-n_rows/2) # for top and bottom 429 c2 = -max(-n_cols, center_col-sheet_cols-n_cols/2) 430 r2 = -max(-n_rows, center_row-sheet_rows-n_rows/2) 431 432 self.set((r1,r2,c1,c2))

433 434 # CBALERT: assumes the user wants the bounds to be centered about 435 # the unit, which might not be true.

436 - def positionedcrop(self,x,y,sheet_coord_system):

437 """ 438 Offset the bounds_template to this cf's location and store the 439 result in the 'bounds' attribute. 440 441 Also stores the input_sheet_slice for access by C. 442 """ 443 # translate to this cf's location 444 cf_row,cf_col = sheet_coord_system.sheet2matrixidx(x,y) 445 446 # should result in same no. of comps of right bounds as before during init, 447 # since i removed one calc from init 448 449 bounds_x,bounds_y=self.compute_bounds(sheet_coord_system).centroid() 450 451 b_row,b_col=sheet_coord_system.sheet2matrixidx(bounds_x,bounds_y) 452 453 row_offset = cf_row-b_row 454 col_offset = cf_col-b_col 455 self.translate(row_offset,col_offset)

456 457 ############### 458 459

460 - def translate(self, r, c):

461 """ 462 Translate the slice by the specified number of rows 463 and columns. 464 """ 465 self+=[r,r,c,c]

466

467 - def set(self,slice_specification):

468 """Set this slice from some iterable that specifies (r1,r2,c1,c2).""" 469 self.put([0,1,2,3],slice_specification) # pylint: disable-msg=E1101

470

471 - def shape_on_sheet(self):

472 """Return the shape of the array that this Slice would give on its sheet.""" 473 return self[1]-self[0],self[3]-self[2]

474

475 - def crop_to_sheet(self,sheet_coord_system):

476 """Crop the slice to the SheetCoordinateSystem's bounds.""" 477 maxrow,maxcol = sheet_coord_system.shape 478 479 self[0] = max(0,self[0]) 480 self[1] = min(maxrow,self[1]) 481 self[2] = max(0,self[2]) 482 self[3] = min(maxcol,self[3]) 483 484 485 ### CB: working on methods below here 486 # (+ not keeping these names) 487 488 489 @staticmethod

490 - def _createoddslicespec(bounds,scs,min_matrix_radius):

491 """ 492 Create the 'odd' Slice that best approximates the specifed 493 sheet-coordinate bounds. 494 495 The supplied bounds are translated to have a center at the 496 center of one of the sheet's units (we arbitrarily use the 497 center unit), and then these bounds are converted to a slice 498 in such a way that the slice exactly includes all units whose 499 centers are within the bounds (see boundsspec2slicespec()). 500 However, to ensure that the bounds are treated symmetrically, 501 we take the right and bottom bounds and reflect these about 502 the center of the slice (i.e. we take the 'xradius' to be 503 right_col-center_col and the 'yradius' to be 504 bottom_col-center_row). Hence, if the bounds happen to go 505 through units, if the units are included on the right and 506 bottom bounds, they will be included on the left and top 507 bounds. This ensures that the slice has odd dimensions. 508 """ 509 bounds_xcenter,bounds_ycenter=bounds.centroid() 510 sheet_rows,sheet_cols = scs.shape 511 512 # arbitrary (e.g. could use 0,0) 513 center_row,center_col = sheet_rows/2,sheet_cols/2 514 unit_xcenter,unit_ycenter=scs.matrixidx2sheet(center_row, 515 center_col) 516 517 bounds.translate(unit_xcenter-bounds_xcenter, 518 unit_ycenter-bounds_ycenter) 519 520 ########## CEBALERT: assumes weights are to be centered about each unit. 521 r1,r2,c1,c2 = Slice._boundsspec2slicespec(bounds.lbrt(),scs) 522 523 # use the calculated radius unless it's smaller than the min 524 xrad=max(c2-center_col-1,min_matrix_radius) 525 yrad=max(r2-center_row-1,min_matrix_radius) 526 527 r2=center_row+yrad+1 528 c2=center_col+xrad+1 529 r1=center_row-yrad 530 c1=center_col-xrad 531 ########## 532 533 # weights matrix must be odd (otherwise this method has an error) 534 # CEBALERT: this test should move to a test file. 535 ## if rows%2!=1 or cols%2!=1: 536 ## raise AssertionError("nominal_bounds_template yielded even-height or even-width weights matrix (%s rows, %s columns) - weights matrix must have odd dimensions."%(rows,cols)) 537 538 return (r1,r2,c1,c2)

539 540 541 ## # CEBALERT: with min_matrix_radius, this test is unnecessary? (check) 542 ## # user-supplied bounds must lead to a weights matrix of at least 1x1 543 ## rows,cols = weights_slice.shape_on_sheet() 544 ## if rows==0 or cols==0: 545 ## raise ValueError("nominal_bounds_template results in a zero-sized weights matrix (%s,%s) - you may need to supply a larger nominal_bounds_template or increase the density of the sheet."%(rows,cols)) 546 547 548 549 550 @staticmethod

551 - def _boundsspec2slicespec(boundsspec,scs):

552 """ 553 Convert an iterable boundsspec (supplying l,b,r,t of a 554 BoundingRegion) into a Slice specification. 555 556 Includes all units whose centers are within the specified 557 sheet-coordinate bounds specified by boundsspec. 558 559 Exact inverse of _slicespec2boundsspec(). 560 """ 561 l,b,r,t = boundsspec 562 563 t_m,l_m = scs.sheet2matrix(l,t) 564 b_m,r_m = scs.sheet2matrix(r,b) 565 566 l_idx = int(ceil(l_m-0.5)) 567 t_idx = int(ceil(t_m-0.5)) 568 # CBENHANCEMENT: Python 2.6's math.trunc()? 569 r_idx = int(floor(r_m+0.5)) 570 b_idx = int(floor(b_m+0.5)) 571 572 return t_idx,b_idx,l_idx,r_idx

573 574 575 @staticmethod

576 - def _slicespec2boundsspec(slicespec,scs):

577 """ 578 Convert an iterable slicespec (supplying r1,r2,c1,c2 of a 579 Slice) into a BoundingRegion specification. 580 581 Exact inverse of _boundsspec2slicespec(). 582 """ 583 r1,r2,c1,c2 = slicespec 584 585 left,bottom = scs.matrix2sheet(r2,c1) 586 right, top = scs.matrix2sheet(r1,c2) 587 588 return ((left,bottom),(right,top))

589

Source Code for Module topo.base.sheetcoords