topo.transferfn.basic

1 """ 2 Simple functions operating on a matrix, potentially modifying it. 3 4 These are useful for neuron output functions, normalization of 5 matrices, etc. 6 7 All of these function objects (callable objects) should work for 8 Numpy array arguments of arbitrary shape. Some may also work for 9 scalars. 10 11 $Id: basic.py 11108 2010-07-05 10:27:54Z ceball $ 12 """ 13 __version__='$Revision: 11108 $' 14 15 import copy 16 17 import param 18 19 import numpy, numpy.random 20 import numpy.oldnumeric as Numeric 21 from numpy import exp,zeros,ones,power 22 23 import topo 24 import topo.base.functionfamily 25 from topo.base.sheet import activity_type 26 from topo.base.arrayutil import clip_lower,clip_upper 27 from topo.base.arrayutil import L2norm, norm 28 from topo.base.functionfamily import TransferFn 29 # Imported here so that all TransferFns will be in the same package 30 from topo.base.functionfamily import IdentityTF 31 32 33 # CEBHACKALERT: these need to respect the mask - which will be passed in. 34 35 # CEBALERT: copy.copy(x) where x is a numpy array should probably be x.copy() 36

37 -class PiecewiseLinear(TransferFn):

38 """ 39 Piecewise-linear TransferFn with lower and upper thresholds. 40 41 Values below the lower_threshold are set to zero, those above 42 the upper threshold are set to 1.0, and those in between are 43 scaled linearly. 44 """ 45 lower_bound = param.Number(default=0.0,softbounds=(0.0,1.0)) 46 upper_bound = param.Number(default=1.0,softbounds=(0.0,1.0)) 47

48 - def __call__(self,x):

49 fact = 1.0/(self.upper_bound-self.lower_bound) 50 x -= self.lower_bound 51 x *= fact 52 x.clip(0.0,1.0,out=x)

53 54 55

56 -class Sigmoid(TransferFn):

57 """ 58 Sigmoidal (logistic) transfer function: 1/(1+exp-(r*x+k)). 59 60 As defined in Jochen Triesch, ICANN 2005, LNCS 3696 pp. 65-70. 61 The parameters control the growth rate (r) and the x position (k) 62 of the exponential. 63 64 This function is a special case of the GeneralizedLogistic 65 function, with parameters r=r, l=0, u=1, m=-k/2r, and b=1. See 66 Richards, F.J. (1959), A flexible growth function for empirical 67 use. J. Experimental Botany 10: 290--300, 1959. 68 http://en.wikipedia.org/wiki/Generalised_logistic_curve 69 """ 70 71 r = param.Number(default=1,doc="Parameter controlling the growth rate") 72 k = param.Number(default=0,doc="Parameter controlling the x-postion") 73

74 - def __call__(self,x):

75 x_orig = copy.copy(x) 76 x *= 0.0 77 x += 1.0 / (1.0 + exp(-(self.r*x_orig+self.k)))

78 79 80

81 -class NakaRushton(TransferFn):

82 #JABALERT: Please write the equation into words in the docstring, as in Sigmoid. 83 """ 84 Naka-Rushton curve. 85 86 From Naka, K. and Rushton, W. (1996), S-potentials from luminosity 87 units in the retina of fish (Cyprinidae). J. Physiology 185:587-599. 88 89 The Naka-Rushton curve has been shown to be a good approximation 90 of constrast gain control in cortical neurons. The input of the 91 curve is usually contrast, but under the assumption that the 92 firing rate of a model neuron is directly proportional to the 93 contrast, it can be used as a TransferFn for a Sheet. 94 95 The parameter c50 corresponds to the contrast at which the half of 96 the maximal output is reached. For a Sheet TransferFn this translates 97 to the input for which a neuron will respond with activity 0.5. 98 """ 99 100 c50 = param.Number(default=0.1, doc=""" 101 The input of the neuron at which it responds at half of its maximal firing rate (1.0).""") 102 103 e = param.Number(default=1.0,doc="""The exponent of the input x.""") 104 105 #JABALERT: (pow(x_orig,self.e) should presumably be done only once, using a temporary

106 - def __call__(self,x):

