| 316 | | OU_generator = _OU_generator_python |
|---|
| | 325 | |
|---|
| | 326 | def _OU_generator_python2(self, dt, tau, sigma, y0, t_start=0.0, t_stop=1000.0, array=False,time_it=False): |
|---|
| | 327 | """ |
|---|
| | 328 | Generates an Orstein Ulbeck process using the forward euler method. The function returns |
|---|
| | 329 | an AnalogSignal object |
|---|
| | 330 | |
|---|
| | 331 | Inputs: |
|---|
| | 332 | dt - the time resolution in milliseconds of th signal |
|---|
| | 333 | tau - the correlation time in milliseconds |
|---|
| | 334 | sigma - std dev of the process |
|---|
| | 335 | y0 - initial value of the process, at t_start |
|---|
| | 336 | t_start - start time in milliseconds |
|---|
| | 337 | t_stop - end time in milliseconds |
|---|
| | 338 | array - if True, a numpy array of the OU signal is returned, |
|---|
| | 339 | rather than an AnalogSignal object. |
|---|
| | 340 | |
|---|
| | 341 | Examples: |
|---|
| | 342 | >> stgen.OU_generator(0.1, 2, 3, 0, 0, 10000) |
|---|
| | 343 | """ |
|---|
| | 344 | |
|---|
| | 345 | import time |
|---|
| | 346 | |
|---|
| | 347 | if time_it: |
|---|
| | 348 | t1 = time.time() |
|---|
| | 349 | |
|---|
| | 350 | t = numpy.arange(t_start,t_stop,dt) |
|---|
| | 351 | N = len(t) |
|---|
| | 352 | y = numpy.zeros(N,float) |
|---|
| | 353 | y[0] = y0 |
|---|
| | 354 | fac = dt/tau |
|---|
| | 355 | gauss = fac*y0+numpy.sqrt(2*fac)*sigma*self.rng.standard_normal(N-1) |
|---|
| | 356 | mfac = 1-fac |
|---|
| | 357 | |
|---|
| | 358 | # python loop... bad+slow! |
|---|
| | 359 | for i in xrange(1,N): |
|---|
| | 360 | idx = i-1 |
|---|
| | 361 | y[i] = y[idx]*mfac+gauss[idx] |
|---|
| | 362 | |
|---|
| | 363 | if time_it: |
|---|
| | 364 | print time.time()-t1 |
|---|
| | 365 | |
|---|
| | 366 | if array: |
|---|
| | 367 | return (y,t) |
|---|
| | 368 | |
|---|
| | 369 | result = AnalogSignal(y,dt,t_start=0,t_stop=t_stop-t_start) |
|---|
| | 370 | result.time_offset(t_start) |
|---|
| | 371 | return result |
|---|
| | 372 | |
|---|
| | 373 | # use slow python implementation for the time being |
|---|
| | 374 | # TODO: provide optimized C/weave implementation if possible |
|---|
| | 375 | |
|---|
| | 376 | |
|---|
| | 377 | def OU_generator_weave1(self, dt,tau,sigma,y0,t_start=0.0,t_stop=1000.0,time_it=False): |
|---|
| | 378 | """ generates an OU process using forward euler |
|---|
| | 379 | dt = resolution |
|---|
| | 380 | tau = correlation time |
|---|
| | 381 | sigma = std dev of process |
|---|
| | 382 | y0 = mean/initial value |
|---|
| | 383 | tsim = how long to simulate |
|---|
| | 384 | """ |
|---|
| | 385 | import scipy.weave |
|---|
| | 386 | |
|---|
| | 387 | import time |
|---|
| | 388 | |
|---|
| | 389 | if time_it: |
|---|
| | 390 | t1 = time.time() |
|---|
| | 391 | |
|---|
| | 392 | |
|---|
| | 393 | t = numpy.arange(t_start,t_stop,dt) |
|---|
| | 394 | N = len(t) |
|---|
| | 395 | y = numpy.zeros(N,float) |
|---|
| | 396 | y[0] = y0 |
|---|
| | 397 | fac = dt/tau |
|---|
| | 398 | gauss = fac*y0+numpy.sqrt(2*fac)*sigma*self.rng.standard_normal(N-1) |
|---|
| | 399 | |
|---|
| | 400 | |
|---|
| | 401 | # python loop... bad+slow! |
|---|
| | 402 | #for i in xrange(1,len(t)): |
|---|
| | 403 | # y[i] = y[i-1]+dt/tau*(y0-y[i-1])+numpy.sqrt(2*dt/tau)*sigma*numpy.random.normal() |
|---|
| | 404 | |
|---|
| | 405 | # use weave instead |
|---|
| | 406 | |
|---|
| | 407 | code = """ |
|---|
| | 408 | |
|---|
| | 409 | double f = 1.0-fac; |
|---|
| | 410 | |
|---|
| | 411 | for(int i=1;i<Ny[0];i++) { |
|---|
| | 412 | y(i) = y(i-1)*f + gauss(i-1); |
|---|
| | 413 | } |
|---|
| | 414 | """ |
|---|
| | 415 | |
|---|
| | 416 | scipy.weave.inline(code,['y', 'gauss', 'fac'], |
|---|
| | 417 | type_converters=scipy.weave.converters.blitz) |
|---|
| | 418 | |
|---|
| | 419 | |
|---|
| | 420 | if time_it: |
|---|
| | 421 | print 'Elapsed ',time.time()-t1,' seconds.' |
|---|
| | 422 | |
|---|
| | 423 | if array: |
|---|
| | 424 | return (y,t) |
|---|
| | 425 | |
|---|
| | 426 | result = AnalogSignal(y,dt,t_start=0,t_stop=t_stop-t_start) |
|---|
| | 427 | result.time_offset(t_start) |
|---|
| | 428 | return result |
|---|
| | 429 | |
|---|
| | 430 | |
|---|
| | 431 | |
|---|
| | 432 | OU_generator = _OU_generator_python2 |
|---|
