1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6'''This module provides convenience objects to handle frequency responses.''' 

7 

8from __future__ import print_function, division, absolute_import 

9 

10import math 

11import logging 

12 

13import numpy as num 

14from scipy import signal 

15 

16from pyrocko import util, evalresp 

17from pyrocko.guts import Object, Float, Int, String, Complex, Tuple, List, \ 

18 StringChoice, Bool 

19from pyrocko.guts_array import Array 

20 

21try: 

22 newstr = unicode 

23except NameError: 

24 newstr = str 

25 

26 

27guts_prefix = 'pf' 

28 

29logger = logging.getLogger('pyrocko.response') 

30 

31 

32def asarray_1d(x, dtype): 

33 if isinstance(x, (list, tuple)) and x and isinstance(x[0], (str, newstr)): 

34 return num.asarray(list(map(dtype, x)), dtype=dtype) 

35 else: 

36 a = num.asarray(x, dtype=dtype) 

37 if not a.ndim == 1: 

38 raise ValueError('could not convert to 1D array') 

39 return a 

40 

41 

42def finalize_construction(breakpoints): 

43 breakpoints.sort() 

44 breakpoints_out = [] 

45 f_last = None 

46 for f, c in breakpoints: 

47 if f_last is not None and f == f_last: 

48 breakpoints_out[-1][1] += c 

49 else: 

50 breakpoints_out.append([f, c]) 

51 

52 f_last = f 

53 

54 breakpoints_out = [(f, c) for (f, c) in breakpoints_out if c != 0] 

55 return breakpoints_out 

56 

57 

58class FrequencyResponseCheckpoint(Object): 

59 frequency = Float.T() 

60 value = Float.T() 

61 

62 

63class IsNotScalar(Exception): 

64 pass 

65 

66 

67class FrequencyResponse(Object): 

68 ''' 

69 Evaluates frequency response at given frequencies. 

70 ''' 

71 

72 checkpoints = List.T(FrequencyResponseCheckpoint.T()) 

73 

74 def evaluate(self, freqs): 

75 return num.ones(freqs.size, dtype=complex) 

76 

77 def evaluate1(self, freq): 

78 return self.evaluate(num.atleast_1d(freq))[0] 

79 

80 def is_scalar(self): 

81 ''' 

82 Check if this is a flat response. 

83 ''' 

84 

85 if type(self) is FrequencyResponse: 

86 return True 

87 else: 

88 return False # default for derived classes 

89 

90 def get_scalar(self): 

91 ''' 

92 Get factor if this is a flat response. 

93 ''' 

94 if type(self) is FrequencyResponse: 

95 return 1.0 

96 else: 

97 raise IsNotScalar() # default for derived classes 

98 

99 def get_fmax(self): 

100 return None 

101 

102 def construction(self): 

103 return [] 

104 

105 

106class Gain(FrequencyResponse): 

107 ''' 

108 A flat frequency response. 

109 ''' 

110 

111 constant = Complex.T(default=1.0+0j) 

112 

113 def evaluate(self, freqs): 

114 return util.num_full_like(freqs, self.constant, dtype=complex) 

115 

116 def is_scalar(self): 

117 return True 

118 

119 def get_scalar(self): 

120 return self.constant 

121 

122 

123class Evalresp(FrequencyResponse): 

124 ''' 

125 Calls evalresp and generates values of the instrument response transfer 

126 function. 

127 

128 :param respfile: response file in evalresp format 

129 :param trace: trace for which the response is to be extracted from the file 

130 :param target: ``'dis'`` for displacement or ``'vel'`` for velocity 

131 ''' 

132 

133 respfile = String.T() 

134 nslc_id = Tuple.T(4, String.T()) 

135 target = String.T(default='dis') 

136 instant = Float.T() 

137 stages = Tuple.T(2, Int.T(), optional=True) 

138 

139 def __init__( 

140 self, 

141 respfile, 

142 trace=None, 

143 target='dis', 

144 nslc_id=None, 

145 time=None, 

146 stages=None, 

147 **kwargs): 

148 

149 if trace is not None: 

150 nslc_id = trace.nslc_id 

151 time = (trace.tmin + trace.tmax) / 2. 

152 

153 FrequencyResponse.__init__( 

154 self, 

155 respfile=respfile, 

156 nslc_id=nslc_id, 

157 instant=time, 

158 target=target, 

159 stages=stages, 

160 **kwargs) 