107 #print 'A:', x 108 #print 'B:', pow(x,self.e) / (pow(x,self.e) + pow(self.c50,self.e)) 109 x_orig = copy.copy(x) 110 x *= 0 111 x += pow(x_orig,self.e) / (pow(x_orig,self.e) + pow(self.c50,self.e))

112 113 114

115 -class GeneralizedLogistic(TransferFn):

116 """ 117 The generalized logistic curve (Richards' curve): y = l + (u /(1 + b * exp(-r*(x-2*m))^(1/b))). 118 119 The logistic curve is a flexible function for specifying a 120 nonlinear growth curve using five parameters: 121 122 * l: the lower asymptote 123 * u: the upper asymptote minus l 124 * m: the time of maximum growth 125 * r: the growth rate 126 * b: affects near which asymptote maximum growth occurs 127 128 From Richards, F.J. (1959), A flexible growth function for empirical 129 use. J. Experimental Botany 10: 290--300. 130 http://en.wikipedia.org/wiki/Generalised_logistic_curve 131 """ 132 133 # JABALERT: Reword these to say what they are, not what they 134 # control, if they are anything that can be expressed naturally. 135 # E.g. is l a parameter controlling the lower asymptote, or is it 136 # simply the lower asymptote? If it's the lower asymptote, say 137 # doc="Lower asymptote.". Only if the parameter's relationship to 138 # what it controls is very indirect should it be worded as below. 139 l = param.Number(default=1,doc="Parameter controlling the lower asymptote.") 140 u = param.Number(default=1,doc="Parameter controlling the upper asymptote (upper asymptote minus lower asymptote.") 141 m = param.Number(default=1,doc="Parameter controlling the time of maximum growth.") 142 r = param.Number(default=1,doc="Parameter controlling the growth rate.") 143 b = param.Number(default=1,doc="Parameter which affects near which asymptote maximum growth occurs.") 144

145 - def __call__(self,x):

146 x_orig = copy.copy(x) 147 x *= 0.0 148 x += self.l + ( self.u /(1 + self.b*exp(-self.r *(x_orig - 2*self.m))**(1 / self.b)) )

149 150 151

152 -class DivisiveNormalizeL1(TransferFn):

153 """ 154 TransferFn that divides an array by its L1 norm. 155 156 This operation ensures that the sum of the absolute values of the 157 array is equal to the specified norm_value, rescaling each value 158 to make this true. The array is unchanged if the sum of absolute 159 values is zero. For arrays of non-negative values where at least 160 one is non-zero, this operation is equivalent to a divisive sum 161 normalization. 162 """ 163 norm_value = param.Number(default=1.0) 164

165 - def __call__(self,x):

166 """L1-normalize the input array, if it has a nonzero sum.""" 167 current_sum = 1.0*Numeric.sum(abs(x.ravel())) 168 if current_sum != 0: 169 factor = (self.norm_value/current_sum) 170 x *= factor

171 172 173

174 -class DivisiveNormalizeL2(TransferFn):

175 """ 176 TransferFn to divide an array by its Euclidean length (aka its L2 norm). 177 178 For a given array interpreted as a flattened vector, keeps the 179 Euclidean length of the vector at a specified norm_value. 180 """ 181 norm_value = param.Number(default=1.0) 182

183 - def __call__(self,x):

184 tot = 1.0*L2norm(x.ravel()) 185 if tot != 0: 186 factor = (self.norm_value/tot) 187 x *= factor

188 189 190

191 -class DivisiveNormalizeLinf(TransferFn):

192 """ 193 TransferFn to divide an array by its L-infinity norm 194 (i.e. the maximum absolute value of its elements). 195 196 For a given array interpreted as a flattened vector, scales the 197 elements divisively so that the maximum absolute value is the 198 specified norm_value. 199 200 The L-infinity norm is also known as the divisive infinity norm 201 and Chebyshev norm. 202 """ 203 norm_value = param.Number(default=1.0) 204

205 - def __call__(self,x):

206 tot = 1.0*(numpy.abs(x)).max() 207 if tot != 0: 208 factor = (self.norm_value/tot) 209 x *= factor

210 211 212

213 -class DivisiveNormalizeLp(TransferFn):

