Source Code for Module topo.pattern.basic

   1  from __future__ import with_statement 
   2  """ 
   3  Simple two-dimensional mathematical or geometrical pattern generators. 
   4   
   5  $Id: basic.py 11327 2010-07-31 10:43:18Z ceball $ 
   6  """ 
   7  __version__='$Revision: 11327 $' 
   8   
   9  from math import pi, sqrt 
  10   
  11  import numpy 
  12  from numpy.oldnumeric import around,bitwise_and,sin,cos,bitwise_or 
  13  from numpy import asarray, float32, nonzero, zeros, shape, hstack,\ 
  14      linspace, abs, round, fft, alltrue 
  15   
  16  import param 
  17  from param.parameterized import ParamOverrides,as_uninitialized 
  18   
  19  import topo 
  20  # Imported here so that all PatternGenerators will be in the same package 
  21  from topo.base.patterngenerator import Constant, PatternGenerator 
  22   
  23  from topo.base.arrayutil import wrap 
  24  from topo.base.sheetcoords import SheetCoordinateSystem 
  25  from topo.misc.patternfn import gaussian,exponential,gabor,line,disk,ring,\ 
  26      sigmoid,arc_by_radian,arc_by_center,smooth_rectangle,float_error_ignore 
  27   
  28  from topo import numbergen 
  29   
  30   
  31  # Could add a Gradient class, where the brightness varies as a 
  32  # function of an equation for a plane.  This could be useful as a 
  33  # background, or to see how sharp a gradient is needed to get a 
  34  # response. 
  35   
  36  # CEBALERT: do we need this? If so, please remove this question. 
  37 -class Null(Constant): 
  38      """ 
  39      A constant pattern of zero activity. 
  40      """ 
  41      scale = param.Number(default=0,constant=True,precedence=-1) 
  42   
  43   
  44 -class HalfPlane(PatternGenerator): 
  45      """ 
  46      Constant pattern on in half of the plane, and off in the rest, 
  47      with optional Gaussian smoothing. 
  48      """ 
  49       
  50      smoothing = param.Number(default=0.02,bounds=(0.0,None),softbounds=(0.0,0.5), 
  51                               precedence=0.61,doc="Width of the Gaussian fall-off.") 
  52   
  53 -    def function(self,p): 
  54          if p.smoothing==0.0: 
  55              falloff=self.pattern_y*0.0 
  56          else: 
  57              with float_error_ignore(): 
  58                  falloff=numpy.exp(numpy.divide(-self.pattern_y*self.pattern_y, 
  59                                                  2*p.smoothing*p.smoothing)) 
  60   
  61          return numpy.where(self.pattern_y>0.0,1.0,falloff) 
  62   
  63   
  64 -class Gaussian(PatternGenerator): 
  65      """ 
  66      2D Gaussian pattern generator. 
  67   
  68      The sigmas of the Gaussian are calculated from the size and 
  69      aspect_ratio parameters: 
  70   
  71        ysigma=size/2 
  72        xsigma=ysigma*aspect_ratio 
  73   
  74      The Gaussian is then computed for the given (x,y) values as:: 
  75       
  76        exp(-x^2/(2*xsigma^2) - y^2/(2*ysigma^2) 
  77      """ 
  78       
  79      aspect_ratio = param.Number(default=1/0.31,bounds=(0.0,None),softbounds=(0.0,6.0), 
  80          precedence=0.31,doc=""" 
  81          Ratio of the width to the height. 
  82          Specifically, xsigma=ysigma*aspect_ratio (see size).""") 
  83       
  84      size = param.Number(default=0.155,doc=""" 
  85          Overall size of the Gaussian, defined by: 
  86          exp(-x^2/(2*xsigma^2) - y^2/(2*ysigma^2) 
  87          where ysigma=size/2 and xsigma=size/2*aspect_ratio.""") 
  88   
  89 -    def function(self,p): 
  90          ysigma = p.size/2.0 
  91          xsigma = p.aspect_ratio*ysigma 
  92           
  93          return gaussian(self.pattern_x,self.pattern_y,xsigma,ysigma) 
  94   
  95   
  96 -class ExponentialDecay(PatternGenerator): 
  97      """ 
  98      2D Exponential pattern generator. 
  99   
 100      Exponential decay based on distance from a central peak, 
 101      i.e. exp(-d), where d is the distance from the center (assuming 
 102      size=1.0 and aspect_ratio==1.0).  More generally, the size and 
 103      aspect ratio determine the scaling of x and y dimensions: 
 104   
 105        yscale=size/2 
 106        xscale=yscale*aspect_ratio 
 107   
 108      The exponential is then computed for the given (x,y) values as:: 
 109       
 110        exp(-sqrt((x/xscale)^2 - (y/yscale)^2)) 
 111      """ 
 112       
 113      aspect_ratio = param.Number(default=1/0.31,bounds=(0.0,None),softbounds=(0.0,2.0), 
 114          precedence=0.31,doc="""Ratio of the width to the height.""") 
 115       
 116      size = param.Number(default=0.155,doc=""" 
 117          Overall scaling of the x and y dimensions.""") 
 118   
 119 -    def function(self,p): 
 120          yscale = p.size/2.0 
 121          xscale = p.aspect_ratio*yscale 
 122           
 123          return exponential(self.pattern_x,self.pattern_y,xscale,yscale) 
 124   
 125   
 126 -class SineGrating(PatternGenerator): 
 127      """2D sine grating pattern generator.""" 
 128       
 129      frequency = param.Number(default=2.4,bounds=(0.0,None),softbounds=(0.0,10.0), 
 130                         precedence=0.50, doc="Frequency of the sine grating.") 
 131       
 132      phase     = param.Number(default=0.0,bounds=(0.0,None),softbounds=(0.0,2*pi), 
 133                         precedence=0.51,doc="Phase of the sine grating.") 
 134   
 135 -    def function(self,p): 
 136          """Return a sine grating pattern (two-dimensional sine wave).""" 
 137          return 0.5 + 0.5*sin(p.frequency*2*pi*self.pattern_y + p.phase)         
 138   
 139   
 140   
 141 -class Gabor(PatternGenerator): 
 142      """2D Gabor pattern generator.""" 
 143       
 144      frequency = param.Number(default=2.4,bounds=(0.0,None),softbounds=(0.0,10.0), 
 145          precedence=0.50,doc="Frequency of the sine grating component.") 
 146       
 147      phase = param.Number(default=0.0,bounds=(0.0,None),softbounds=(0.0,2*pi), 
 148          precedence=0.51,doc="Phase of the sine grating component.") 
 149       
 150      aspect_ratio = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,2.0), 
 151          precedence=0.31,doc= 
 152          """ 
 153          Ratio of pattern width to height. 
 154          The width of the Gaussian component is size*aspect_ratio (see Gaussian). 
 155          """) 
 156       
 157      size = param.Number(default=0.25,doc=""" 
 158          Determines the height of the Gaussian component (see Gaussian).""") 
 159   
 160 -    def function(self,p): 
 161          height = p.size/2.0 
 162          width = p.aspect_ratio*height 
 163           
 164          return gabor(self.pattern_x,self.pattern_y,width,height, 
 165                       p.frequency,p.phase) 
 166   
 167   
 168 -class Line(PatternGenerator): 
 169      """2D line pattern generator.""" 
 170   
 171      thickness   = param.Number(default=0.006,bounds=(0.0,None),softbounds=(0.0,1.0), 
 172                           precedence=0.60, 
 173                           doc="Thickness (width) of the solid central part of the line.") 
 174      smoothing = param.Number(default=0.05,bounds=(0.0,None),softbounds=(0.0,0.5), 
 175                         precedence=0.61, 
 176                         doc="Width of the Gaussian fall-off.") 
 177   
 178 -    def function(self,p): 
 179          return line(self.pattern_y,p.thickness,p.smoothing) 
 180   
 181   
 182 -class Disk(PatternGenerator): 
 183      """ 
 184      2D disk pattern generator. 
 185   
 186      An elliptical disk can be obtained by adjusting the aspect_ratio of a circular 
 187      disk; this transforms a circle into an ellipse by stretching the circle in the 
 188      y (vertical) direction. 
 189   
 190      The Gaussian fall-off at a point P is an approximation for non-circular disks, 
 191      since the point on the ellipse closest to P is taken to be the same point as 
 192      the point on the circle before stretching that was closest to P. 
 193      """ 
 194   
 195      aspect_ratio  = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,2.0), 
 196          precedence=0.31,doc= 
 197          "Ratio of width to height; size*aspect_ratio gives the width of the disk.") 
 198   
 199      size  = param.Number(default=0.5,doc="Top to bottom height of the disk") 
 200       
 201      smoothing = param.Number(default=0.1,bounds=(0.0,None),softbounds=(0.0,0.5), 
 202                         precedence=0.61,doc="Width of the Gaussian fall-off") 
 203       
 204 -    def function(self,p): 
 205          height = p.size 
 206   
 207          if p.aspect_ratio==0.0: 
 208              return self.pattern_x*0.0 
 209   
 210          return disk(self.pattern_x/p.aspect_ratio,self.pattern_y,height, 
 211                      p.smoothing) 
 212   
 213   
 214 -class Ring(PatternGenerator): 
 215      """ 
 216      2D ring pattern generator. 
 217   
 218      See the Disk class for a note about the Gaussian fall-off. 
 219      """ 
 220   
 221      thickness = param.Number(default=0.015,bounds=(0.0,None),softbounds=(0.0,0.5), 
 222          precedence=0.60,doc="Thickness (line width) of the ring.") 
 223       
 224      smoothing = param.Number(default=0.1,bounds=(0.0,None),softbounds=(0.0,0.5), 
 225          precedence=0.61,doc="Width of the Gaussian fall-off inside and outside the ring.") 
 226       
 227      aspect_ratio = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,2.0), 
 228          precedence=0.31,doc= 
 229          "Ratio of width to height; size*aspect_ratio gives the overall width.") 
 230       
 231      size = param.Number(default=0.5) 
 232   
 233 -    def function(self,p): 
 234          height = p.size 
 235          if p.aspect_ratio==0.0: 
 236              return self.pattern_x*0.0 
 237   
 238          return ring(self.pattern_x/p.aspect_ratio,self.pattern_y,height, 
 239                      p.thickness,p.smoothing) 
 240       
 241   
 242 -class OrientationContrast(SineGrating): 
 243      """ 
 244      Circular pattern for testing responses to differences in contrast. 
 245   
 246      The pattern contains a sine grating ring surrounding a sine grating disk, each 
 247      with parameters (orientation, size, scale and offset) that can be 
 248      changed independently. 
 249      """ 
 250    
 251      orientationcenter   = param.Number(default=0.0,bounds=(0.0,2*pi), doc="Orientation of the center grating.") 
 252      orientationsurround = param.Number(default=0.0,bounds=(0.0,2*pi), doc="Orientation of the surround grating.") 
 253      sizecenter     = param.Number(default=0.5,bounds=(0.0,None),softbounds=(0.0,10.0), doc="Size of the center grating.") 
 254      sizesurround   = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,10.0), doc="Size of the surround grating.") 
 255      scalecenter    = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,10.0), doc="Scale of the center grating.") 
 256      scalesurround  = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,10.0), doc="Scale of the surround grating.") 
 257      offsetcenter   = param.Number(default=0.0,bounds=(0.0,None),softbounds=(0.0,10.0), doc="Offset of the center grating.") 
 258      offsetsurround = param.Number(default=0.0,bounds=(0.0,None),softbounds=(0.0,10.0), doc="Offset of the surround grating.") 
 259      smoothing      = param.Number(default=0.0,bounds=(0.0,None),softbounds=(0.0,0.5),  doc="Width of the Gaussian fall-off inside and outside the ring.") 
 260      thickness      = param.Number(default=0.015,bounds=(0.0,None),softbounds=(0.0,0.5),doc="Thickness (line width) of the ring.") 
 261      aspect_ratio   = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,2.0),  doc="Ratio of width to height; size*aspect_ratio gives the overall width.") 
 262      size           = param.Number(default=0.5) 
 263   
 264 -    def __call__(self,**params_to_override): 
 265          p = ParamOverrides(self,params_to_override) 
 266          input_1=SineGrating(mask_shape=Disk(smoothing=0,size=1.0),phase=p.phase, frequency=p.frequency, 
 267                              orientation=p.orientationcenter, 
 268                              scale=p.scalecenter, offset=p.offsetcenter, 
 269                              x=p.x, y=p.y,size=p.sizecenter) 
 270          input_2=SineGrating(mask_shape=Ring(thickness=p.thickness,smoothing=0,size=1.0),phase=p.phase, frequency=p.frequency, 
 271                              orientation=p.orientationsurround, scale=p.scalesurround, offset=p.offsetsurround, 
 272                              x=p.x, y=p.y, size=p.sizesurround) 
 273   
 274          patterns = [input_1(xdensity=p.xdensity,ydensity=p.ydensity,bounds=p.bounds), 
 275                      input_2(xdensity=p.xdensity,ydensity=p.ydensity,bounds=p.bounds)] 
 276                         
 277          image_array = numpy.add.reduce(patterns) 
 278          return image_array 
 279   
 280   
 281   
 282 -class RawRectangle(PatternGenerator): 
 283      """ 
 284      2D rectangle pattern generator with no smoothing, for use when drawing 
 285      patterns pixel by pixel. 
 286      """ 
 287       
 288      aspect_ratio   = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,2.0), 
 289          precedence=0.31,doc= 
 290          "Ratio of width to height; size*aspect_ratio gives the width of the rectangle.") 
 291       
 292      size  = param.Number(default=0.5,doc="Height of the rectangle.") 
 293   
 294 -    def function(self,p): 
 295          height = p.size 
 296          width = p.aspect_ratio*height 
 297          return bitwise_and(abs(self.pattern_x)<=width/2.0, 
 298                             abs(self.pattern_y)<=height/2.0) 
 299   
 300   
 301   
 302 -class Rectangle(PatternGenerator): 
 303      """2D rectangle pattern, with Gaussian smoothing around the edges.""" 
 304   
 305      aspect_ratio = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,6.0), 
 306          precedence=0.31,doc= 
 307          "Ratio of width to height; size*aspect_ratio gives the width of the rectangle.") 
 308       
 309      size = param.Number(default=0.5,doc="Height of the rectangle.") 
 310   
 311      smoothing = param.Number(default=0.05,bounds=(0.0,None),softbounds=(0.0,0.5), 
 312          precedence=0.61,doc="Width of the Gaussian fall-off outside the rectangle.") 
 313   
 314 -    def function(self,p): 
 315          height=p.size 
 316          width=p.aspect_ratio*height 
 317   
 318          return smooth_rectangle(self.pattern_x, self.pattern_y,  
 319                                  width, height, p.smoothing, p.smoothing) 
 320   
 321   
 322   
 323 -class Arc(PatternGenerator): 
 324      """ 
 325      2D arc pattern generator. 
 326       
 327      Draws an arc (partial ring) of the specified size (radius*2), 
 328      starting at radian 0.0 and ending at arc_length.  The orientation 
 329      can be changed to choose other start locations.  The pattern is 
 330      centered at the center of the ring. 
 331   
 332      See the Disk class for a note about the Gaussian fall-off. 
 333      """ 
 334   
 335      aspect_ratio = param.Number(default=1.0,bounds=(0.0,None),softbounds=(0.0,6.0), 
 336          precedence=0.31,doc=""" 
 337          Ratio of width to height; size*aspect_ratio gives the overall width.""") 
 338   
 339      thickness = param.Number(default=0.015,bounds=(0.0,None),softbounds=(0.0,0.5), 
 340          precedence=0.60,doc="Thickness (line width) of the ring.") 
 341       
 342      smoothing = param.Number(default=0.05,bounds=(0.0,None),softbounds=(0.0,0.5), 
 343          precedence=0.61,doc="Width of the Gaussian fall-off inside and outside the ring.") 
 344   
 345      arc_length = param.Number(default=pi,bounds=(0.0,None),softbounds=(0.0,2.0*pi), 
 346                                inclusive_bounds=(True,False),precedence=0.62, doc=""" 
 347          Length of the arc, in radians, starting from orientation 0.0.""") 
 348   
 349      size = param.Number(default=0.5) 
 350   
 351 -    def function(self,p): 
 352          if p.aspect_ratio==0.0: 
 353              return self.pattern_x*0.0 
 354   
 355          return arc_by_radian(self.pattern_x/p.aspect_ratio, self.pattern_y, p.size, 
 356                               (2*pi-p.arc_length, 0.0), p.thickness, p.smoothing) 
 357   
 358   
 359 -class Curve(Arc): 
 360      """ 
 361      2D curve pattern generator. 
 362   
 363      Based on Arc, but centered on a tangent point midway through the 
 364      arc, rather than at the center of a ring, and with curvature 
 365      controlled directly rather than through the overall size of the 
 366      pattern. 
 367   
 368      Depending on the size_type, the size parameter can control either 
 369      the width of the pattern, keeping this constant regardless of 
 370      curvature, or the length of the curve, keeping that constant 
 371      instead (as for a long thin object being bent). 
 372   
 373      Specifically, for size_type=='constant_length', the curvature 
 374      parameter determines the ratio of height to width of the arc, with 
 375      positive curvature for concave shape and negative for convex. The 
 376      size parameter determines the width of the curve. 
 377   
 378      For size_type=='constant_width', the curvature parameter 
 379      determines the portion of curve radian to 2pi, and the curve 
 380      radius is changed accordingly following the formula:: 
 381   
 382        size=2pi*radius*curvature 
 383   
 384      Thus, the size parameter determines the total length of the 
 385      curve. Positive curvature stands for concave shape, and negative 
 386      for convex. 
 387   
 388      See the Disk class for a note about the Gaussian fall-off. 
 389      """ 
 390   
 391      # Hide unused parameters 
 392      arc_length = param.Number(precedence=-1.0) 
 393      aspect_ratio = param.Number(default=1.0, precedence=-1.0) 
 394   
 395      size_type = param.ObjectSelector(default='constant_length', 
 396          objects=['constant_length','constant_width'],precedence=0.61,doc=""" 
 397          For a given size, whether to draw a curve with that total length,  
 398          or with that width, keeping it constant as curvature is varied.""") 
 399   
 400      curvature = param.Number(default=0.5, bounds=(-0.5, 0.5), precedence=0.62, doc=""" 
 401          Ratio of height to width of the arc, with positive value giving 
 402          a concave shape and negative value giving convex.""") 
 403   
 404 -    def function(self,p): 
 405          return arc_by_center(self.pattern_x/p.aspect_ratio,self.pattern_y,  
 406                               (p.size,p.size*p.curvature), 
 407                               (p.size_type=='constant_length'),  
 408                               p.thickness, p.smoothing) 
 409   
 410   
 411   
 412  #JABALERT: Can't this be replaced with a Composite? 
 413 -class TwoRectangles(Rectangle): 
 414      """Two 2D rectangle pattern generator.""" 
 415   
 416      x1 = param.Number(default=-0.15,bounds=(-1.0,1.0),softbounds=(-0.5,0.5), 
 417                  doc="X center of rectangle 1.") 
 418       
 419      y1 = param.Number(default=-0.15,bounds=(-1.0,1.0),softbounds=(-0.5,0.5), 
 420                  doc="Y center of rectangle 1.") 
 421       
 422      x2 = param.Number(default=0.15,bounds=(-1.0,1.0),softbounds=(-0.5,0.5), 
 423                  doc="X center of rectangle 2.") 
 424       
 425      y2 = param.Number(default=0.15,bounds=(-1.0,1.0),softbounds=(-0.5,0.5), 
 426                  doc="Y center of rectangle 2.") 
 427   
 428      # YC: Maybe this can be implemented much more cleanly by calling 
 429      # the parent's function() twice, but it's hard to see how to  
 430      # set the (x,y) offset for the parent. 
 431 -    def function(self,p): 
 432          height = p.size 
 433          width = p.aspect_ratio*height 
 434   
 435          return bitwise_or( 
 436                 bitwise_and(bitwise_and( 
 437                          (self.pattern_x-p.x1)<=p.x1+width/4.0, 
 438                          (self.pattern_x-p.x1)>=p.x1-width/4.0), 
 439                        bitwise_and( 
 440                          (self.pattern_y-p.y1)<=p.y1+height/4.0, 
 441                          (self.pattern_y-p.y1)>=p.y1-height/4.0)), 
 442                 bitwise_and(bitwise_and( 
 443                          (self.pattern_x-p.x2)<=p.x2+width/4.0, 
 444                          (self.pattern_x-p.x2)>=p.x2-width/4.0), 
 445                        bitwise_and( 
 446                          (self.pattern_y-p.y2)<=p.y2+height/4.0, 
 447                          (self.pattern_y-p.y2)>=p.y2-height/4.0))) 
 448   
 449   
 450 -class SquareGrating(PatternGenerator): 
 451      """2D squarewave grating pattern generator.""" 
 452       
 453      frequency = param.Number(default=2.4,bounds=(0.0,None),softbounds=(0.0,10.0), 
 454          precedence=0.50,doc="Frequency of the square grating.") 
 455       
 456      phase     = param.Number(default=0.0,bounds=(0.0,None),softbounds=(0.0,2*pi), 
 457          precedence=0.51,doc="Phase of the square grating.") 
 458   
 459      # We will probably want to add anti-aliasing to this, 
 460      # and there might be an easier way to do it than by 
 461      # cropping a sine grating. 
 462   
 463 -    def function(self,p): 
 464          """ 
 465          Return a square-wave grating (alternating black and white bars). 
 466          """ 
 467          return around(0.5 + 0.5*sin(p.frequency*2*pi*self.pattern_y + p.phase)) 
 468   
 469   
 470  # CB: I removed motion_sign from this class because I think it is 
 471  # unnecessary. But maybe I misunderstood the original author's 
 472  # intention? 
 473  # 
 474  # In any case, the original implementation was incorrect - it was not 
 475  # possible to get some motion directions (directions in one whole 
 476  # quadrant were missed out). 
 477  # 
 478  # Note that to get a 2pi range of directions, one must use a 2pi range 
 479  # of orientations (there are two directions for any given 
 480  # orientation).  Alternatively, we could generate a random sign, and 
 481  # use an orientation restricted to a pi range. 
 482   
 483 -class Sweeper(PatternGenerator): 
 484      """ 
 485      PatternGenerator that sweeps a supplied PatternGenerator in a direction 
 486      perpendicular to its orientation. 
 487      """ 
 488   
 489      generator = param.Parameter(default=Gaussian(),precedence=0.97, doc="Pattern to sweep.") 
 490   
 491      speed = param.Number(default=0.25,bounds=(0.0,None),doc=""" 
 492          Sweep speed: number of sheet coordinate units per unit time.""") 
 493   
 494      step = param.Number(default=1,doc=""" 
 495          Number of steps at the given speed to move in the sweep direction. 
 496          The distance moved is speed*step.""") 
 497   
 498      # Provide access to value needed for measuring maps 
 499 -    def __get_phase(self): return self.generator.phase 
 500 -    def __set_phase(self,new_val): self.generator.phase = new_val 
 501      phase = property(__get_phase,__set_phase) 
 502   
 503 -    def function(self,p): 
 504          """Selects and returns one of the patterns in the list.""" 
 505          pg = p.generator 
 506          motion_orientation=p.orientation+pi/2.0 
 507   
 508          new_x = p.x+p.size*pg.x 
 509          new_y = p.y+p.size*pg.y 
 510           
 511          image_array = pg(xdensity=p.xdensity,ydensity=p.ydensity,bounds=p.bounds, 
 512                           x=new_x + p.speed*p.step*cos(motion_orientation), 
 513                           y=new_y + p.speed*p.step*sin(motion_orientation), 
 514                           orientation=p.orientation, 
 515                           scale=pg.scale*p.scale,offset=pg.offset+p.offset) 
 516           
 517          return image_array 
 518   
 519   
 520 -class Composite(PatternGenerator): 
 521      """ 
 522      PatternGenerator that accepts a list of other PatternGenerators. 
 523      To create a new pattern, asks each of the PatternGenerators in the 
 524      list to create a pattern, then it combines the patterns to create a  
 525      single pattern that it returns. 
 526      """ 
 527   
 528      # The Accum_Replace operator from LISSOM is not yet supported, 
 529      # but it should be added once PatternGenerator bounding boxes 
 530      # are respected and/or GenericImage patterns support transparency. 
 531      operator = param.Parameter(numpy.maximum,precedence=0.98,doc=""" 
 532          Binary Numpy function used to combine the individual patterns. 
 533   
 534          Any binary Numpy array "ufunc" returning the same 
 535          type of array as the operands and supporting the reduce 
 536          operator is allowed here.  Supported ufuncs include:: 
 537   
 538            add 
 539            subtract 
 540            multiply 
 541            divide 
 542            maximum 
 543            minimum 
 544            remainder 
 545            power 
 546            logical_and 
 547            logical_or 
 548            logical_xor 
 549   
 550          The most useful ones are probably add and maximum, but there 
 551          are uses for at least some of the others as well (e.g. to 
 552          remove pieces of other patterns). 
 553   
 554          You can also write your own operators, by making a class that 
 555          has a static method named "reduce" that returns an array of the 
 556          same size and type as the arrays in the list.  For example:: 
 557           
 558            class return_first(object): 
 559                @staticmethod 
 560                def reduce(x): 
 561                    return x[0] 
 562                 
 563          """) 
 564       
 565      generators = param.List(default=[Constant(scale=0.0)],precedence=0.97, 
 566          class_=PatternGenerator,doc=""" 
 567          List of patterns to use in the composite pattern.  The default is 
 568          a blank pattern, and should thus be overridden for any useful work.""") 
 569   
 570      size  = param.Number(default=1.0,doc="Scaling factor applied to all sub-patterns.") 
 571   
 572       
 573 -    def _advance_pattern_generators(self,p): 
 574          """ 
 575          Subclasses can override this method to provide constraints on 
 576          the values of generators' parameters and/or eliminate 
 577          generators from this list if necessary. 
 578          """ 
 579          return p.generators 
 580   
 581           
 582      # JABALERT: To support large numbers of patterns on a large input region, 
 583      # should be changed to evaluate each pattern in a small box, and then 
 584      # combine them at the full Composite Bounding box size. 
 585 -    def function(self,p): 
 586          """Constructs combined pattern out of the individual ones.""" 
 587          generators = self._advance_pattern_generators(p) 
 588   
 589          assert hasattr(p.operator,'reduce'),repr(p.operator)+" does not support 'reduce'." 
 590   
 591          # CEBALERT: mask gets applied by all PGs including the Composite itself 
 592          # (leads to redundant calculations in current lissom_oo_or usage, but 
 593          # will lead to problems/limitations in the future). 
 594          patterns = [pg(xdensity=p.xdensity,ydensity=p.ydensity, 
 595                         bounds=p.bounds,mask=p.mask, 
 596                         x=p.x+p.size*(pg.x*cos(p.orientation)- pg.y*sin(p.orientation)), 
 597                         y=p.y+p.size*(pg.x*sin(p.orientation)+ pg.y*cos(p.orientation)), 
 598                         orientation=pg.orientation+p.orientation, 
 599                         size=pg.size*p.size) 
 600                      for pg in generators] 
 601          image_array = p.operator.reduce(patterns) 
 602          return image_array 
 603   
 604   
 605   
 606 -class SeparatedComposite(Composite): 
 607      """ 
 608      Generalized version of the Composite PatternGenerator that enforces spacing constraints 
 609      between pattern centers. 
 610   
 611      Currently supports minimum spacing, but can be generalized to 
 612      support maximum spacing also (and both at once). 
 613      """ 
 614   
 615      min_separation = param.Number(default=0.0, bounds = (0,None), 
 616                              softbounds = (0.0,1.0), doc=""" 
 617          Minimum distance to enforce between all pairs of pattern centers. 
 618   
 619          Useful for ensuring that multiple randomly generated patterns 
 620          do not overlap spatially.  Note that as this this value is 
 621          increased relative to the area in which locations are chosen, 
 622          the likelihood of a pattern appearing near the center of the 
 623          area will decrease.  As this value approaches the available 
 624          area, the corners become far more likely to be chosen, due to 
 625          the distances being greater along the diagonals. 
 626          """) 
 627          ### JABNOTE: Should provide a mechanism for collecting and 
 628          ### plotting the training pattern center distribution, so that 
 629          ### such issues can be checked. 
 630   
 631      max_trials = param.Integer(default = 50, bounds = (0,None), 
 632                           softbounds = (0,100), precedence=-1, doc=""" 
 633          Number of times to try for a new pattern location that meets the criteria. 
 634           
 635          This is an essentially arbitrary timeout value that helps 
 636          prevent an endless loop in case the requirements cannot be 
 637          met.""") 
 638   
 639           
 640 -    def __distance_valid(self, g0, g1, p): 
 641          """ 
 642          Returns true if the distance between the (x,y) locations of two generators 
 643          g0 and g1 is greater than a minimum separation. 
 644   
 645          Can be extended easily to support other criteria. 
 646          """ 
 647          dist = sqrt((g1.x - g0.x) ** 2 + 
 648                      (g1.y - g0.y) ** 2) 
 649          return dist >= p.min_separation 
 650   
 651   
 652 -    def _advance_pattern_generators(self,p): 
 653          """ 
 654          Advance the parameters for each generator for this presentation. 
 655   
 656          Picks a position for each generator that is accepted by __distance_valid 
 657          for all combinations.  Returns a new list of the generators, with 
 658          some potentially omitted due to failure to meet the constraints. 
 659          """ 
 660           
 661          valid_generators = [] 
 662          for g in p.generators: 
 663               
 664              for trial in xrange(self.max_trials): 
 665                  # Generate a new position and add generator if it's ok 
 666                   
 667                  if alltrue([self.__distance_valid(g,v,p) for v in valid_generators]): 
 668                      valid_generators.append(g) 
 669                      break 
 670                   
 671                  vals = (g.force_new_dynamic_value('x'), g.force_new_dynamic_value('y')) 
 672                   
 673              else: 
 674                  self.warning("Unable to place pattern %s subject to given constraints" % 
 675                               g.name) 
 676   
 677          return valid_generators 
 678   
 679   
 680   
 681 -class Selector(PatternGenerator): 
 682      """ 
 683      PatternGenerator that selects from a list of other PatternGenerators. 
 684      """ 
 685   
 686      generators = param.List(precedence=0.97,class_=PatternGenerator,bounds=(1,None), 
 687          default=[Disk(x=-0.3,aspect_ratio=0.5), Rectangle(x=0.3,aspect_ratio=0.5)], 
 688          doc="List of patterns from which to select.") 
 689   
 690      size = param.Number(default=1.0,doc="Scaling factor applied to all sub-patterns.") 
 691   
 692      # CB: needs to have time_fn=None 
 693      index = param.Number(default=numbergen.UniformRandom(lbound=0,ubound=1.0,seed=76), 
 694          bounds=(-1.0,1.0),precedence=0.20,doc=""" 
 695          Index into the list of pattern generators, on a scale from 0 
 696          (start of the list) to 1.0 (end of the list).  Typically a 
 697          random value or other number generator, to allow a different item 
 698          to be selected each time.""") 
 699   
 700   
 701 -    def function(self,p): 
 702          """Selects and returns one of the patterns in the list.""" 
 703          int_index=int(len(p.generators)*wrap(0,1.0,p.index)) 
 704          pg=p.generators[int_index] 
 705   
 706          image_array = pg(xdensity=p.xdensity,ydensity=p.ydensity,bounds=p.bounds, 
 707                           x=p.x+p.size*(pg.x*cos(p.orientation)-pg.y*sin(p.orientation)), 
 708                           y=p.y+p.size*(pg.x*sin(p.orientation)+pg.y*cos(p.orientation)), 
 709                           orientation=pg.orientation+p.orientation,size=pg.size*p.size, 
 710                           scale=pg.scale*p.scale,offset=pg.offset+p.offset) 
 711                          
 712          return image_array 
 713   
 714 -    def get_current_generator(self): 
 715          """Return the current generator (as specified by self.index).""" 
 716          int_index=int(len(self.generators)*wrap(0,1.0,self.inspect_value('index'))) 
 717          return self.generators[int_index] 
 718   
 719   
 720   
 721   
 722  ### JABALERT: This class should be eliminated if at all possible; it 
 723  ### is just a specialized version of Composite, and should be 
 724  ### implementable directly using what is already in Composite.     
 725 -class GaussiansCorner(PatternGenerator): 
 726      """ 
 727      Two Gaussian pattern generators with a variable intersection point,  
 728      appearing as a corner or cross. 
 729      """ 
 730       
 731      x = param.Number(default=-0.15,bounds=(-1.0,1.0),softbounds=(-0.5,0.5), 
 732                  doc="X center of the corner") 
 733       
 734      y = param.Number(default=-0.15,bounds=(-1.0,1.0),softbounds=(-0.5,0.5), 
 735                  doc="Y center of the corner") 
 736                   
 737      size = param.Number(default=0.5,bounds=(0,None), softbounds=(0.1,1), 
 738                  doc="The size of the corner") 
 739       
 740      aspect_ratio = param.Number(default=1/0.31, bounds=(0,None), softbounds=(1,10), 
 741                  doc="Ratio of the width to the height for both Gaussians") 
 742                   
 743      angle = param.Number(default=0.5*pi,bounds=(0,pi), softbounds=(0.01*pi,0.99*pi), 
 744                  doc="The angle of the corner") 
 745                   
 746      cross = param.Number(default=0.4, bounds=(0,1), softbounds=(0,1), 
 747                  doc="Where the two Gaussians cross, as a fraction of their half length") 
 748       
 749   
 750 -    def __call__(self,**params_to_override): 
 751          p = ParamOverrides(self,params_to_override) 
 752           
 753          g_1 = Gaussian() 
 754          g_2 = Gaussian() 
 755           
 756          x_1 = g_1(orientation = p.orientation, bounds = p.bounds, xdensity = p.xdensity, 
 757                              ydensity = p.ydensity, offset = p.offset, size = p.size, 
 758                              aspect_ratio = p.aspect_ratio, 
 759                              x = p.x + 0.7 * cos(p.orientation) * p.cross * p.size * p.aspect_ratio, 
 760                              y = p.y + 0.7 * sin(p.orientation) * p.cross * p.size * p.aspect_ratio) 
 761          x_2 = g_2(orientation = p.orientation+p.angle, bounds = p.bounds, xdensity = p.xdensity, 
 762                              ydensity = p.ydensity, offset = p.offset, size = p.size, 
 763                              aspect_ratio = p.aspect_ratio, 
 764                              x = p.x + 0.7 * cos(p.orientation+p.angle) * p.cross * p.size * p.aspect_ratio, 
 765                              y = p.y + 0.7 * sin(p.orientation+p.angle) * p.cross * p.size * p.aspect_ratio) 
 766           
 767          return numpy.maximum( x_1, x_2 ) 
 768   
 769   
 770   
 771 -class Translator(PatternGenerator): 
 772      """ 
 773      PatternGenerator that translates another PatternGenerator over 
 774      time. 
 775   
 776      This PatternGenerator will create a series of episodes, where in 
 777      each episode the underlying generator is moved in a fixed 
 778      direction at a fixed speed.  To begin an episode, the Translator's 
 779      x, y, and direction are evaluated (e.g. from random 
 780      distributions), and the underlying generator is then drawn at 
 781      those values plus changes over time that are determined by the 
 782      speed.  The orientation of the underlying generator should be set 
 783      to 0 to get motion perpendicular to the generator's orientation 
 784      (which is typical). 
 785   
 786      Note that at present the parameter values for x, y, and direction 
 787      cannot be passed in when the instance is called; only the values 
 788      set on the instance are used. 
 789      """ 
 790      generator = param.ClassSelector(default=Gaussian(), 
 791          class_=PatternGenerator,doc="""Pattern to be translated.""") 
 792   
 793      direction = param.Number(default=0.0,softbounds=(-pi,pi),doc=""" 
 794          The direction in which the pattern should move, in radians.""") 
 795       
 796      speed = param.Number(default=0.01,bounds=(0.0,None),doc=""" 
 797          The speed with which the pattern should move, 
 798          in sheet coordinates per simulation time unit.""") 
 799   
 800      reset_period = param.Number(default=1,bounds=(0.0,None),doc=""" 
 801          Period between generating each new translation episode.""") 
 802   
 803      episode_interval = param.Number(default=0,doc=""" 
 804          Interval between successive translation episodes. 
 805           
 806          If nonzero, the episode_separator pattern is presented for 
 807          this amount of simulation time after each episode, e.g. to 
 808          allow processing of the previous episode to complete.""") 
 809   
 810      episode_separator = param.ClassSelector(default=Constant(scale=0.0), 
 811           class_=PatternGenerator,doc=""" 
 812           Pattern to display during the episode_interval, if any. 
 813           The default is a blank pattern.""") 
 814                                                                                 
 815   
 816 -    def _advance_params(self): 
 817          """ 
 818          Explicitly generate new values for these parameters only 
 819          when appropriate. 
 820          """ 
 821          for param in ['x','y','direction']: 
 822              self.force_new_dynamic_value(param) 
 823          self.last_time = topo.sim.time() 
 824   
 825          
 826 -    def __init__(self,**params): 
 827          super(Translator,self).__init__(**params) 
 828          self._advance_params() 
 829   
 830           
 831 -    def __call__(self,**params_to_override): 
 832          p=ParamOverrides(self,params_to_override) 
 833           
 834          if topo.sim.time() >= self.last_time + p.reset_period: 
 835              ## Returns early if within episode interval 
 836              if topo.sim.time()<self.last_time+p.reset_period+p.episode_interval: 
 837                  return p.episode_separator(xdensity=p.xdensity, 
 838                                             ydensity=p.ydensity, 
 839                                             bounds=p.bounds) 
 840              else: 
 841                  self._advance_params() 
 842   
 843          # JABALERT: Does not allow x, y, or direction to be passed in 
 844          # to the call; fixing this would require implementing 
 845          # inspect_value and force_new_dynamic_value (for 
 846          # use in _advance_params) for ParamOverrides. 
 847          # 
 848          # Access parameter values without giving them new values 
 849          assert ('x' not in params_to_override and 
 850                  'y' not in params_to_override and 
 851                  'direction' not in params_to_override) 
 852          x = self.inspect_value('x') 
 853          y = self.inspect_value('y') 
 854          direction = self.inspect_value('direction') 
 855   
 856          # compute how much time elapsed from the last reset 
 857          # float(t) required because time could be e.g. gmpy.mpq 
 858          t = float(topo.sim.time()-self.last_time) 
 859   
 860          ## CEBALERT: mask gets applied twice, both for the underlying 
 861          ## generator and for this one.  (leads to redundant 
 862          ## calculations in current lissom_oo_or usage, but will lead 
 863          ## to problems/limitations in the future). 
 864          return p.generator( 
 865              xdensity=p.xdensity,ydensity=p.ydensity,bounds=p.bounds, 
 866              x=x+t*cos(direction)*p.speed+p.generator.x, 
 867              y=y+t*sin(direction)*p.speed+p.generator.y, 
 868              orientation=(direction-pi/2)+p.generator.orientation) 
 869   
 870   
 871 -class DifferenceOfGaussians(PatternGenerator): 
 872      """ 
 873      Two-dimensional difference of gaussians pattern. 
 874      """ 
 875   
 876      positive_size = param.Number(default=0.5, bounds=(0.0,None), softbounds=(0.0,5.0), 
 877          precedence=(1), doc="""size parameter for the center Gaussian.""") 
 878       
 879      positive_aspect_ratio = param.Number(default=2.0, bounds=(0.0,None), softbounds=(0.0,5.0), 
 880          precedence=(2), doc="""aspect_ratio parameter for the center Gaussian.""") 
 881       
 882      positive_x = param.Number(default=0.0, bounds=(None,None), softbounds=(-2.0,2.0), 
 883          precedence=(3), doc="""x position for the central peak of the positive gaussian.""") 
 884       
 885      positive_y = param.Number(default=0.0, bounds=(None,None), softbounds=(-2.0,2.0), 
 886          precedence=(4), doc="""y position for the central peak of the positive gaussian.""") 
 887   
 888      negative_size = param.Number(default=1.0, bounds=(0.0,None), softbounds=(0.0,5.0), 
 889          precedence=(5), doc="""size parameter for the surround Gaussian.""") 
 890       
 891      negative_aspect_ratio = param.Number(default=2.0, bounds=(0.0,None), softbounds=(0.0,5.0), 
 892          precedence=(6), doc="""aspect_ratio parameter for the surround Gaussian.""") 
 893       
 894      negative_x = param.Number(default=0.0, bounds=(None,None), softbounds=(-2.0,2.0), 
 895          precedence=(7), doc="""x position for the central peak of the negative gaussian.""") 
 896       
 897      negative_y = param.Number(default=0.0, bounds=(None,None), softbounds=(-2.0,2.0), 
 898          precedence=(8), doc="""y position for the central peak of the negative gaussian.""") 
 899       
 900 -    def function(self, p): 
 901          center   = Gaussian(size=p.positive_size*p.size, aspect_ratio=p.positive_aspect_ratio, 
 902                              orientation=p.orientation, x=p.positive_x+p.x, y=p.positive_y+p.y, 
 903                              output_fns=[topo.transferfn.DivisiveNormalizeL1()]) 
 904                                       
 905          surround = Gaussian(size=p.negative_size*p.size, aspect_ratio=p.negative_aspect_ratio, 
 906                              orientation=p.orientation, x=p.negative_x+p.x, y=p.negative_y+p.y, 
 907                              output_fns=[topo.transferfn.DivisiveNormalizeL1()]) 
 908           
 909          return Composite(generators=[center,surround], operator=numpy.subtract, 
 910                           xdensity=p.xdensity, ydensity=p.ydensity, bounds=p.bounds)() 
 911                           
 912                           
 913 -class Sigmoid(PatternGenerator): 
 914      """ 
 915      Two-dimensional sigmoid pattern, dividing the plane into positive 
 916      and negative halves with a smoothly sloping transition between them. 
 917      """ 
 918       
 919      slope = param.Number(default=10.0, bounds=(None,None), softbounds=(-100.0,100.0),doc=""" 
 920          Multiplicative parameter controlling the smoothness of the transition  
 921          between the two regions; high values give a sharp transition.""") 
 922   
 923 -    def function(self, p): 
 924          return sigmoid(self.pattern_y, p.slope) 
 925            
 926   
 927 -class SigmoidedDoG(PatternGenerator): 
 928      """ 
 929      Sigmoid multiplicatively combined with a difference of Gaussians, 
 930      such that one part of the plane can be the mirror image of the other. 
 931      """ 
 932           
 933      center_size = param.Number(default=0.5, bounds=(0.0,None), softbounds=(0.0,5.0), 
 934          precedence=(1), doc="""size parameter for the center Gaussian.""") 
 935       
 936      center_aspect_ratio = param.Number(default=2.0, bounds=(0.0,None), softbounds=(0.0,5.0), 
 937          precedence=(2), doc="""aspect_ratio parameter for the center Gaussian.""") 
 938       
 939      surround_size = param.Number(default=1.0, bounds=(0.0,None), softbounds=(0.0,5.0), 
 940          precedence=(3), doc="""size parameter for the surround Gaussian.""") 
 941       
 942      surround_aspect_ratio = param.Number(default=1.0, bounds=(0.0,None), softbounds=(0.0,5.0), 
 943          precedence=(4), doc="""aspect_ratio parameter for the surround Gaussian.""") 
 944       
 945      sigmoid_slope = param.Number(default=10.0, bounds=(None,None), softbounds=(-100.0,100.0), 
 946          precedence=(5), doc="""slope parameter for the Sigmoid.""") 
 947       
 948      sigmoid_x = param.Number(default=0.0, bounds=(None,None), softbounds=(-1.0,1.0), 
 949          precedence=(6), doc="""x parameter for the Sigmoid.""") 
 950   
 951      sigmoid_y = param.Number(default=0.0, bounds=(None,None), softbounds=(-1.0,1.0), 
 952          precedence=(7), doc="""y parameter for the Sigmoid.""") 
 953   
 954                                                                                          
 955 -    def function(self, p): 
 956          dog = DifferenceOfGaussians(positive_size=p.center_size*p.size, positive_aspect_ratio=p.center_aspect_ratio, 
 957                                      negative_size=p.surround_size*p.size, negative_aspect_ratio=p.surround_aspect_ratio, 
 958                                      positive_x=p.x, positive_y=p.y, negative_x=p.x, negative_y=p.y) 
 959   
 960          sigmoid = Sigmoid(slope=p.sigmoid_slope, orientation=p.orientation+pi/2, 
 961                            x=p.sigmoid_x+p.x, y=p.sigmoid_y+p.y) 
 962   
 963          return Composite(generators=[dog, sigmoid], operator=numpy.multiply, 
 964                           xdensity=p.xdensity, ydensity=p.ydensity, bounds=p.bounds)() 
 965      
 966   
 967                                
 968 -def rectangular(signal_size): 
 969      """ 
 970      Generates a Rectangular signal smoothing window, 
 971      """ 
 972      return [1.0]*int(signal_size) 
 973       
 974   
 975       
 976 -class PowerSpectrum(PatternGenerator): 
 977      """ 
 978      Outputs the spectral density of a rolling window of the input 
 979      signal each time it is called. Over time, the results could be 
 980      arranged into a spectrogram, e.g. for an audio signal. 
 981      """ 
 982      # Can be instantiated by hand, but it makes little sense to do so 
 983      # in e.g. a GUI, since it requires a "signal" array parameter to 
 984      # do anything useful.  Alternatively, could provide a useful 
 985      # default for "signal" (and presumably make it a parameter), so 
 986      # that it could be instantiated at least as an example. 
 987      __abstract=True  
 988       
 989      window_increment = param.Number(default=1,constant=True,doc=""" 
 990          The most recent portion of the signal on which to perform the Fourier 
 991          transform, in units of 1/sample_rate, i.e., the length of a 
 992          sliding window on which to operate. 
 993   
 994          Note that the Fourier transform algorithm is most efficient 
 995          for matrix sizes that are powers of 2, or that can be 
 996          decomposed into small prime factors; see numpy.fft.rfft.""" ) 
 997           
 998      window_length = param.Number(default=0.0001,constant=True,doc=""" 
 999          The amount of overlap between each window, in units of 1/sample_rate.""") 