161 

162 def evaluate(self, freqs): 

163 network, station, location, channel = self.nslc_id 

164 if self.stages is None: 

165 stages = (-1, 0) 

166 else: 

167 stages = self.stages[0]+1, self.stages[1] 

168 

169 x = evalresp.evalresp( 

170 sta_list=station, 

171 cha_list=channel, 

172 net_code=network, 

173 locid=location, 

174 instant=self.instant, 

175 freqs=freqs, 

176 units=self.target.upper(), 

177 file=self.respfile, 

178 start_stage=stages[0], 

179 stop_stage=stages[1], 

180 rtype='CS') 

181 

182 transfer = x[0][4] 

183 return transfer 

184 

185 

186class InverseEvalresp(FrequencyResponse): 

187 ''' 

188 Calls evalresp and generates values of the inverse instrument response for 

189 deconvolution of instrument response. 

190 

191 :param respfile: response file in evalresp format 

192 :param trace: trace for which the response is to be extracted from the file 

193 :param target: ``'dis'`` for displacement or ``'vel'`` for velocity 

194 ''' 

195 

196 respfile = String.T() 

197 nslc_id = Tuple.T(4, String.T()) 

198 target = String.T(default='dis') 

199 instant = Float.T() 

200 

201 def __init__(self, respfile, trace, target='dis', **kwargs): 

202 FrequencyResponse.__init__( 

203 self, 

204 respfile=respfile, 

205 nslc_id=trace.nslc_id, 

206 instant=(trace.tmin + trace.tmax)/2., 

207 target=target, 

208 **kwargs) 

209 

210 def evaluate(self, freqs): 

211 network, station, location, channel = self.nslc_id 

212 x = evalresp.evalresp(sta_list=station, 

213 cha_list=channel, 

214 net_code=network, 

215 locid=location, 

216 instant=self.instant, 

217 freqs=freqs, 

218 units=self.target.upper(), 

219 file=self.respfile, 

220 rtype='CS') 

221 

222 transfer = x[0][4] 

223 return 1./transfer 

224 

225 

226def aslist(x): 

227 if x is None: 

228 return [] 

229 

230 try: 

231 return list(x) 

232 except TypeError: 

233 return [x] 

234 

235 

236class PoleZeroResponse(FrequencyResponse): 

237 ''' 

238 Evaluates frequency response from pole-zero representation. 

239 

240 :param zeros: positions of zeros 

241 :type zeros: list of complex 

242 :param poles: positions of poles 

243 :type poles: list of complex 

244 :param constant: gain factor 

245 :type constant: complex 

246 

247 :: 

248 

249 (j*2*pi*f - zeros[0]) * (j*2*pi*f - zeros[1]) * ... 

250 T(f) = constant * ---------------------------------------------------- 

251 (j*2*pi*f - poles[0]) * (j*2*pi*f - poles[1]) * ... 

252 

253 

254 The poles and zeros should be given as angular frequencies, not in Hz. 

255 ''' 

256 

257 zeros = List.T(Complex.T()) 

258 poles = List.T(Complex.T()) 

259 constant = Complex.T(default=1.0+0j) 

260 

261 def __init__( 

262 self, 

263 zeros=None, 

264 poles=None, 

265 constant=1.0+0j, 

266 **kwargs): 

267 

268 if zeros is None: 

269 zeros = [] 

270 if poles is None: 

271 poles = [] 

272 

273 FrequencyResponse.__init__( 

274 self, 

275 zeros=aslist(zeros), 

276 poles=aslist(poles), 

277 constant=constant, 

278 **kwargs) 

279 

280 def evaluate(self, freqs): 

281 if hasattr(signal, 'freqs_zpk'): # added in scipy 0.19.0 

282 return signal.freqs_zpk( 

283 self.zeros, self.poles, self.constant, freqs*2.*num.pi)[1] 

284 else: 

285 jomeg = 1.0j * 2.*num.pi*freqs 

286 

287 a = num.ones(freqs.size, dtype=complex)*self.constant 

288 for z in self.zeros: 

289 a *= jomeg-z 

290 for p in self.poles: 

291 a /= jomeg-p 

292 

293 return a 

294 

295 def is_scalar(self): 

296 return len(self.zeros) == 0 and len(self.poles) == 0 

297 

298 def get_scalar(self): 