214 """ 215 TransferFn to divide an array by its Lp-Norm, where p is specified. 216 217 For a parameter p and a given array interpreted as a flattened 218 vector, keeps the Lp-norm of the vector at a specified norm_value. 219 Faster versions are provided separately for the typical L1-norm 220 and L2-norm cases. Defaults to be the same as an L2-norm, i.e., 221 DivisiveNormalizeL2. 222 """ 223 p = param.Number(default=2) 224 norm_value = param.Number(default=1.0) 225

226 - def __call__(self,x):

227 tot = 1.0*norm(x.ravel(),self.p) 228 if tot != 0: 229 factor = (self.norm_value/tot) 230 x *=factor

231 232 233

234 -class HalfRectifyAndSquare(TransferFn):

235 """ 236 Transfer function that applies a half-wave rectification (clips at zero) 237 and then squares the values. 238 """ 239 t = param.Number(default=0.0,doc=""" 240 The threshold at which output becomes non-zero.""") 241

242 - def __call__(self,x):

243 x -= self.t 244 clip_lower(x,0) 245 x *= x

246 247

248 -class HalfRectifyAndPower(TransferFn):

249 """ 250 Transfer function that applies a half-wave rectification (i.e., 251 clips at zero), and then raises the result to the e-th power 252 (where the exponent e can be selected arbitrarily). 253 """ 254 e = param.Number(default=2.0,doc=""" 255 The exponent to which the thresholded value is raised.""") 256 257 t = param.Number(default=0.0,doc=""" 258 The threshold level subtracted from x.""") 259

260 - def __call__(self,x):

261 x -= self.t 262 clip_lower(x,0) 263 a = power(x,self.e) 264 x*=0 265 x+=a

266

267 -class ExpLinear(TransferFn):

268 """ 269 Transfer function that is exponential until t from which point it is linear. 270 """ 271 e = param.Number(default=1.0,doc=""" 272 The exponent of the exponetial part of the curve""") 273 t1 = param.Number(default=0.5,doc=""" 274 The threshold level where function becomes non-zero""") 275 t2 = param.Number(default=1.0,doc=""" 276 The threshold level at which curve becomes linear""") 277 278 a = param.Number(default=1.0,doc=""" 279 The overall scaling of the function""") 280

281 - def __call__(self,x):

282 x-=self.t1 283 clip_lower(x,0) 284 z = (x>=self.t2)*(1.0*exp(self.e*(-self.t2))+(x-self.t2)) + (x<self.t2)*(exp(self.e*(x-self.t2))-exp(self.e*(-self.t2))) 285 x*=0 286 x+=self.a*z

287 288 289

290 -class Square(TransferFn):

291 """Transfer function that applies a squaring nonlinearity.""" 292

293 - def __call__(self,x):

294 x *= x

295 296 297

298 -class BinaryThreshold(TransferFn):

299 """ 300 Forces all values below a threshold to zero, and above it to 1.0. 301 """ 302 threshold = param.Number(default=0.25, doc="Decision point for determining binary value.") 303

304 - def __call__(self,x):

305 above_threshold = x>=self.threshold 306 x *= 0.0 307 x += above_threshold

308 309

310 -class Threshold(TransferFn):

311 """ 312 Forces all values below a threshold to zero, and leaves others unchanged. 313 """ 314 threshold = param.Number(default=0.25, doc="Decision point for determining values to clip.") 315

316 - def __call__(self,x):

317 clip_upper(x,self.threshold)

318 319 320 # JAALERT: rename to something like PlasticTransferFn

321 -class TransferFnWithState(TransferFn):

322 """ 323 Abstract base class for TransferFns that need to maintain a self.plastic parameter. 324 325 These TransferFns typically maintain some form of internal history 326 or other state from previous calls, which can be disabled by 327 override_plasticity_state(). 328 """ 329 330 plastic = param.Boolean(default=True, doc=""" 331 Whether or not to update the internal state on each call. 332 Allows plasticity to be turned off during analysis, and then re-enabled.""") 333 334 __abstract = True 335

336 - def __init__(self,**params):

337 super(TransferFnWithState,self).__init__(**params) 338 self._plasticity_setting_stack = []

339 340

341 - def override_plasticity_state(self, new_plasticity_state):