1000   
1001      sample_rate = param.Number(default=44100,constant=True,doc=""" 
1002          Number of samples per second, which defines the range for frequency.""") 
1003   
1004      windowing_function = param.Parameter(default=rectangular,constant=True,doc=""" 
1005          This function is multiplied with the current window, i.e. the 
1006          most recent portion of the waveform interval of a signal, before 
1007          performing the Fourier transform.  It thus shapes the 
1008          interval, which would otherwise always be rectangular. 
1009   
1010          The function chosen here dictates the tradeoff between 
1011          resolving comparable signal strengths with similar 
1012          frequencies, and resolving disparate signal strengths with 
1013          dissimilar frequencies. 
1014   
1015          numpy provides a number of options, e.g. bartlett, blackman, 
1016          hamming, hanning, kaiser; see 
1017          http://docs.scipy.org/doc/numpy/reference/routines.window.html 
1018          You can also supply your own.""") 
1019   
1020      min_frequency = param.Number(default=1,doc=""" 
1021          Smallest frequency for which to return an amplitude.""") 
1022   
1023      max_frequency = param.Number(default=20000,doc=""" 
1024          Largest frequency for which to return an amplitude.""") 
1025   
1026                   
1027 -    def __init__(self, signal, **params): 
1028          super(PowerSpectrum, self).__init__(**params) 
1029          self._initialize_window_parameters(signal, **params) 
1030   
1031      @as_uninitialized 
1032 -    def _initialize_window_parameters(self, signal, **params): 
1033          # BKALERT: how much preprocessing to offer? E.g. Offer to remove DC? Etc. 
1034                   
1035          # For subclasses: to specify the values of parameters on this, the parent class, 
1036          # subclasses might first need access to their own parameter values.  
1037          # Having the window initialization in this separate method allows subclasses 
1038          # to make the usual super.__init__(**params) call. 
1039          for parameter,value in params.items(): 
1040              setattr(self,parameter,value) 
1041           
1042          self.signal = asarray(signal, dtype=float32)         
1043          assert len(self.signal) > 0 
1044           
1045          self._window_start = 0 
1046          self._samples_per_window = int(self.window_length*self.sample_rate) 
1047          self._smoothing_window = self.windowing_function(self._samples_per_window)   
1048   
1049          assert self._samples_per_window > 0 
1050           
1051          # calculate the discrete frequency bins possible for the given sample rate. 
1052          # (discarding the negative frequencies) 
1053          self._all_frequencies = fft.fftfreq(self._samples_per_window, d=1.0/self.sample_rate)[0:self._samples_per_window/2] 
1054          assert self._all_frequencies.min() >= 0 
1055        
1056 -    def _create_spacing(self, mini, maxi):  
1057          """ 
1058          Overload if custom frequency spacing is required. 
1059          """ 
1060          # frequency spacing to use, i.e. mapping of frequencies to sheet rows,  
1061          self._frequency_indices = round(linspace(maxi, mini, num=(maxi-mini), endpoint=True)).astype(int) 
1062                 
1063      # CEBALERT: given all the constant params, could do some caching 
1064 -    def _create_indices(self, p): 
1065          if not self._all_frequencies.min() <= p.min_frequency \ 
1066              or not self._all_frequencies.max() >= p.max_frequency: 
1067              raise ValueError("Specified frequency interval [%s,%s] is unavailable \ 
1068                              (actual interval is [%s,%s]. Adjust sample_rate and/or window_length." 
1069                               %( p.min_frequency, p.max_frequency, \ 
1070                               self._all_frequencies.min(), self._all_frequencies.max() )) 
1071           
1072          # index of first nonzero element whose value is above or equal to the min frequency. 
1073          mini = nonzero(self._all_frequencies >= p.min_frequency)[0][0] 
1074          # index of first nonzero element whose value is above or equal to the max frequency. 
1075          maxi = nonzero(self._all_frequencies <= p.max_frequency)[0][-1] 
1076          self._create_spacing(mini, maxi) 
1077           
1078 -    def _extract_sample_window(self, p): 
1079          """ 
1080          Overload if special behaviour is required when a signal ends. 
1081          """ 
1082          start = self._window_start 
1083          end = start+self._samples_per_window 
1084           
1085          # move window forward for next cycle 
1086          self._window_start += int(self.window_increment * self.sample_rate)  
1087           
1088          if end > self.signal.size: 
1089              raise ValueError("Reached the end of the signal.") 
1090          return self.signal[start:end] 
1091               
1092 -    def _get_amplitudes(self, p): 
1093          """ 
1094          Perform a real Discrete Fourier Transform (DFT; implemented 
1095          using a Fast Fourier Transform algorithm, FFT) of the current 
1096          sample from the signal multiplied by the smoothing window. 
1097   
1098          See numpy.rfft for information about the Fourier transform. 
1099          """ 
1100          # perform a fft on it for amplitudes, discarding the negative frequency values. 
1101          signal_sample = self._extract_sample_window(p) 
1102          assert shape(signal_sample)[0] == shape(self._smoothing_window)[0] 
1103          all_amplitudes = abs(fft.rfft(signal_sample * self._smoothing_window))[0 : len(signal_sample)/2] 
1104                         
1105          # number of discrete units defined by our sheet size, into which we need to  
1106          # partition our frequency index space (and thereby frequency space). 
1107          indices_per_unit = len(self._frequency_indices)/self._sheet_dimensions[0] 
1108           
1109          # list to store summed amplitudes per unit 
1110          amplitudes = [0.0]*self._sheet_dimensions[0] 
1111   
1112          # for each unit in the total number of sheet units 
1113          for unit in range(0, self._sheet_dimensions[0]): 
1114           
1115              # find the largest frequency that falls in it        
1116              frequency_end_index = self._frequency_indices[unit*indices_per_unit]  
1117              # and the smallest frequency that falls in it 
1118              frequency_start_index = self._frequency_indices[(unit*indices_per_unit)+indices_per_unit] 
1119               
1120              # Useful when calibrating: prints out frequencies that fall each unit 
1121              #print (self._all_frequencies[frequency_start_index], self._all_frequencies[frequency_end_index]) 
1122               
1123              # store the sum of the amplitudes for each of those frequencies 
1124              for frequency in range(frequency_start_index, frequency_end_index+1): 
1125                  # normalise amplitudes: divide by number of unique frequencies per unit (partition) 
1126                  amplitudes[unit] += all_amplitudes[frequency]/indices_per_unit 
1127           
1128          return (asarray(amplitudes, float32).reshape(len(amplitudes),1))         
1129           
1130 -    def __call__(self, **params_to_override): 
1131          # override defaults with user defined parameters 
1132          p = ParamOverrides(self, params_to_override) 
1133           
1134          # get the dimensions of the generator sheet. 
1135          self._sheet_dimensions = SheetCoordinateSystem(p.bounds, p.xdensity, p.ydensity).shape 
1136           
1137          # calculate frequency bin divisions. 
1138          self._create_indices(p) 
1139           
1140          # perform a fft to get amplitudes of the composite frequencies. 
1141          return self._get_amplitudes(p) 
1142       
1143       
1144   
1145 -class Spectrogram(PowerSpectrum): 
1146      """ 
1147      Extends PowerSpectrum to provide a temporal buffer, yielding 
1148      a 2D representation of a fixed-width spectrogram. 
1149      """ 
1150      # See comments on PowerSpectrum; this could be instantiated by 
1151      # hand, but in a GUI it would not be useful at present due to the 
1152      # signal parameter. 
1153      __abstract=True 
1154       
1155      seconds_per_timestep=param.Number(default=1.0,doc=""" 
1156          Number of seconds represented by 1 simulation time step.""") 
1157   
1158      sample_window=param.Number(default=1.0,doc=""" 
1159          The length of interval of the signal (in seconds) on which to 
1160          perform the Fourier transform. 
1161   
1162          How much history of the signal to include in the window. 
1163          sample_window > seconds_per_timestep -> window overlap 
1164                                            
1165          The Fourier transform algorithm is most efficient if the 
1166          resulting window_length(sample_window * sample_rate) is a 
1167          power of 2, or can be decomposed into small prime factors; see 
1168          numpy.fft.""") 
1169        
1170 -    def __init__(self, signal, **params): 
1171          # this is used by _create_indices when initialising the spectrogram array. 
1172          self._first_run = True 
1173    
1174          # overload default sample_window and seconds_per_timestep, 
1175          # if a new ones are supplied. 
1176          for parameter,value in params.items(): 
1177              if parameter == "sample_window" or \ 
1178                 parameter == "seconds_per_timestep": 
1179                  setattr(self,parameter,value) 
1180           
1181          if self.sample_window < self.seconds_per_timestep: 
1182              self.warning("sample_window < seconds_per_timestep; some signal will be skipped.") 
1183           
1184          super(Spectrogram, self).__init__(signal, window_increment=self.seconds_per_timestep, 
1185                                             window_length=self.sample_window, **params) 
1186           
1187 -    def _create_indices(self, p): 
1188          super(Spectrogram, self)._create_indices(p) 
1189   
1190          # initalise a new, empty, spectrogram of the size of the generator sheet 
1191          # i.e. allow user to implicitly control spectrogram time history through  
1192          # sheet size. 
1193          if self._first_run: 
1194              self._spectrogram = zeros(self._sheet_dimensions, dtype=float32) 
1195              self._first_run = False 
1196               
1197 -    def __call__(self, **params_to_override): 
1198          # override defaults with user defined parameters. 
1199          p = ParamOverrides(self, params_to_override) 
1200           
1201          # get the dimensions of the generator sheet. 
1202          self._sheet_dimensions = SheetCoordinateSystem(p.bounds, p.xdensity, p.ydensity).shape 
1203           
1204          # calculate frequency bin divisions. 
1205          self._create_indices(p) 
1206                   
1207          # perform a fft to get amplitudes of the composite frequencies. 
1208          amplitudes = self._get_amplitudes(p) 
1209           
1210          # first make sure arrays are of compatible size, then add on  
1211          # latest spectral information to the spectrograms leftmost edge. 
1212          assert shape(amplitudes)[0] == shape(self._spectrogram)[0] 
1213          self._spectrogram = hstack((amplitudes, self._spectrogram)) 
1214                   
1215          # knock off the column on the spectrograms right-most edge, 
1216          # i.e. the oldest spectral information. 
1217          self._spectrogram = self._spectrogram[0:, 0:self._spectrogram.shape[1]-1] 
1218           
1219          # the following print statements are very useful when calibrating sheets, 
1220          # allowing you to calculate how much time history a particular generator 
1221          # sheets x dimension corresponds to. 
1222          #print shape(amplitudes); print shape(self._spectrogram) 
1223           
1224          return self._spectrogram 
1225           
1226                         
1227  __all__ = list(set([k for k,v in locals().items() if isinstance(v,type) and issubclass(v,PatternGenerator)])) 
1228