299 ''' 

300 Get factor if this is a flat response. 

301 ''' 

302 if self.is_scalar(): 

303 return self.constant 

304 else: 

305 raise IsNotScalar() 

306 

307 def inverse(self): 

308 return PoleZeroResponse( 

309 poles=list(self.zeros), 

310 zeros=list(self.poles), 

311 constant=1.0/self.constant) 

312 

313 def to_analog(self): 

314 b, a = signal.zpk2tf(self.zeros, self.poles, self.constant) 

315 return AnalogFilterResponse(aslist(b), aslist(a)) 

316 

317 def to_digital(self, deltat, method='bilinear'): 

318 from scipy.signal import cont2discrete, zpk2tf 

319 

320 z, p, k, _ = cont2discrete( 

321 (self.zeros, self.poles, self.constant), 

322 deltat, method=method) 

323 

324 b, a = zpk2tf(z, p, k) 

325 

326 return DigitalFilterResponse(b, a, deltat) 

327 

328 def to_digital_polezero(self, deltat, method='bilinear'): 

329 from scipy.signal import cont2discrete 

330 

331 z, p, k, _ = cont2discrete( 

332 (self.zeros, self.poles, self.constant), 

333 deltat, method=method) 

334 

335 return DigitalPoleZeroResponse(z, p, k, deltat) 

336 

337 def construction(self): 

338 breakpoints = [] 

339 for zero in self.zeros: 

340 f = abs(zero) / (2.*math.pi) 

341 breakpoints.append((f, 1)) 

342 

343 for pole in self.poles: 

344 f = abs(pole) / (2.*math.pi) 

345 breakpoints.append((f, -1)) 

346 

347 return finalize_construction(breakpoints) 

348 

349 

350class DigitalPoleZeroResponse(FrequencyResponse): 

351 ''' 

352 Evaluates frequency response from digital filter pole-zero representation. 

353 

354 :param zeros: positions of zeros 

355 :type zeros: list of complex 

356 :param poles: positions of poles 

357 :type poles: list of complex 

358 :param constant: gain factor 

359 :type constant: complex 

360 :param deltat: sampling interval 

361 :type deltat: float 

362 

363 The poles and zeros should be given as angular frequencies, not in Hz. 

364 ''' 

365 

366 zeros = List.T(Complex.T()) 

367 poles = List.T(Complex.T()) 

368 constant = Complex.T(default=1.0+0j) 

369 deltat = Float.T() 

370 

371 def __init__( 

372 self, 

373 zeros=None, 

374 poles=None, 

375 constant=1.0+0j, 

376 deltat=None, 

377 **kwargs): 

378 

379 if zeros is None: 

380 zeros = [] 

381 if poles is None: 

382 poles = [] 

383 if deltat is None: 

384 raise ValueError( 

385 'Sampling interval `deltat` must be given for ' 

386 'DigitalPoleZeroResponse') 

387 

388 FrequencyResponse.__init__( 

389 self, zeros=aslist(zeros), poles=aslist(poles), constant=constant, 

390 deltat=deltat, **kwargs) 

391 

392 def check_sampling_rate(self): 

393 if self.deltat == 0.0: 

394 raise InvalidResponseError( 

395 'Invalid digital response: sampling rate undefined') 

396 

397 def get_fmax(self): 

398 self.check_sampling_rate() 

399 return 0.5 / self.deltat 

400 

401 def evaluate(self, freqs): 

402 self.check_sampling_rate() 

403 return signal.freqz_zpk( 

404 self.zeros, self.poles, self.constant, 

405 freqs*(2.*math.pi*self.deltat))[1] 

406 

407 def is_scalar(self): 

408 return len(self.zeros) == 0 and len(self.poles) == 0 

409 

410 def get_scalar(self): 

411 ''' 

412 Get factor if this is a flat response. 

413 ''' 

414 if self.is_scalar(): 

415 return self.constant 

416 else: 

417 raise IsNotScalar() 

418 

419 def to_digital(self, deltat): 

420 self.check_sampling_rate() 

421 from scipy.signal import zpk2tf 

422 

423 b, a = zpk2tf(self.zeros, self.poles, self.constant) 

424 return DigitalFilterResponse(b, a, deltat) 

425 

426 

427class ButterworthResponse(FrequencyResponse): 