342 """ 343 Temporarily disable plasticity of internal state. 344 345 This function should be implemented by all subclasses so that 346 after a call, the output should always be the same for any 347 given input pattern (apart from true randomness or other 348 differences that do not depend on an internal state), and no 349 call should have any effect that persists after a subsequent 350 restore_plasticity_state() call. 351 352 By default, simply saves a copy of the 'plastic' parameter to 353 an internal stack (so that it can be restored by 354 restore_plasticity_state()), and then sets the plastic 355 parameter to the given value (True or False). 356 """ 357 self._plasticity_setting_stack.append(self.plastic) 358 self.plastic=new_plasticity_state

359 360

361 - def restore_plasticity_state(self):

362 """ 363 Re-enable plasticity of internal state after an override_plasticity_state call. 364 365 This function should be implemented by all subclasses to 366 remove the effect of the most recent override_plasticity_state call, 367 i.e. to reenable changes to the internal state, without any 368 lasting effect from the time during which plasticity was disabled. 369 370 By default, simply restores the last saved value of the 371 'plastic' parameter. 372 """ 373 self.plastic = self._plasticity_setting_stack.pop()

374

375 - def state_push(self):

376 """ 377 Save the current state onto a stack, to be restored using state_pop. 378 379 Subclasses must implement state_push and state_pop if they 380 store any lasting state across invocations, so that the result 381 of state_pop will be the state that was present at the 382 previous state_push. 383 """ 384 pass

385

386 - def state_pop(self):

387 """ 388 Restore the state saved by the most recent state_push call. 389 """ 390 pass

391 392 393 # CB: it's not ideal that all TransferFnWithRandomState fns have 394 # the plastic stuff (from TransferFnWithState).

395 -class TransferFnWithRandomState(TransferFnWithState):

396 """ 397 Abstract base class for TransferFns that use a random number generator. 398 """ 399 400 random_generator = param.Parameter( 401 default=numpy.random.RandomState(seed=(10,10)),doc= 402 """ 403 numpy's RandomState provides methods for generating random 404 numbers (see RandomState's help for more information). 405 406 Note that all instances of subclasses of 407 TransferFnWithRandomState will share this RandomState object, 408 and hence its state. To create an instance of an 409 TransferFnWithRandomState subclass that has its own state, set 410 this parameter on the instance to a new RandomState instance. 411 """) 412 413 __abstract = True 414

415 - def __init__(self,**params):

416 super(TransferFnWithRandomState,self).__init__(**params) 417 self.__random_generators_stack = []

418

419 - def state_push(self):

420 """ 421 Save the current random number generator (onto the stack), 422 replacing it with a copy. 423 """ 424 self.__random_generators_stack.append(self.random_generator) 425 self.random_generator=copy.copy(self.random_generator) 426 super(TransferFnWithRandomState,self).state_push()

427

428 - def state_pop(self):

429 """ 430 Retrieve the previous random number generator from the stack. 431 """ 432 self.random_generator = self.__random_generators_stack.pop() 433 super(TransferFnWithRandomState,self).state_push()

434 435 436

437 -class PoissonSample(TransferFnWithRandomState):

438 """ 439 Simulate Poisson-distributed activity with specified mean values. 440 441 This transfer function interprets each matrix value as the 442 (potentially scaled) rate of a Poisson process and replaces it 443 with a sample from the appropriate Poisson distribution. 444 445 To allow the matrix to contain values in a suitable range (such as 446 [0.0,1.0]), the input matrix is scaled by the parameter in_scale, 447 and the baseline_rate is added before sampling. After sampling, 448 the output value is then scaled by out_scale. The function thus 449 performs this transformation:: 450 451 x <- P(in_scale * x + baseline_rate) * out_scale 452 453 where x is a matrix value and P(r) samples from a Poisson 454 distribution with rate r. 455 """ 456 457 in_scale = param.Number(default=1.0,doc=""" 458 Amount by which to scale the input.""") 459 460 baseline_rate = param.Number(default=0.0,doc=""" 461 Constant to add to the input after scaling, resulting in a baseline 462 Poisson process rate.""") 463 464 out_scale = param.Number(default=1.0,doc=""" 465 Amount by which to scale the output (e.g. 1.0/in_scale).""") 466

467 - def __call__(self,x):

468 469 x *= self.in_scale 470 x += self.baseline_rate 471 sample = self.random_generator.poisson(x,x.shape) 472 x *= 0.0 473 x += sample 474 x *= self.out_scale

475 476 477 478 # CEBALERT: I think this TF is a bit misleading because it does not 479 # alter the input array. Is a TransferFn the right thing to use to 480 # record an average activity? Could consider replacing this with my 481 # "attribute tracking" code (some of which is not yet in SVN).

482 -class ActivityAveragingTF(TransferFnWithState):

483 """ 484 Calculates the average of the input activity. 485 486 The average is calculated as an exponential moving average, where 487 the weighting for each older data point decreases exponentially. 488 The degree of weighing for the previous values is expressed as a 489 constant smoothing factor. 490 491 The plastic parameter allows the updating of the average values 492 to be disabled temporarily, e.g. while presenting test patterns. 493 """ 494 495 step = param.Number(default=1, doc=""" 496 How often to update the average. 497 498 For instance, step=1 means to update it every time this OF is 499 called; step=2 means to update it every other time.""") 500 501 smoothing = param.Number(default=0.9997, doc=""" 502 The degree of weighting for the previous average, when calculating the new average.""") 503 504 initial_average=param.Number(default=0, doc="Starting value for the average activity.") 505 506

507 - def __init__(self,**params):

508 super(ActivityAveragingTF,self).__init__(**params) 509 self.n_step = 0 510 self.x_avg=None

511 512

513 - def __call__(self,x):

514 if self.x_avg is None: 515 self.x_avg=self.initial_average*ones(x.shape, activity_type) 516 517 # Collect values on each appropriate step 518 519 self.n_step += 1 520 if self.n_step == self.step: 521 self.n_step = 0 522 if self.plastic: 523 self.x_avg = (1.0-self.smoothing)*x + self.smoothing*self.x_avg

524 525

526 -class HomeostaticMaxEnt(TransferFnWithRandomState):

527 """ 528 Implementation of homeostatic intrinsic plasticity from Jochen Triesch, 529 ICANN 2005, LNCS 3696 pp.65-70. 530 531 A sigmoid activation function is adapted automatically to achieve 532 desired average firing rate and approximately exponential 533 distribution of firing rates (for the maximum possible entropy). 534 535 Note that this TransferFn has state, so the history of calls to it 536 will affect future behavior. The plastic parameter can be used 537 to disable changes to the state. 538 539 Also calculates average activity as useful debugging information, 540 for use with ValueTrackingOutoutFn Average activity is calculated as 541 an exponential moving average with a smoothing factor (smoothing). 542 For more information see: 543 NIST/SEMATECH e-Handbook of Statistical Methods, Single Exponential Smoothing 544 http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc431.htm 545 """ 546 547 eta = param.Number(default=0.0002,doc="Learning rate for homeostatic plasticity.") 548 549 mu = param.Number(default=0.01,doc="Target average firing rate.") 550 551 smoothing = param.Number(default=0.9997, doc=""" 552 Weighting of previous activity vs. current activity when calculating the average.""") 553 554 a_init = param.Parameter(default=None,doc="Multiplicative parameter controlling the exponential.") 555 556 b_init = param.Parameter(default=None,doc="Additive parameter controlling the exponential.") 557 558 step = param.Number(default=1, doc=""" How often to update the a and b parameters. For instance, step=1 means to update it every time this OF is 559 called; step=2 means to update it every other time.""") 560

561 - def __init__(self,**params):

562 super(HomeostaticMaxEnt,self).__init__(**params) 563 self.first_call = True 564 self.n_step=0 565 self.__current_state_stack=[] 566 self.a=None 567 self.b=None 568 self.y_avg=None

569

570 - def __call__(self,x):

571 if self.first_call: 572 self.first_call = False 573 if self.a_init==None: 574 self.a = self.random_generator.uniform(low=10, high=20,size=x.shape) 575 else: 576 self.a = ones(x.shape, x.dtype.char) * self.a_init 577 if self.b_init==None: 578 self.b = self.random_generator.uniform(low=-8.0, high=-4.0,size=x.shape) 579 else: 580 self.b = ones(x.shape, x.dtype.char) * self.b_init 581 self.y_avg = zeros(x.shape, x.dtype.char) 582 583 # Apply sigmoid function to x, resulting in what Triesch calls y 584 x_orig = copy.copy(x) 585 586 x *= 0.0 587 x += 1.0 / (1.0 + exp(-(self.a*x_orig + self.b))) 588 589 590 self.n_step += 1 591 if self.n_step == self.step: 592 self.n_step = 0 593 if self.plastic: 594 self.y_avg = (1.0-self.smoothing)*x + self.smoothing*self.y_avg #Calculate average for use in debugging only 595 596 # Update a and b 597 self.a += self.eta * (1.0/self.a + x_orig - (2.0 + 1.0/self.mu)*x_orig*x + x_orig*x*x/self.mu) 598 self.b += self.eta * (1.0 - (2.0 + 1.0/self.mu)*x + x*x/self.mu)