428 ''' 

429 Butterworth frequency response. 

430 

431 :param corner: corner frequency of the response 

432 :param order: order of the response 

433 :param type: either ``high`` or ``low`` 

434 ''' 

435 

436 corner = Float.T(default=1.0) 

437 order = Int.T(default=4) 

438 type = StringChoice.T(choices=['low', 'high'], default='low') 

439 

440 def to_polezero(self): 

441 z, p, k = signal.butter( 

442 self.order, self.corner*2.*math.pi, 

443 btype=self.type, analog=True, output='zpk') 

444 

445 return PoleZeroResponse( 

446 zeros=aslist(z), 

447 poles=aslist(p), 

448 constant=float(k)) 

449 

450 def to_digital(self, deltat): 

451 b, a = signal.butter( 

452 self.order, self.corner*2.*deltat, 

453 self.type, analog=False) 

454 

455 return DigitalFilterResponse(b, a, deltat) 

456 

457 def to_analog(self): 

458 b, a = signal.butter( 

459 self.order, self.corner*2.*math.pi, 

460 self.type, analog=True) 

461 

462 return AnalogFilterResponse(b, a) 

463 

464 def to_digital_polezero(self, deltat): 

465 z, p, k = signal.butter( 

466 self.order, self.corner*2*deltat, 

467 btype=self.type, analog=False, output='zpk') 

468 

469 return DigitalPoleZeroResponse(z, p, k, deltat) 

470 

471 def evaluate(self, freqs): 

472 b, a = signal.butter( 

473 self.order, self.corner*2.*math.pi, 

474 self.type, analog=True) 

475 

476 return signal.freqs(b, a, freqs*2.*math.pi)[1] 

477 

478 

479class SampledResponse(FrequencyResponse): 

480 ''' 

481 Interpolates frequency response given at a set of sampled frequencies. 

482 

483 :param frequencies,values: frequencies and values of the sampled response 

484 function. 

485 :param left,right: values to return when input is out of range. If set to 

486 ``None`` (the default) the endpoints are returned. 

487 ''' 

488 

489 frequencies = Array.T(shape=(None,), dtype=float, serialize_as='list') 

490 values = Array.T(shape=(None,), dtype=complex, serialize_as='list') 

491 left = Complex.T(optional=True) 

492 right = Complex.T(optional=True) 

493 

494 def __init__(self, frequencies, values, left=None, right=None, **kwargs): 

495 FrequencyResponse.__init__( 

496 self, 

497 frequencies=asarray_1d(frequencies, float), 

498 values=asarray_1d(values, complex), 

499 **kwargs) 

500 

501 def evaluate(self, freqs): 

502 ereal = num.interp( 

503 freqs, self.frequencies, num.real(self.values), 

504 left=self.left, right=self.right) 

505 eimag = num.interp( 

506 freqs, self.frequencies, num.imag(self.values), 

507 left=self.left, right=self.right) 

508 transfer = ereal + 1.0j*eimag 

509 return transfer 

510 

511 def inverse(self): 

512 ''' 

513 Get inverse as a new :py:class:`SampledResponse` object. 

514 ''' 

515 

516 def inv_or_none(x): 

517 if x is not None: 

518 return 1./x 

519 

520 return SampledResponse( 

521 self.frequencies, 1./self.values, 

522 left=inv_or_none(self.left), 

523 right=inv_or_none(self.right)) 

524 

525 

526class IntegrationResponse(FrequencyResponse): 

527 ''' 

528 The integration response, optionally multiplied by a constant gain. 

529 

530 :param n: exponent (integer) 

531 :param gain: gain factor (float) 

532 

533 :: 

534 

535 gain 

536 T(f) = -------------- 

537 (j*2*pi * f)^n 

538 ''' 

539 

540 n = Int.T(optional=True, default=1) 

541 gain = Float.T(optional=True, default=1.0) 

542 

543 def __init__(self, n=1, gain=1.0, **kwargs): 

544 FrequencyResponse.__init__(self, n=n, gain=gain, **kwargs) 

545 

546 def evaluate(self, freqs): 

547 nonzero = freqs != 0.0 

548 resp = num.zeros(freqs.size, dtype=complex) 

549 resp[nonzero] = self.gain / (1.0j * 2. * num.pi*freqs[nonzero])**self.n 

550 return resp 

551 

552 

553class DifferentiationResponse(FrequencyResponse): 