599 600

601 - def state_push(self):

602 self.__current_state_stack.append((copy.copy(self.a), copy.copy(self.b), copy.copy(self.y_avg), copy.copy(self.first_call))) 603 super(HomeostaticMaxEnt,self).state_push()

604 605

606 - def state_pop(self):

607 self.a, self.b, self.y_avg, self.first_call = self.__current_state_stack.pop() 608 super(HomeostaticMaxEnt,self).state_pop()

609 610

611 -class ScalingTF(TransferFnWithState):

612 """ 613 Scales input activity based on the current average activity (x_avg). 614 615 The scaling is calculated to bring x_avg for each unit closer to a 616 specified target average. Calculates a scaling factor that is 617 greater than 1 if x_avg is less than the target and less than 1 if 618 x_avg is greater than the target, and multiplies the input 619 activity by this scaling factor. 620 621 The plastic parameter allows the updating of the average values 622 to be disabled temporarily, e.g. while presenting test patterns. 623 """ 624 625 target = param.Number(default=0.01, doc=""" 626 Target average activity for each unit.""") 627 628 step=param.Number(default=1, doc=""" 629 How often to calculate the average activity and scaling factor.""") 630 631 smoothing = param.Number(default=0.9997, doc=""" 632 Determines the degree of weighting of previous activity vs. 633 current activity when calculating the average.""") 634 635

636 - def __init__(self,**params):

637 super(ScalingTF,self).__init__(**params) 638 self.n_step = 0 639 self.x_avg=None 640 self.sf=None

641

642 - def __call__(self,x):

643 644 if self.x_avg is None: 645 self.x_avg=self.target*ones(x.shape, activity_type) 646 if self.sf is None: 647 self.sf=ones(x.shape, activity_type) 648 649 # Collect values on each appropriate step 650 651 self.n_step += 1 652 if self.n_step == self.step: 653 self.n_step = 0 654 if self.plastic: 655 self.sf *= 0.0 656 self.sf += self.target/self.x_avg 657 self.x_avg = (1.0-self.smoothing)*x + self.smoothing*self.x_avg 658 659 x *= self.sf

660 661 662 663 664

665 -class Hysteresis(TransferFnWithState):

666 """ 667 Smoothly interpolates a matrix between simulation time steps, with 668 exponential falloff. 669 """ 670 671 time_constant = param.Number(default=0.3,doc=""" 672 Controls the time scale of the interpolation.""") 673

674 - def __init__(self,**params):

675 super(Hysteresis,self).__init__(**params) 676 self.first_call = True 677 self.__current_state_stack=[] 678 self.old_a = 0 679 topo.base.functionfamily.PatternDrivenAnalysis.pre_presentation_hooks.append(self.reset)

680

681 - def __call__(self,x):

682 if self.first_call is True: 683 self.old_a = x.copy() * 0.0 684 self.first_call = False 685 686 #if (float(topo.sim.time()) %1.0) <= 0.15: self.old_a = self.old_a* 0 687 new_a = x.copy() 688 self.old_a = self.old_a + (new_a - self.old_a)*self.time_constant 689 x*=0 690 x += self.old_a

691

692 - def reset(self):

693 self.old_a *= 0

694

695 - def state_push(self):

696 self.__current_state_stack.append((copy.copy(self.old_a), copy.copy(self.first_call))) 697 super(Hysteresis,self).state_push()

698

699 - def state_pop(self):

700 self.old_a,self.first_call = self.__current_state_stack.pop() 701 super(Hysteresis,self).state_pop()

702 703 704 __all__ = list(set([k for k,v in locals().items() if isinstance(v,type) and issubclass(v,TransferFn)])) 705

Source Code for Module topo.transferfn.basic