554 ''' 

555 The differentiation response, optionally multiplied by a constant gain. 

556 

557 :param n: exponent (integer) 

558 :param gain: gain factor (float) 

559 

560 :: 

561 

562 T(f) = gain * (j*2*pi * f)^n 

563 ''' 

564 

565 n = Int.T(optional=True, default=1) 

566 gain = Float.T(optional=True, default=1.0) 

567 

568 def __init__(self, n=1, gain=1.0, **kwargs): 

569 FrequencyResponse.__init__(self, n=n, gain=gain, **kwargs) 

570 

571 def evaluate(self, freqs): 

572 return self.gain * (1.0j * 2. * num.pi * freqs)**self.n 

573 

574 

575class DigitalFilterResponse(FrequencyResponse): 

576 ''' 

577 Frequency response of an analog filter. 

578 

579 (see :py:func:`scipy.signal.freqz`). 

580 ''' 

581 

582 b = List.T(Float.T()) 

583 a = List.T(Float.T()) 

584 deltat = Float.T() 

585 drop_phase = Bool.T(default=False) 

586 

587 def __init__(self, b, a, deltat, drop_phase=False, **kwargs): 

588 FrequencyResponse.__init__( 

589 self, b=aslist(b), a=aslist(a), deltat=float(deltat), 

590 drop_phase=drop_phase, **kwargs) 

591 

592 def check_sampling_rate(self): 

593 if self.deltat == 0.0: 

594 raise InvalidResponseError( 

595 'Invalid digital response: sampling rate undefined') 

596 

597 def is_scalar(self): 

598 return len(self.a) == 1 and len(self.b) == 1 

599 

600 def get_scalar(self): 

601 if self.is_scalar(): 

602 return self.b[0] / self.a[0] 

603 else: 

604 raise IsNotScalar() 

605 

606 def get_fmax(self): 

607 if not self.is_scalar(): 

608 self.check_sampling_rate() 

609 return 0.5 / self.deltat 

610 else: 

611 return None 

612 

613 def evaluate(self, freqs): 

614 if self.is_scalar(): 

615 return util.num_full_like(freqs, self.get_scalar(), dtype=complex) 

616 

617 self.check_sampling_rate() 

618 

619 ok = freqs <= 0.5/self.deltat 

620 coeffs = num.zeros(freqs.size, dtype=complex) 

621 

622 coeffs[ok] = signal.freqz( 

623 self.b, self.a, freqs[ok]*2.*math.pi * self.deltat)[1] 

624 

625 coeffs[num.logical_not(ok)] = None 

626 if self.drop_phase: 

627 return num.abs(coeffs) 

628 else: 

629 return coeffs 

630 

631 def filter(self, tr): 

632 self.check_sampling_rate() 

633 

634 from pyrocko import trace 

635 trace.assert_same_sampling_rate(self, tr) 

636 tr_new = tr.copy(data=False) 

637 tr_new.set_ydata(signal.lfilter(self.b, self.a, tr.get_ydata())) 

638 return tr_new 

639 

640 

641class AnalogFilterResponse(FrequencyResponse): 

642 ''' 

643 Frequency response of an analog filter. 

644 

645 (see :py:func:`scipy.signal.freqs`). 

646 ''' 

647 

648 b = List.T(Float.T()) 

649 a = List.T(Float.T()) 

650 

651 def __init__(self, b, a, **kwargs): 

652 FrequencyResponse.__init__( 

653 self, b=aslist(b), a=aslist(a), **kwargs) 

654 

655 def evaluate(self, freqs): 

656 return signal.freqs(self.b, self.a, freqs*2.*math.pi)[1] 

657 

658 def to_digital(self, deltat, method='bilinear'): 

659 from scipy.signal import cont2discrete 

660 b, a, _ = cont2discrete((self.b, self.a), deltat, method=method) 

661 if b.ndim == 2: 

662 b = b[0] 

663 return DigitalFilterResponse(b.tolist(), a.tolist(), deltat) 

664 

665 

666class MultiplyResponse(FrequencyResponse): 

667 ''' 

668 Multiplication of several :py:class:`FrequencyResponse` objects. 

669 ''' 

670 

671 responses = List.T(FrequencyResponse.T()) 

672 

673 def __init__(self, responses=None, **kwargs): 

674 if responses is None: 

675 responses = [] 

676 FrequencyResponse.__init__(self, responses=responses, **kwargs) 

677 

678 def get_fmax(self): 

679 fmaxs = [resp.get_fmax() for resp in self.responses] 

680 fmaxs = [fmax for fmax in fmaxs if fmax is not None] 

681 if not fmaxs: 

682 return None 

683 else: 

684 return min(fmaxs) 

685 

686 def evaluate(self, freqs): 

687 a = num.ones(freqs.size, dtype=complex) 

688 for resp in self.responses: 

689 a *= resp.evaluate(freqs) 

690 

691 return a 

692 

693 def is_scalar(self): 

694 return all(resp.is_scalar() for resp in self.responses) 

695 

696 def get_scalar(self): 

697 ''' 

698 Get factor if this is a flat response. 

699 ''' 

700 if self.is_scalar(): 

701 return num.product(resp.get_scalar() for resp in self.responses) 

702 else: 

703 raise IsNotScalar() 

704 

705 def simplify(self): 

706 poles = [] 

707 zeros = [] 

708 constant = 1.0 

709 responses = [] 

710 for resp in self.responses: 

711 if isinstance(resp, PoleZeroResponse): 

712 poles.extend(resp.poles) 

713 zeros.extend(resp.zeros) 

714 constant *= resp.constant 

715 else: 

716 responses.append(resp) 

717 

718 if poles or zeros or constant != 1.0: 

719 responses[0:0] = [ 

720 PoleZeroResponse(poles=poles, zeros=zeros, constant=constant)] 

721 

722 self.responses = responses 

723 

724 def construction(self): 

725 breakpoints = [] 

726 for resp in self.responses: 

727 breakpoints.extend(resp.construction()) 

728 

729 return finalize_construction(breakpoints) 

730 

731 

732class DelayResponse(FrequencyResponse): 

733 

734 delay = Float.T() 

735 

736 def evaluate(self, freqs): 

737 return num.exp(-2.0J * self.delay * num.pi * freqs) 

738 

739 

740class InvalidResponseError(Exception): 

741 pass 

742 

743 

744class InvalidResponse(FrequencyResponse): 

745 

746 message = String.T() 

747 

748 def __init__(self, message): 

749 FrequencyResponse.__init__(self, message=message) 

750 self.have_warned = False 

751 

752 def evaluate(self, freqs): 

753 if not self.have_warned: 

754 logger.warning('Invalid response: %s' % self.message) 

755 self.have_warned = True 

756 

757 return util.num_full_like(freqs, None, dtype=num.complex) 

758 

759 

760def simplify_responses(responses): 

761 

762 def unpack_multi(responses): 

763 for resp in responses: 

764 if isinstance(resp, MultiplyResponse): 

765 for sub in unpack_multi(resp.responses): 

766 yield sub 

767 else: 

768 yield resp 

769 

770 def cancel_pzs(poles, zeros): 

771 poles_new = [] 

772 zeros_new = list(zeros) 

773 for p in poles: 

774 try: 

775 zeros_new.pop(zeros_new.index(p)) 

776 except ValueError: 

777 poles_new.append(p) 

778 

779 return poles_new, zeros_new 

780 

781 def combine_pzs(responses): 

782 poles = [] 

783 zeros = [] 

784 constant = 1.0 

785 out = [] 

786 for resp in responses: 

787 if isinstance(resp, PoleZeroResponse): 

788 poles.extend(resp.poles) 

789 zeros.extend(resp.zeros) 

790 constant *= resp.constant 

791 else: 

792 out.append(resp) 

793 

794 poles, zeros = cancel_pzs(poles, zeros) 

795 if poles or zeros: 

796 out.insert(0, PoleZeroResponse( 

797 poles=poles, zeros=zeros, constant=constant)) 

798 elif constant != 1.0: 

799 out.insert(0, Gain(constant=constant)) 

800 

801 return out 

802 

803 def split(xs, condition): 

804 out = [], [] 

805 for x in xs: 

806 out[condition(x)].append(x) 

807 

808 return out 

809 

810 def combine_gains(responses): 

811 non_scalars, scalars = split(responses, lambda resp: resp.is_scalar()) 

812 if scalars: 

813 factor = num.prod([resp.get_scalar() for resp in scalars]) 

814 yield Gain(constant=factor) 

815 

816 for resp in non_scalars: 

817 yield resp 

818 

819 return list(combine_gains(combine_pzs(unpack_multi(responses))))