Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/util.py: 78%

1229 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-04 09:52 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5''' 

6Utility functions for Pyrocko. 

7 

8.. _time-handling-mode: 

9 

10High precision time handling mode 

11................................. 

12 

13Pyrocko can treat timestamps either as standard double precision (64 bit) 

14floating point values, or as high precision floats (``numpy.float128`` or 

15``numpy.float96``, whichever is available, `see NumPy Scalars 

16<https://numpy.org/doc/stable/reference/arrays.scalars.html>`_), aliased here 

17as :py:class:`~pyrocko.util.hpfloat`. High precision time stamps are required 

18when handling data with sub-millisecond precision, i.e. kHz/MHz data streams 

19and event catalogs derived from such data. 

20 

21Not all functions in Pyrocko and in programs depending on Pyrocko may work 

22correctly with high precision times. Therefore, Pyrocko's high precision time 

23handling mode has to be actively activated by user config, command line option 

24or enforced within a certain script/program. 

25 

26The default high precision time handling mode can be configured globally with 

27the user configuration variable 

28:py:gattr:`~pyrocko.config.PyrockoConfig.use_high_precision_time`. Calling the 

29function :py:func:`use_high_precision_time` overrides the default from the 

30config file. This function may be called at startup of a program/script 

31requiring a specific time handling mode. 

32 

33To create a valid time stamp for use in Pyrocko (e.g. in 

34:py:class:`~pyrocko.model.event.Event` or 

35:py:class:`~pyrocko.trace.Trace` objects), use: 

36 

37.. code-block :: python 

38 

39 import time 

40 from pyrocko import util 

41 

42 # By default using mode selected in user config, override with: 

43 # util.use_high_precision_time(True) # force high precision mode 

44 # util.use_high_precision_time(False) # force low precision mode 

45 

46 t1 = util.str_to_time('2020-08-27 10:22:00') 

47 t2 = util.str_to_time('2020-08-27 10:22:00.111222') 

48 t3 = util.to_time_float(time.time()) 

49 

50 # To get the appropriate float class, use: 

51 

52 time_float = util.get_time_float() 

53 # -> float, numpy.float128 or numpy.float96 

54 [isinstance(t, time_float) for t in [t1, t2, t3]] 

55 # -> [True, True, True] 

56 

57 # Shortcut: 

58 util.check_time_class(t1) 

59 

60Module content 

61.............. 

62 

63.. py:class:: hpfloat 

64 

65 Alias for NumPy's high precision float data type ``float128`` or 

66 ``float96``, if available. 

67 

68 On platforms lacking support for high precision floats, an attempt to 

69 create a ``hpfloat`` instance, raises :py:exc:`HPFloatUnavailable`. 

70 

71''' 

72 

73import time 

74import logging 

75import os 

76import sys 

77import re 

78import calendar 

79import math 

80import fnmatch 

81import inspect 

82import weakref 

83try: 

84 import fcntl 

85except ImportError: 

86 fcntl = None 

87import optparse 

88import os.path as op 

89import errno 

90 

91import numpy as num 

92from scipy import signal 

93import pyrocko 

94from pyrocko import dummy_progressbar 

95 

96 

97from urllib.parse import urlencode, quote, unquote # noqa 

98from urllib.request import Request, build_opener, HTTPDigestAuthHandler # noqa 

99from urllib.request import urlopen as _urlopen # noqa 

100from urllib.error import HTTPError, URLError # noqa 

101 

102 

103try: 

104 import certifi 

105 import ssl 

106 g_ssl_context = ssl.create_default_context(cafile=certifi.where()) 

107except ImportError: 

108 g_ssl_context = None 

109 

110 

111class URLErrorSSL(URLError): 

112 

113 def __init__(self, urlerror): 

114 self.__dict__ = urlerror.__dict__.copy() 

115 

116 def __str__(self): 

117 return ( 

118 'Requesting web resource failed and the problem could be ' 

119 'related to SSL. Python standard libraries on some older ' 

120 'systems (like Ubuntu 14.04) are known to have trouble ' 

121 "with some SSL setups of today's servers: %s" 

122 % URLError.__str__(self)) 

123 

124 

125def urlopen_ssl_check(*args, **kwargs): 

126 try: 

127 return urlopen(*args, **kwargs) 

128 except URLError as e: 

129 if str(e).find('SSL') != -1: 

130 raise URLErrorSSL(e) 

131 else: 

132 raise 

133 

134 

135def urlopen(*args, **kwargs): 

136 

137 if 'context' not in kwargs and g_ssl_context is not None: 

138 kwargs['context'] = g_ssl_context 

139 

140 return _urlopen(*args, **kwargs) 

141 

142 

143try: 

144 long 

145except NameError: 

146 long = int 

147 

148 

149force_dummy_progressbar = False 

150 

151 

152try: 

153 from pyrocko import util_ext 

154except ImportError: 

155 util_ext = None 

156 

157 

158logger = logging.getLogger('pyrocko.util') 

159 

160 

161try: 

162 num_full = num.full 

163except AttributeError: 

164 def num_full(shape, fill_value, dtype=None, order='C'): 

165 a = num.empty(shape, dtype=dtype, order=order) 

166 a.fill(fill_value) 

167 return a 

168 

169try: 

170 num_full_like = num.full_like 

171except AttributeError: 

172 def num_full_like(arr, fill_value, dtype=None, order='K', subok=True): 

173 a = num.empty_like(arr, dtype=dtype, order=order, subok=subok) 

174 a.fill(fill_value) 

175 return a 

176 

177 

178g_setup_logging_args = 'pyrocko', 'warning' 

179 

180 

181def setup_logging(programname='pyrocko', levelname='warning'): 

182 ''' 

183 Initialize logging. 

184 

185 :param programname: program name to be written in log 

186 :param levelname: string indicating the logging level ('debug', 'info', 

187 'warning', 'error', 'critical') 

188 

189 This function is called at startup by most pyrocko programs to set up a 

190 consistent logging format. This is simply a shortcut to a call to 

191 :py:func:`logging.basicConfig()`. 

192 ''' 

193 

194 global g_setup_logging_args 

195 g_setup_logging_args = (programname, levelname) 

196 

197 levels = {'debug': logging.DEBUG, 

198 'info': logging.INFO, 

199 'warning': logging.WARNING, 

200 'error': logging.ERROR, 

201 'critical': logging.CRITICAL} 

202 

203 logging.basicConfig( 

204 level=levels[levelname], 

205 format=programname+':%(name)-25s - %(levelname)-8s - %(message)s') 

206 

207 

208def subprocess_setup_logging_args(): 

209 ''' 

210 Get arguments from previous call to setup_logging. 

211 

212 These can be sent down to a worker process so it can setup its logging 

213 in the same way as the main process. 

214 ''' 

215 return g_setup_logging_args 

216 

217 

218def data_file(fn): 

219 return os.path.join(os.path.split(__file__)[0], 'data', fn) 

220 

221 

222class DownloadError(Exception): 

223 ''' 

224 Raised when a download failed. 

225 ''' 

226 pass 

227 

228 

229class PathExists(DownloadError): 

230 ''' 

231 Raised when the download target file already exists. 

232 ''' 

233 pass 

234 

235 

236def _download(url, fpath, username=None, password=None, 

237 force=False, method='download', stats=None, 

238 status_callback=None, entries_wanted=None, 

239 recursive=False, header=None): 

240 

241 import requests 

242 from requests.auth import HTTPBasicAuth 

243 from requests.exceptions import HTTPError as req_HTTPError 

244 

245 requests.adapters.DEFAULT_RETRIES = 5 

246 urljoin = requests.compat.urljoin 

247 

248 session = requests.Session() 

249 session.header = header 

250 session.auth = None if username is None\ 

251 else HTTPBasicAuth(username, password) 

252 

253 status = { 

254 'ntotal_files': 0, 

255 'nread_files': 0, 

256 'ntotal_bytes_all_files': 0, 

257 'nread_bytes_all_files': 0, 

258 'ntotal_bytes_current_file': 0, 

259 'nread_bytes_current_file': 0, 

260 'finished': False 

261 } 

262 

263 try: 

264 url_to_size = {} 

265 

266 if callable(status_callback): 

267 status_callback(status) 

268 

269 if not recursive and url.endswith('/'): 

270 raise DownloadError( 

271 'URL: %s appears to be a directory' 

272 ' but recurvise download is False' % url) 

273 

274 if recursive and not url.endswith('/'): 

275 url += '/' 

276 

277 def parse_directory_tree(url, subdir=''): 

278 r = session.get(urljoin(url, subdir)) 

279 r.raise_for_status() 

280 

281 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text) 

282 

283 files = sorted(set(subdir + fn for fn in entries 

284 if not fn.endswith('/'))) 

285 

286 if entries_wanted is not None: 

287 files = [fn for fn in files 

288 if (fn in wanted for wanted in entries_wanted)] 

289 

290 status['ntotal_files'] += len(files) 

291 

292 dirs = sorted(set(subdir + dn for dn in entries 

293 if dn.endswith('/') 

294 and dn not in ('./', '../'))) 

295 

296 for dn in dirs: 

297 files.extend(parse_directory_tree( 

298 url, subdir=dn)) 

299 

300 return files 

301 

302 def get_content_length(url): 

303 if url not in url_to_size: 

304 r = session.head(url, headers={'Accept-Encoding': ''}) 

305 

306 content_length = r.headers.get('content-length', None) 

307 if content_length is None: 

308 logger.debug('Could not get HTTP header ' 

309 'Content-Length for %s' % url) 

310 

311 content_length = None 

312 

313 else: 

314 content_length = int(content_length) 

315 status['ntotal_bytes_all_files'] += content_length 

316 if callable(status_callback): 

317 status_callback(status) 

318 

319 url_to_size[url] = content_length 

320 

321 return url_to_size[url] 

322 

323 def download_file(url, fn): 

324 logger.info('starting download of %s...' % url) 

325 ensuredirs(fn) 

326 

327 fsize = get_content_length(url) 

328 status['ntotal_bytes_current_file'] = fsize 

329 status['nread_bytes_current_file'] = 0 

330 if callable(status_callback): 

331 status_callback(status) 

332 

333 r = session.get(url, stream=True, timeout=5) 

334 r.raise_for_status() 

335 

336 frx = 0 

337 fn_tmp = fn + '.%i.temp' % os.getpid() 

338 with open(fn_tmp, 'wb') as f: 

339 for d in r.iter_content(chunk_size=1024): 

340 f.write(d) 

341 frx += len(d) 

342 

343 status['nread_bytes_all_files'] += len(d) 

344 status['nread_bytes_current_file'] += len(d) 

345 if callable(status_callback): 

346 status_callback(status) 

347 

348 os.rename(fn_tmp, fn) 

349 

350 if fsize is not None and frx != fsize: 

351 logger.warning( 

352 'HTTP header Content-Length: %i bytes does not match ' 

353 'download size %i bytes' % (fsize, frx)) 

354 

355 logger.info('finished download of %s' % url) 

356 

357 status['nread_files'] += 1 

358 if callable(status_callback): 

359 status_callback(status) 

360 

361 if recursive: 

362 if op.exists(fpath) and not force: 

363 raise PathExists('path %s already exists' % fpath) 

364 

365 files = parse_directory_tree(url) 

366 

367 dsize = 0 

368 for fn in files: 

369 file_url = urljoin(url, fn) 

370 dsize += get_content_length(file_url) or 0 

371 

372 if method == 'calcsize': 

373 return dsize 

374 

375 else: 

376 for fn in files: 

377 file_url = urljoin(url, fn) 

378 download_file(file_url, op.join(fpath, fn)) 

379 

380 else: 

381 status['ntotal_files'] += 1 

382 if callable(status_callback): 

383 status_callback(status) 

384 

385 fsize = get_content_length(url) 

386 if method == 'calcsize': 

387 return fsize 

388 else: 

389 download_file(url, fpath) 

390 

391 except req_HTTPError as e: 

392 logging.warning('http error: %s' % e) 

393 raise DownloadError('could not download file(s) from: %s' % url) 

394 

395 finally: 

396 status['finished'] = True 

397 if callable(status_callback): 

398 status_callback(status) 

399 session.close() 

400 

401 

402def download_file( 

403 url, fpath, username=None, password=None, status_callback=None, 

404 **kwargs): 

405 return _download( 

406 url, fpath, username, password, 

407 recursive=False, 

408 status_callback=status_callback, 

409 **kwargs) 

410 

411 

412def download_dir( 

413 url, fpath, username=None, password=None, status_callback=None, 

414 **kwargs): 

415 

416 return _download( 

417 url, fpath, username, password, 

418 recursive=True, 

419 status_callback=status_callback, 

420 **kwargs) 

421 

422 

423class HPFloatUnavailable(Exception): 

424 ''' 

425 Raised when a high precision float type would be required but is not 

426 available. 

427 ''' 

428 pass 

429 

430 

431class dummy_hpfloat(object): 

432 def __init__(self, *args, **kwargs): 

433 raise HPFloatUnavailable( 

434 'NumPy lacks support for float128 or float96 data type on this ' 

435 'platform.') 

436 

437 

438if hasattr(num, 'float128'): 

439 hpfloat = num.float128 

440 

441elif hasattr(num, 'float96'): 

442 hpfloat = num.float96 

443 

444else: 

445 hpfloat = dummy_hpfloat 

446 

447 

448g_time_float = None 

449g_time_dtype = None 

450 

451 

452class TimeFloatSettingError(Exception): 

453 pass 

454 

455 

456def use_high_precision_time(enabled): 

457 ''' 

458 Globally force a specific time handling mode. 

459 

460 See :ref:`High precision time handling mode <time-handling-mode>`. 

461 

462 :param enabled: enable/disable use of high precision time type 

463 :type enabled: bool 

464 

465 This function should be called before handling/reading any time data. 

466 It can only be called once. 

467 

468 Special attention is required when using multiprocessing on a platform 

469 which does not use fork under the hood. In such cases, the desired setting 

470 must be set also in the subprocess. 

471 ''' 

472 _setup_high_precision_time_mode(enabled_app=enabled) 

473 

474 

475def _setup_high_precision_time_mode(enabled_app=False): 

476 global g_time_float 

477 global g_time_dtype 

478 

479 if not (g_time_float is None and g_time_dtype is None): 

480 raise TimeFloatSettingError( 

481 'Cannot set time handling mode: too late, it has already been ' 

482 'fixed by an earlier call.') 

483 

484 from pyrocko import config 

485 

486 conf = config.config() 

487 enabled_config = conf.use_high_precision_time 

488 

489 enabled_env = os.environ.get('PYROCKO_USE_HIGH_PRECISION_TIME', None) 

490 if enabled_env is not None: 

491 try: 

492 enabled_env = int(enabled_env) == 1 

493 except ValueError: 

494 raise TimeFloatSettingError( 

495 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME ' 

496 'should be set to 0 or 1.') 

497 

498 enabled = enabled_config 

499 mode_from = 'config variable `use_high_precision_time`' 

500 notify = enabled 

501 

502 if enabled_env is not None: 

503 if enabled_env != enabled: 

504 notify = True 

505 enabled = enabled_env 

506 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`' 

507 

508 if enabled_app is not None: 

509 if enabled_app != enabled: 

510 notify = True 

511 enabled = enabled_app 

512 mode_from = 'application override' 

513 

514 logger.debug(''' 

515Pyrocko high precision time mode selection (latter override earlier): 

516 config: %s 

517 env: %s 

518 app: %s 

519 -> enabled: %s'''.lstrip() % ( 

520 enabled_config, enabled_env, enabled_app, enabled)) 

521 

522 if notify: 

523 logger.info('Pyrocko high precision time mode %s by %s.' % ( 

524 'activated' if enabled else 'deactivated', 

525 mode_from)) 

526 

527 if enabled: 

528 g_time_float = hpfloat 

529 g_time_dtype = hpfloat 

530 else: 

531 g_time_float = float 

532 g_time_dtype = num.float64 

533 

534 

535def get_time_float(): 

536 ''' 

537 Get the effective float class for timestamps. 

538 

539 See :ref:`High precision time handling mode <time-handling-mode>`. 

540 

541 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the 

542 current time handling mode 

543 ''' 

544 global g_time_float 

545 

546 if g_time_float is None: 

547 _setup_high_precision_time_mode() 

548 

549 return g_time_float 

550 

551 

552def get_time_dtype(): 

553 ''' 

554 Get effective NumPy float class to handle timestamps. 

555 

556 See :ref:`High precision time handling mode <time-handling-mode>`. 

557 ''' 

558 

559 global g_time_dtype 

560 

561 if g_time_dtype is None: 

562 _setup_high_precision_time_mode() 

563 

564 return g_time_dtype 

565 

566 

567def to_time_float(t): 

568 ''' 

569 Convert float to valid time stamp in the current time handling mode. 

570 

571 See :ref:`High precision time handling mode <time-handling-mode>`. 

572 ''' 

573 return get_time_float()(t) 

574 

575 

576class TimestampTypeError(ValueError): 

577 pass 

578 

579 

580def check_time_class(t, error='raise'): 

581 ''' 

582 Type-check variable against current time handling mode. 

583 

584 See :ref:`High precision time handling mode <time-handling-mode>`. 

585 ''' 

586 

587 if t == 0.0: 

588 return 

589 

590 if not isinstance(t, get_time_float()): 

591 message = ( 

592 'Timestamp %g is of type %s but should be of type %s with ' 

593 "Pyrocko's currently selected time handling mode.\n\n" 

594 'See https://pyrocko.org/docs/current/library/reference/util.html' 

595 '#high-precision-time-handling-mode' % ( 

596 t, type(t), get_time_float())) 

597 

598 if error == 'raise': 

599 raise TimestampTypeError(message) 

600 elif error == 'warn': 

601 logger.warning(message) 

602 else: 

603 assert False 

604 

605 

606class Stopwatch(object): 

607 ''' 

608 Simple stopwatch to measure elapsed wall clock time. 

609 

610 Usage:: 

611 

612 s = Stopwatch() 

613 time.sleep(1) 

614 print s() 

615 time.sleep(1) 

616 print s() 

617 ''' 

618 

619 def __init__(self): 

620 self.start = time.time() 

621 

622 def __call__(self): 

623 return time.time() - self.start 

624 

625 

626def wrap(text, line_length=80): 

627 ''' 

628 Paragraph and list-aware wrapping of text. 

629 ''' 

630 

631 text = text.strip('\n') 

632 at_lineend = re.compile(r' *\n') 

633 at_para = re.compile(r'((^|(\n\s*)?\n)(\s+[*] )|\n\s*\n)') 

634 

635 paragraphs = at_para.split(text)[::5] 

636 listindents = at_para.split(text)[4::5] 

637 newlist = at_para.split(text)[3::5] 

638 

639 listindents[0:0] = [None] 

640 listindents.append(True) 

641 newlist.append(None) 

642 

643 det_indent = re.compile(r'^ *') 

644 

645 outlines = [] 

646 for ip, p in enumerate(paragraphs): 

647 if not p: 

648 continue 

649 

650 if listindents[ip] is None: 

651 _indent = det_indent.findall(p)[0] 

652 findent = _indent 

653 else: 

654 findent = listindents[ip] 

655 _indent = ' ' * len(findent) 

656 

657 ll = line_length - len(_indent) 

658 llf = ll 

659 

660 oldlines = [s.strip() for s in at_lineend.split(p.rstrip())] 

661 p1 = ' '.join(oldlines) 

662 possible = re.compile(r'(^.{1,%i}|.{1,%i})( |$)' % (llf, ll)) 

663 for imatch, match in enumerate(possible.finditer(p1)): 

664 parout = match.group(1) 

665 if imatch == 0: 

666 outlines.append(findent + parout) 

667 else: 

668 outlines.append(_indent + parout) 

669 

670 if ip != len(paragraphs)-1 and ( 

671 listindents[ip] is None 

672 or newlist[ip] is not None 

673 or listindents[ip+1] is None): 

674 

675 outlines.append('') 

676 

677 return outlines 

678 

679 

680def ewrap(lines, width=80, indent=''): 

681 lines = list(lines) 

682 if not lines: 

683 return '' 

684 fwidth = max(len(s) for s in lines) 

685 nx = max(1, (80-len(indent)) // (fwidth+1)) 

686 i = 0 

687 rows = [] 

688 while i < len(lines): 

689 rows.append(indent + ' '.join(x.ljust(fwidth) for x in lines[i:i+nx])) 

690 i += nx 

691 

692 return '\n'.join(rows) 

693 

694 

695class BetterHelpFormatter(optparse.IndentedHelpFormatter): 

696 

697 def __init__(self, *args, **kwargs): 

698 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs) 

699 

700 def format_option(self, option): 

701 ''' 

702 From IndentedHelpFormatter but using a different wrap method. 

703 ''' 

704 

705 help_text_position = 4 + self.current_indent 

706 help_text_width = self.width - help_text_position 

707 

708 opts = self.option_strings[option] 

709 opts = '%*s%s' % (self.current_indent, '', opts) 

710 if option.help: 

711 help_text = self.expand_default(option) 

712 

713 if self.help_position + len(help_text) + 1 <= self.width: 

714 lines = [ 

715 '', 

716 '%-*s %s' % (self.help_position, opts, help_text), 

717 ''] 

718 else: 

719 lines = [''] 

720 lines.append(opts) 

721 lines.append('') 

722 if option.help: 

723 help_lines = wrap(help_text, help_text_width) 

724 lines.extend(['%*s%s' % (help_text_position, '', line) 

725 for line in help_lines]) 

726 lines.append('') 

727 

728 return '\n'.join(lines) 

729 

730 def format_description(self, description): 

731 if not description: 

732 return '' 

733 

734 if self.current_indent == 0: 

735 lines = [] 

736 else: 

737 lines = [''] 

738 

739 lines.extend(wrap(description, self.width - self.current_indent)) 

740 if self.current_indent == 0: 

741 lines.append('\n') 

742 

743 return '\n'.join( 

744 ['%*s%s' % (self.current_indent, '', line) for line in lines]) 

745 

746 

747class ProgressBar: 

748 def __init__(self, label, n): 

749 from pyrocko.progress import progress 

750 self._context = progress.view() 

751 self._context.__enter__() 

752 self._task = progress.task(label, n) 

753 

754 def update(self, i): 

755 self._task.update(i) 

756 

757 def finish(self): 

758 self._task.done() 

759 if self._context: 

760 self._context.__exit__() 

761 self._context = None 

762 

763 

764def progressbar(label, maxval): 

765 if force_dummy_progressbar: 

766 return dummy_progressbar.ProgressBar(maxval=maxval).start() 

767 

768 return ProgressBar(label, maxval) 

769 

770 

771def progress_beg(label): 

772 ''' 

773 Notify user that an operation has started. 

774 

775 :param label: name of the operation 

776 

777 To be used in conjuction with :py:func:`progress_end`. 

778 ''' 

779 

780 sys.stderr.write(label) 

781 sys.stderr.flush() 

782 

783 

784def progress_end(label=''): 

785 ''' 

786 Notify user that an operation has ended. 

787 

788 :param label: name of the operation 

789 

790 To be used in conjuction with :py:func:`progress_beg`. 

791 ''' 

792 

793 sys.stderr.write(' done. %s\n' % label) 

794 sys.stderr.flush() 

795 

796 

797class ArangeError(ValueError): 

798 ''' 

799 Raised by :py:func:`arange2` for inconsistent range specifications. 

800 ''' 

801 pass 

802 

803 

804def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'): 

805 ''' 

806 Return evenly spaced numbers over a specified interval. 

807 

808 Like :py:func:`numpy.arange` but returning floating point numbers by 

809 default and with defined behaviour when stepsize is inconsistent with 

810 interval bounds. It is considered inconsistent if the difference between 

811 the closest multiple of ``step`` and ``stop`` is larger than ``epsilon * 

812 step``. Inconsistencies are handled according to the ``error`` parameter. 

813 If it is set to ``'raise'`` an exception of type :py:exc:`ArangeError` is 

814 raised. If it is set to ``'round'``, ``'floor'``, or ``'ceil'``, ``stop`` 

815 is silently changed to the closest, the next smaller, or next larger 

816 multiple of ``step``, respectively. 

817 ''' 

818 

819 assert error in ('raise', 'round', 'floor', 'ceil') 

820 

821 start = dtype(start) 

822 stop = dtype(stop) 

823 step = dtype(step) 

824 

825 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round) 

826 

827 n = int(rnd((stop - start) / step)) + 1 

828 stop_check = start + (n-1) * step 

829 

830 if error == 'raise' and abs(stop_check - stop) > step * epsilon: 

831 raise ArangeError( 

832 'inconsistent range specification: start=%g, stop=%g, step=%g' 

833 % (start, stop, step)) 

834 

835 x = num.arange(n, dtype=dtype) 

836 x *= step 

837 x += start 

838 return x 

839 

840 

841def polylinefit(x, y, n_or_xnodes): 

842 ''' 

843 Fit piece-wise linear function to data. 

844 

845 :param x,y: arrays with coordinates of data 

846 :param n_or_xnodes: int, number of segments or x coordinates of polyline 

847 

848 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of 

849 polyline, root-mean-square error 

850 ''' 

851 

852 x = num.asarray(x) 

853 y = num.asarray(y) 

854 

855 if isinstance(n_or_xnodes, int): 

856 n = n_or_xnodes 

857 xmin = x.min() 

858 xmax = x.max() 

859 xnodes = num.linspace(xmin, xmax, n+1) 

860 else: 

861 xnodes = num.asarray(n_or_xnodes) 

862 n = xnodes.size - 1 

863 

864 assert len(x) == len(y) 

865 assert n > 0 

866 

867 ndata = len(x) 

868 a = num.zeros((ndata+(n-1), n*2)) 

869 for i in range(n): 

870 xmin_block = xnodes[i] 

871 xmax_block = xnodes[i+1] 

872 if i == n-1: # don't loose last point 

873 indices = num.where( 

874 num.logical_and(xmin_block <= x, x <= xmax_block))[0] 

875 else: 

876 indices = num.where( 

877 num.logical_and(xmin_block <= x, x < xmax_block))[0] 

878 

879 a[indices, i*2] = x[indices] 

880 a[indices, i*2+1] = 1.0 

881 

882 w = float(ndata)*100. 

883 if i < n-1: 

884 a[ndata+i, i*2] = xmax_block*w 

885 a[ndata+i, i*2+1] = 1.0*w 

886 a[ndata+i, i*2+2] = -xmax_block*w 

887 a[ndata+i, i*2+3] = -1.0*w 

888 

889 d = num.concatenate((y, num.zeros(n-1))) 

890 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2)) 

891 

892 ynodes = num.zeros(n+1) 

893 ynodes[:n] = model[:, 0]*xnodes[:n] + model[:, 1] 

894 ynodes[1:] += model[:, 0]*xnodes[1:] + model[:, 1] 

895 ynodes[1:n] *= 0.5 

896 

897 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2)) 

898 

899 return xnodes, ynodes, rms_error 

900 

901 

902def plf_integrate_piecewise(x_edges, x, y): 

903 ''' 

904 Calculate definite integral of piece-wise linear function on intervals. 

905 

906 Use trapezoidal rule to calculate definite integral of a piece-wise linear 

907 function for a series of consecutive intervals. ``x_edges`` and ``x`` must 

908 be sorted. 

909 

910 :param x_edges: array with edges of the intervals 

911 :param x,y: arrays with coordinates of piece-wise linear function's 

912 control points 

913 ''' 

914 

915 x_all = num.concatenate((x, x_edges)) 

916 ii = num.argsort(x_all) 

917 y_edges = num.interp(x_edges, x, y) 

918 y_all = num.concatenate((y, y_edges)) 

919 xs = x_all[ii] 

920 ys = y_all[ii] 

921 y_all[ii[1:]] = num.cumsum((xs[1:] - xs[:-1]) * 0.5 * (ys[1:] + ys[:-1])) 

922 return num.diff(y_all[-len(y_edges):]) 

923 

924 

925def diff_fd_1d_4o(dt, data): 

926 ''' 

927 Approximate first derivative of an array (forth order, central FD). 

928 

929 :param dt: sampling interval 

930 :param data: NumPy array with data samples 

931 

932 :returns: NumPy array with same shape as input 

933 

934 Interior points are approximated to fourth order, edge points to first 

935 order right- or left-sided respectively, points next to edge to second 

936 order central. 

937 ''' 

938 import scipy.signal 

939 

940 ddata = num.empty_like(data, dtype=float) 

941 

942 if data.size >= 5: 

943 ddata[2:-2] = scipy.signal.lfilter( 

944 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt) 

945 

946 if data.size >= 3: 

947 ddata[1] = (data[2] - data[0]) / (2. * dt) 

948 ddata[-2] = (data[-1] - data[-3]) / (2. * dt) 

949 

950 if data.size >= 2: 

951 ddata[0] = (data[1] - data[0]) / dt 

952 ddata[-1] = (data[-1] - data[-2]) / dt 

953 

954 if data.size == 1: 

955 ddata[0] = 0.0 

956 

957 return ddata 

958 

959 

960def diff_fd_1d_2o(dt, data): 

961 ''' 

962 Approximate first derivative of an array (second order, central FD). 

963 

964 :param dt: sampling interval 

965 :param data: NumPy array with data samples 

966 

967 :returns: NumPy array with same shape as input 

968 

969 Interior points are approximated to second order, edge points to first 

970 order right- or left-sided respectively. 

971 

972 Uses :py:func:`numpy.gradient`. 

973 ''' 

974 

975 return num.gradient(data, dt) 

976 

977 

978def diff_fd_2d_4o(dt, data): 

979 ''' 

980 Approximate second derivative of an array (forth order, central FD). 

981 

982 :param dt: sampling interval 

983 :param data: NumPy array with data samples 

984 

985 :returns: NumPy array with same shape as input 

986 

987 Interior points are approximated to fourth order, next-to-edge points to 

988 second order, edge points repeated. 

989 ''' 

990 import scipy.signal 

991 

992 ddata = num.empty_like(data, dtype=float) 

993 

994 if data.size >= 5: 

995 ddata[2:-2] = scipy.signal.lfilter( 

996 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2) 

997 

998 if data.size >= 3: 

999 ddata[:2] = (data[2] - 2.0 * data[1] + data[0]) / dt**2 

1000 ddata[-2:] = (data[-1] - 2.0 * data[-2] + data[-3]) / dt**2 

1001 

1002 if data.size < 3: 

1003 ddata[:] = 0.0 

1004 

1005 return ddata 

1006 

1007 

1008def diff_fd_2d_2o(dt, data): 

1009 ''' 

1010 Approximate second derivative of an array (second order, central FD). 

1011 

1012 :param dt: sampling interval 

1013 :param data: NumPy array with data samples 

1014 

1015 :returns: NumPy array with same shape as input 

1016 

1017 Interior points are approximated to second order, edge points repeated. 

1018 ''' 

1019 import scipy.signal 

1020 

1021 ddata = num.empty_like(data, dtype=float) 

1022 

1023 if data.size >= 3: 

1024 ddata[1:-1] = scipy.signal.lfilter( 

1025 [1., -2., 1.], [1.], data)[2:] / (dt**2) 

1026 

1027 ddata[0] = ddata[1] 

1028 ddata[-1] = ddata[-2] 

1029 

1030 if data.size < 3: 

1031 ddata[:] = 0.0 

1032 

1033 return ddata 

1034 

1035 

1036def diff_fd(n, order, dt, data): 

1037 ''' 

1038 Approximate 1st or 2nd derivative of an array. 

1039 

1040 :param n: 1 for first derivative, 2 for second 

1041 :param order: order of the approximation 2 and 4 are supported 

1042 :param dt: sampling interval 

1043 :param data: NumPy array with data samples 

1044 

1045 :returns: NumPy array with same shape as input 

1046 

1047 This is a frontend to the functions :py:func:`diff_fd_1d_2o`, 

1048 :py:func:`diff_fd_1d_4o`, :py:func:`diff_fd_2d_2o`, and 

1049 :py:func:`diff_fd_2d_4o`. 

1050 

1051 Raises :py:exc:`ValueError` for unsupported `n` or `order`. 

1052 ''' 

1053 

1054 funcs = { 

1055 1: {2: diff_fd_1d_2o, 4: diff_fd_1d_4o}, 

1056 2: {2: diff_fd_2d_2o, 4: diff_fd_2d_4o}} 

1057 

1058 try: 

1059 funcs_n = funcs[n] 

1060 except KeyError: 

1061 raise ValueError( 

1062 'pyrocko.util.diff_fd: ' 

1063 'Only 1st and 2sd derivatives are supported.') 

1064 

1065 try: 

1066 func = funcs_n[order] 

1067 except KeyError: 

1068 raise ValueError( 

1069 'pyrocko.util.diff_fd: ' 

1070 'Order %i is not supported for %s derivative. Supported: %s' % ( 

1071 order, ['', '1st', '2nd'][n], 

1072 ', '.join('%i' % order for order in sorted(funcs_n.keys())))) 

1073 

1074 return func(dt, data) 

1075 

1076 

1077class GlobalVars(object): 

1078 reuse_store = dict() 

1079 decitab_nmax = 0 

1080 decitab = {} 

1081 decimate_fir_coeffs = {} 

1082 decimate_fir_remez_coeffs = {} 

1083 decimate_iir_coeffs = {} 

1084 re_frac = None 

1085 

1086 

1087def decimate_coeffs(q, n=None, ftype='iir'): 

1088 

1089 q = int(q) 

1090 

1091 if n is None: 

1092 if ftype == 'fir': 

1093 n = 30 

1094 elif ftype == 'fir-remez': 

1095 n = 45*q 

1096 else: 

1097 n = 8 

1098 

1099 if ftype == 'fir': 

1100 coeffs = GlobalVars.decimate_fir_coeffs 

1101 if (n, 1./q) not in coeffs: 

1102 coeffs[n, 1./q] = signal.firwin(n+1, .75/q, window='hamming') 

1103 

1104 b = coeffs[n, 1./q] 

1105 return b, [1.], n 

1106 

1107 elif ftype == 'fir-remez': 

1108 coeffs = GlobalVars.decimate_fir_remez_coeffs 

1109 if (n, 1./q) not in coeffs: 

1110 coeffs[n, 1./q] = signal.remez( 

1111 n+1, (0., .75/q, 1./q, 1.), 

1112 (1., 0.), fs=2, weight=(1, 50)) 

1113 b = coeffs[n, 1./q] 

1114 return b, [1.], n 

1115 

1116 else: 

1117 coeffs = GlobalVars.decimate_iir_coeffs 

1118 if (n, 0.05, 0.8/q) not in coeffs: 

1119 coeffs[n, 0.05, 0.8/q] = signal.cheby1(n, 0.05, 0.8/q) 

1120 

1121 b, a = coeffs[n, 0.05, 0.8/q] 

1122 return b, a, n 

1123 

1124 

1125def decimate(x, q, n=None, ftype='iir', zi=None, ioff=0): 

1126 ''' 

1127 Downsample the signal x by an integer factor q, using an order n filter 

1128 

1129 By default, an order 8 Chebyshev type I filter is used or a 30 point FIR 

1130 filter with hamming window if ftype is 'fir'. 

1131 

1132 :param x: the signal to be downsampled (1D :class:`numpy.ndarray`) 

1133 :param q: the downsampling factor 

1134 :param n: order of the filter (1 less than the length of the filter for a 

1135 `fir` filter) 

1136 :param ftype: type of the filter; can be `iir`, `fir` or `fir-remez` 

1137 

1138 :returns: the downsampled signal (1D :class:`numpy.ndarray`) 

1139 

1140 ''' 

1141 

1142 b, a, n = decimate_coeffs(q, n, ftype) 

1143 

1144 if zi is None or zi is True: 

1145 zi_ = num.zeros(max(len(a), len(b))-1, dtype=float) 

1146 else: 

1147 zi_ = zi 

1148 

1149 y, zf = signal.lfilter(b, a, x, zi=zi_) 

1150 

1151 if zi is not None: 

1152 return y[n//2+ioff::q].copy(), zf 

1153 else: 

1154 return y[n//2+ioff::q].copy() 

1155 

1156 

1157class UnavailableDecimation(Exception): 

1158 ''' 

1159 Exception raised by :py:func:`decitab` for unavailable decimation factors. 

1160 ''' 

1161 

1162 pass 

1163 

1164 

1165def gcd(a, b, epsilon=1e-7): 

1166 ''' 

1167 Greatest common divisor. 

1168 ''' 

1169 

1170 while b > epsilon*a: 

1171 a, b = b, a % b 

1172 

1173 return a 

1174 

1175 

1176def lcm(a, b): 

1177 ''' 

1178 Least common multiple. 

1179 ''' 

1180 

1181 return a*b // gcd(a, b) 

1182 

1183 

1184def mk_decitab(nmax=100): 

1185 ''' 

1186 Make table with decimation sequences. 

1187 

1188 Decimation from one sampling rate to a lower one is achieved by a 

1189 successive application of :py:func:`decimate` with small integer 

1190 downsampling factors (because using large downsampling factors can make the 

1191 decimation unstable or slow). This function sets up a table with downsample 

1192 sequences for factors up to ``nmax``. 

1193 ''' 

1194 

1195 tab = GlobalVars.decitab 

1196 for i in range(1, 10): 

1197 for j in range(1, i+1): 

1198 for k in range(1, j+1): 

1199 for l_ in range(1, k+1): 

1200 for m in range(1, l_+1): 

1201 p = i*j*k*l_*m 

1202 if p > nmax: 

1203 break 

1204 if p not in tab: 

1205 tab[p] = (i, j, k, l_, m) 

1206 if i*j*k*l_ > nmax: 

1207 break 

1208 if i*j*k > nmax: 

1209 break 

1210 if i*j > nmax: 

1211 break 

1212 if i > nmax: 

1213 break 

1214 

1215 GlobalVars.decitab_nmax = nmax 

1216 

1217 

1218def zfmt(n): 

1219 return '%%0%ii' % (int(math.log10(n - 1)) + 1) 

1220 

1221 

1222def _year_to_time(year): 

1223 tt = (year, 1, 1, 0, 0, 0) 

1224 return to_time_float(calendar.timegm(tt)) 

1225 

1226 

1227def _working_year(year): 

1228 try: 

1229 tt = (year, 1, 1, 0, 0, 0) 

1230 t = calendar.timegm(tt) 

1231 tt2_ = time.gmtime(t) 

1232 tt2 = tuple(tt2_)[:6] 

1233 if tt != tt2: 

1234 return False 

1235 

1236 s = '%i-01-01 00:00:00' % year 

1237 s2 = time.strftime('%Y-%m-%d %H:%M:%S', tt2_) 

1238 if s != s2: 

1239 return False 

1240 

1241 t3 = str_to_time(s2, format='%Y-%m-%d %H:%M:%S') 

1242 s3 = time_to_str(t3, format='%Y-%m-%d %H:%M:%S') 

1243 if s3 != s2: 

1244 return False 

1245 

1246 except Exception: 

1247 return False 

1248 

1249 return True 

1250 

1251 

1252def working_system_time_range(year_min_lim=None, year_max_lim=None): 

1253 ''' 

1254 Check time range supported by the systems's time conversion functions. 

1255 

1256 Returns system time stamps of start of year of first/last fully supported 

1257 year span. If this is before 1900 or after 2100, return first/last century 

1258 which is fully supported. 

1259 

1260 :returns: ``(tmin, tmax, year_min, year_max)`` 

1261 ''' 

1262 

1263 year0 = 2000 

1264 year_min = year0 

1265 year_max = year0 

1266 

1267 itests = list(range(101)) 

1268 for i in range(19): 

1269 itests.append(200 + i*100) 

1270 

1271 for i in itests: 

1272 year = year0 - i 

1273 if year_min_lim is not None and year < year_min_lim: 

1274 year_min = year_min_lim 

1275 break 

1276 elif not _working_year(year): 

1277 break 

1278 else: 

1279 year_min = year 

1280 

1281 for i in itests: 

1282 year = year0 + i 

1283 if year_max_lim is not None and year > year_max_lim: 

1284 year_max = year_max_lim 

1285 break 

1286 elif not _working_year(year + 1): 

1287 break 

1288 else: 

1289 year_max = year 

1290 

1291 return ( 

1292 _year_to_time(year_min), 

1293 _year_to_time(year_max), 

1294 year_min, year_max) 

1295 

1296 

1297g_working_system_time_range = None 

1298 

1299 

1300def get_working_system_time_range(): 

1301 ''' 

1302 Caching variant of :py:func:`working_system_time_range`. 

1303 ''' 

1304 

1305 global g_working_system_time_range 

1306 if g_working_system_time_range is None: 

1307 g_working_system_time_range = working_system_time_range() 

1308 

1309 return g_working_system_time_range 

1310 

1311 

1312def is_working_time(t): 

1313 tmin, tmax, _, _ = get_working_system_time_range() 

1314 return tmin <= t <= tmax 

1315 

1316 

1317def julian_day_of_year(timestamp): 

1318 ''' 

1319 Get the day number after the 1st of January of year in ``timestamp``. 

1320 

1321 :returns: day number as int 

1322 ''' 

1323 

1324 return time.gmtime(int(timestamp)).tm_yday 

1325 

1326 

1327def hour_start(timestamp): 

1328 ''' 

1329 Get beginning of hour for any point in time. 

1330 

1331 :param timestamp: time instant as system timestamp (in seconds) 

1332 

1333 :returns: instant of hour start as system timestamp 

1334 ''' 

1335 

1336 tt = time.gmtime(int(timestamp)) 

1337 tts = tt[0:4] + (0, 0) 

1338 return to_time_float(calendar.timegm(tts)) 

1339 

1340 

1341def day_start(timestamp): 

1342 ''' 

1343 Get beginning of day for any point in time. 

1344 

1345 :param timestamp: time instant as system timestamp (in seconds) 

1346 

1347 :returns: instant of day start as system timestamp 

1348 ''' 

1349 

1350 tt = time.gmtime(int(timestamp)) 

1351 tts = tt[0:3] + (0, 0, 0) 

1352 return to_time_float(calendar.timegm(tts)) 

1353 

1354 

1355def month_start(timestamp): 

1356 ''' 

1357 Get beginning of month for any point in time. 

1358 

1359 :param timestamp: time instant as system timestamp (in seconds) 

1360 

1361 :returns: instant of month start as system timestamp 

1362 ''' 

1363 

1364 tt = time.gmtime(int(timestamp)) 

1365 tts = tt[0:2] + (1, 0, 0, 0) 

1366 return to_time_float(calendar.timegm(tts)) 

1367 

1368 

1369def year_start(timestamp): 

1370 ''' 

1371 Get beginning of year for any point in time. 

1372 

1373 :param timestamp: time instant as system timestamp (in seconds) 

1374 

1375 :returns: instant of year start as system timestamp 

1376 ''' 

1377 

1378 tt = time.gmtime(int(timestamp)) 

1379 tts = tt[0:1] + (1, 1, 0, 0, 0) 

1380 return to_time_float(calendar.timegm(tts)) 

1381 

1382 

1383def iter_days(tmin, tmax): 

1384 ''' 

1385 Yields begin and end of days until given time span is covered. 

1386 

1387 :param tmin,tmax: input time span 

1388 

1389 :yields: tuples with (begin, end) of days as system timestamps 

1390 ''' 

1391 

1392 t = day_start(tmin) 

1393 while t < tmax: 

1394 tend = day_start(t + 26*60*60) 

1395 yield t, tend 

1396 t = tend 

1397 

1398 

1399def iter_months(tmin, tmax): 

1400 ''' 

1401 Yields begin and end of months until given time span is covered. 

1402 

1403 :param tmin,tmax: input time span 

1404 

1405 :yields: tuples with (begin, end) of months as system timestamps 

1406 ''' 

1407 

1408 t = month_start(tmin) 

1409 while t < tmax: 

1410 tend = month_start(t + 24*60*60*33) 

1411 yield t, tend 

1412 t = tend 

1413 

1414 

1415def iter_years(tmin, tmax): 

1416 ''' 

1417 Yields begin and end of years until given time span is covered. 

1418 

1419 :param tmin,tmax: input time span 

1420 

1421 :yields: tuples with (begin, end) of years as system timestamps 

1422 ''' 

1423 

1424 t = year_start(tmin) 

1425 while t < tmax: 

1426 tend = year_start(t + 24*60*60*369) 

1427 yield t, tend 

1428 t = tend 

1429 

1430 

1431def today(): 

1432 return day_start(time.time()) 

1433 

1434 

1435def tomorrow(): 

1436 return day_start(time.time() + 24*60*60) 

1437 

1438 

1439def decitab(n): 

1440 ''' 

1441 Get integer decimation sequence for given downampling factor. 

1442 

1443 :param n: target decimation factor 

1444 

1445 :returns: tuple with downsampling sequence 

1446 ''' 

1447 

1448 if n > GlobalVars.decitab_nmax: 

1449 mk_decitab(n*2) 

1450 if n not in GlobalVars.decitab: 

1451 raise UnavailableDecimation('ratio = %g' % n) 

1452 return GlobalVars.decitab[n] 

1453 

1454 

1455def ctimegm(s, format='%Y-%m-%d %H:%M:%S'): 

1456 ''' 

1457 Convert string representing UTC time to system time. 

1458 

1459 :param s: string to be interpreted 

1460 :param format: format string passed to :py:func:`time.strptime` 

1461 

1462 :returns: system time stamp 

1463 

1464 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime. 

1465 

1466 .. note:: 

1467 This function is to be replaced by :py:func:`str_to_time`. 

1468 ''' 

1469 

1470 return calendar.timegm(time.strptime(s, format)) 

1471 

1472 

1473def gmctime(t, format='%Y-%m-%d %H:%M:%S'): 

1474 ''' 

1475 Get string representation from system time, UTC. 

1476 

1477 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime. 

1478 

1479 .. note:: 

1480 This function is to be repaced by :py:func:`time_to_str`. 

1481 ''' 

1482 

1483 return time.strftime(format, time.gmtime(t)) 

1484 

1485 

1486def gmctime_v(t, format='%a, %d %b %Y %H:%M:%S'): 

1487 ''' 

1488 Get string representation from system time, UTC. Same as 

1489 :py:func:`gmctime` but with a more verbose default format. 

1490 

1491 .. note:: 

1492 This function is to be replaced by :py:func:`time_to_str`. 

1493 ''' 

1494 

1495 return time.strftime(format, time.gmtime(t)) 

1496 

1497 

1498def gmctime_fn(t, format='%Y-%m-%d_%H-%M-%S'): 

1499 ''' 

1500 Get string representation from system time, UTC. Same as 

1501 :py:func:`gmctime` but with a default usable in filenames. 

1502 

1503 .. note:: 

1504 This function is to be replaced by :py:func:`time_to_str`. 

1505 ''' 

1506 

1507 return time.strftime(format, time.gmtime(t)) 

1508 

1509 

1510class TimeStrError(Exception): 

1511 ''' 

1512 Raised for invalid time strings. 

1513 ''' 

1514 pass 

1515 

1516 

1517class FractionalSecondsMissing(TimeStrError): 

1518 ''' 

1519 Exception raised by :py:func:`str_to_time` when the given string lacks 

1520 fractional seconds. 

1521 ''' 

1522 

1523 pass 

1524 

1525 

1526class FractionalSecondsWrongNumberOfDigits(TimeStrError): 

1527 ''' 

1528 Exception raised by :py:func:`str_to_time` when the given string has an 

1529 incorrect number of digits in the fractional seconds part. 

1530 ''' 

1531 

1532 pass 

1533 

1534 

1535def _endswith_n(s, endings): 

1536 for ix, x in enumerate(endings): 

1537 if s.endswith(x): 

1538 return ix 

1539 return -1 

1540 

1541 

1542def str_to_time(s, format='%Y-%m-%d %H:%M:%S.OPTFRAC'): 

1543 ''' 

1544 Convert string representing UTC time to floating point system time. 

1545 

1546 :param s: string representing UTC time 

1547 :param format: time string format 

1548 :returns: system time stamp as floating point value 

1549 

1550 Uses the semantics of :py:func:`time.strptime` but allows for fractional 

1551 seconds. If the format ends with ``'.FRAC'``, anything after a dot is 

1552 interpreted as fractional seconds. If the format ends with ``'.OPTFRAC'``, 

1553 the fractional part, including the dot is made optional. The latter has the 

1554 consequence, that the time strings and the format may not contain any other 

1555 dots. If the format ends with ``'.xFRAC'`` where x is 1, 2, or 3, it is 

1556 ensured, that exactly that number of digits are present in the fractional 

1557 seconds. 

1558 ''' 

1559 

1560 time_float = get_time_float() 

1561 

1562 if util_ext is not None: 

1563 try: 

1564 t, tfrac = util_ext.stt(s, format) 

1565 except util_ext.UtilExtError as e: 

1566 raise TimeStrError( 

1567 '%s, string=%s, format=%s' % (str(e), s, format)) 

1568 

1569 return time_float(t)+tfrac 

1570 

1571 fracsec = 0. 

1572 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC' 

1573 

1574 iend = _endswith_n(format, fixed_endings) 

1575 if iend != -1: 

1576 dotpos = s.rfind('.') 

1577 if dotpos == -1: 

1578 raise FractionalSecondsMissing( 

1579 'string=%s, format=%s' % (s, format)) 

1580 

1581 if iend > 0 and iend != (len(s)-dotpos-1): 

1582 raise FractionalSecondsWrongNumberOfDigits( 

1583 'string=%s, format=%s' % (s, format)) 

1584 

1585 format = format[:-len(fixed_endings[iend])] 

1586 fracsec = float(s[dotpos:]) 

1587 s = s[:dotpos] 

1588 

1589 elif format.endswith('.OPTFRAC'): 

1590 dotpos = s.rfind('.') 

1591 format = format[:-8] 

1592 if dotpos != -1 and len(s[dotpos:]) > 1: 

1593 fracsec = float(s[dotpos:]) 

1594 

1595 if dotpos != -1: 

1596 s = s[:dotpos] 

1597 

1598 try: 

1599 return time_float(calendar.timegm(time.strptime(s, format))) \ 

1600 + fracsec 

1601 except ValueError as e: 

1602 raise TimeStrError('%s, string=%s, format=%s' % (str(e), s, format)) 

1603 

1604 

1605stt = str_to_time 

1606 

1607 

1608def str_to_time_fillup(s): 

1609 ''' 

1610 Default :py:func:`str_to_time` with filling in of missing values. 

1611 

1612 Allows e.g. `'2010-01-01 00:00:00'` as `'2010-01-01 00:00'`, 

1613 `'2010-01-01 00'`, ..., or `'2010'`. 

1614 ''' 

1615 

1616 if s == 'now': 

1617 return time.time() 

1618 

1619 if len(s) in (4, 7, 10, 13, 16): 

1620 s += '0000-01-01 00:00:00'[len(s):] 

1621 

1622 return str_to_time(s) 

1623 

1624 

1625def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'): 

1626 ''' 

1627 Get string representation for floating point system time. 

1628 

1629 :param t: floating point system time 

1630 :param format: time string format 

1631 :returns: string representing UTC time 

1632 

1633 Uses the semantics of :py:func:`time.strftime` but additionally allows for 

1634 fractional seconds. If ``format`` contains ``'.xFRAC'``, where ``x`` is a 

1635 digit between 1 and 9, this is replaced with the fractional part of ``t`` 

1636 with ``x`` digits precision. 

1637 ''' 

1638 

1639 if pyrocko.grumpy > 0: 

1640 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise') 

1641 

1642 if isinstance(format, int): 

1643 if format > 0: 

1644 format = '%Y-%m-%d %H:%M:%S.' + '%iFRAC' % format 

1645 else: 

1646 format = '%Y-%m-%d %H:%M:%S' 

1647 

1648 if util_ext is not None: 

1649 t0 = num.floor(t) 

1650 try: 

1651 return util_ext.tts(int(t0), float(t - t0), format) 

1652 except util_ext.UtilExtError as e: 

1653 raise TimeStrError( 

1654 '%s, timestamp=%f, format=%s' % (str(e), t, format)) 

1655 

1656 if not GlobalVars.re_frac: 

1657 GlobalVars.re_frac = re.compile(r'\.[1-9]FRAC') 

1658 GlobalVars.frac_formats = dict( 

1659 [('.%sFRAC' % x, '%.'+x+'f') for x in '123456789']) 

1660 

1661 ts = float(num.floor(t)) 

1662 tfrac = t-ts 

1663 

1664 m = GlobalVars.re_frac.search(format) 

1665 if m: 

1666 sfrac = (GlobalVars.frac_formats[m.group(0)] % tfrac) 

1667 if sfrac[0] == '1': 

1668 ts += 1. 

1669 

1670 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1) 

1671 

1672 return time.strftime(format, time.gmtime(ts)) 

1673 

1674 

1675tts = time_to_str 

1676_abbr_weekday = 'Mon Tue Wed Thu Fri Sat Sun'.split() 

1677_abbr_month = 'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split() 

1678 

1679 

1680def mystrftime(fmt=None, tt=None, milliseconds=0): 

1681 # Needed by Snuffler for the time axis. In other cases `time_to_str` 

1682 # should be used. 

1683 

1684 if fmt is None: 

1685 fmt = '%Y-%m-%d %H:%M:%S .%r' 

1686 

1687 # Get these two locale independent, needed by Snuffler. 

1688 # Setting the locale seems to have no effect. 

1689 fmt = fmt.replace('%a', _abbr_weekday[tt.tm_wday]) 

1690 fmt = fmt.replace('%b', _abbr_month[tt.tm_mon-1]) 

1691 

1692 fmt = fmt.replace('%r', '%03i' % int(round(milliseconds))) 

1693 fmt = fmt.replace('%u', '%06i' % int(round(milliseconds*1000))) 

1694 fmt = fmt.replace('%n', '%09i' % int(round(milliseconds*1000000))) 

1695 return time.strftime(fmt, tt) 

1696 

1697 

1698def gmtime_x(timestamp): 

1699 etimestamp = float(num.floor(timestamp)) 

1700 tt = time.gmtime(etimestamp) 

1701 ms = (timestamp-etimestamp)*1000 

1702 return tt, ms 

1703 

1704 

1705def plural_s(n): 

1706 if not isinstance(n, int): 

1707 n = len(n) 

1708 

1709 return 's' if n != 1 else '' 

1710 

1711 

1712def ensuredirs(dst): 

1713 ''' 

1714 Create all intermediate path components for a target path. 

1715 

1716 :param dst: target path 

1717 

1718 The leaf part of the target path is not created (use :py:func:`ensuredir` 

1719 if a the target path is a directory to be created). 

1720 ''' 

1721 

1722 d, x = os.path.split(dst.rstrip(os.sep)) 

1723 dirs = [] 

1724 while d and not os.path.exists(d): 

1725 dirs.append(d) 

1726 d, x = os.path.split(d) 

1727 

1728 dirs.reverse() 

1729 

1730 for d in dirs: 

1731 try: 

1732 os.mkdir(d) 

1733 except OSError as e: 

1734 if not e.errno == errno.EEXIST: 

1735 raise 

1736 

1737 

1738def ensuredir(dst): 

1739 ''' 

1740 Create directory and all intermediate path components to it as needed. 

1741 

1742 :param dst: directory name 

1743 

1744 Nothing is done if the given target already exists. 

1745 ''' 

1746 

1747 if os.path.exists(dst): 

1748 return 

1749 

1750 dst.rstrip(os.sep) 

1751 

1752 ensuredirs(dst) 

1753 try: 

1754 os.mkdir(dst) 

1755 except OSError as e: 

1756 if not e.errno == errno.EEXIST: 

1757 raise 

1758 

1759 

1760def reuse(x): 

1761 ''' 

1762 Get unique instance of an object. 

1763 

1764 :param x: hashable object 

1765 :returns: reference to x or an equivalent object 

1766 

1767 Cache object ``x`` in a global dict for reuse, or if x already 

1768 is in that dict, return a reference to it. 

1769 ''' 

1770 

1771 grs = GlobalVars.reuse_store 

1772 if x not in grs: 

1773 grs[x] = x 

1774 return grs[x] 

1775 

1776 

1777def deuse(x): 

1778 grs = GlobalVars.reuse_store 

1779 if x in grs: 

1780 del grs[x] 

1781 

1782 

1783class Anon(object): 

1784 ''' 

1785 Dict-to-object utility. 

1786 

1787 Any given arguments are stored as attributes. 

1788 

1789 Example:: 

1790 

1791 a = Anon(x=1, y=2) 

1792 print a.x, a.y 

1793 ''' 

1794 

1795 def __init__(self, **dict): 

1796 for k in dict: 

1797 self.__dict__[k] = dict[k] 

1798 

1799 

1800def iter_select_files( 

1801 paths, 

1802 include=None, 

1803 exclude=None, 

1804 selector=None, 

1805 show_progress=True, 

1806 pass_through=None): 

1807 

1808 ''' 

1809 Recursively select files (generator variant). 

1810 

1811 See :py:func:`select_files`. 

1812 ''' 

1813 

1814 if show_progress: 

1815 progress_beg('selecting files...') 

1816 

1817 ngood = 0 

1818 check_include = None 

1819 if include is not None: 

1820 rinclude = re.compile(include) 

1821 

1822 def check_include(path): 

1823 m = rinclude.search(path) 

1824 if not m: 

1825 return False 

1826 

1827 if selector is None: 

1828 return True 

1829 

1830 infos = Anon(**m.groupdict()) 

1831 return selector(infos) 

1832 

1833 check_exclude = None 

1834 if exclude is not None: 

1835 rexclude = re.compile(exclude) 

1836 

1837 def check_exclude(path): 

1838 return not bool(rexclude.search(path)) 

1839 

1840 if check_include and check_exclude: 

1841 

1842 def check(path): 

1843 return check_include(path) and check_exclude(path) 

1844 

1845 elif check_include: 

1846 check = check_include 

1847 

1848 elif check_exclude: 

1849 check = check_exclude 

1850 

1851 else: 

1852 check = None 

1853 

1854 if isinstance(paths, str): 

1855 paths = [paths] 

1856 

1857 for path in paths: 

1858 if pass_through and pass_through(path): 

1859 if check is None or check(path): 

1860 yield path 

1861 

1862 elif os.path.isdir(path): 

1863 for (dirpath, dirnames, filenames) in os.walk(path): 

1864 dirnames.sort() 

1865 filenames.sort() 

1866 for filename in filenames: 

1867 path = op.join(dirpath, filename) 

1868 if check is None or check(path): 

1869 yield os.path.abspath(path) 

1870 ngood += 1 

1871 else: 

1872 if check is None or check(path): 

1873 yield os.path.abspath(path) 

1874 ngood += 1 

1875 

1876 if show_progress: 

1877 progress_end('%i file%s selected.' % (ngood, plural_s(ngood))) 

1878 

1879 

1880def select_files( 

1881 paths, include=None, exclude=None, selector=None, show_progress=True, 

1882 regex=None): 

1883 

1884 ''' 

1885 Recursively select files. 

1886 

1887 :param paths: entry path names 

1888 :param include: pattern for conditional inclusion 

1889 :param exclude: pattern for conditional exclusion 

1890 :param selector: callback for conditional inclusion 

1891 :param show_progress: if True, indicate start and stop of processing 

1892 :param regex: alias for ``include`` (backwards compatibility) 

1893 :returns: list of path names 

1894 

1895 Recursively finds all files under given entry points ``paths``. If 

1896 parameter ``include`` is a regular expression, only files with matching 

1897 path names are included. If additionally parameter ``selector`` is given a 

1898 callback function, only files for which the callback returns ``True`` are 

1899 included. The callback should take a single argument. The callback is 

1900 called with a single argument, an object, having as attributes, any named 

1901 groups given in ``include``. 

1902 

1903 Examples 

1904 

1905 To find all files ending in ``'.mseed'`` or ``'.msd'``:: 

1906 

1907 select_files(paths, 

1908 include=r'\\.(mseed|msd)$') 

1909 

1910 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for 

1911 the year:: 

1912 

1913 select_files(paths, 

1914 include=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$', 

1915 selector=(lambda x: int(x.year) == 2009)) 

1916 ''' 

1917 if regex is not None: 

1918 assert include is None 

1919 include = regex 

1920 

1921 return list(iter_select_files( 

1922 paths, include, exclude, selector, show_progress)) 

1923 

1924 

1925def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'): 

1926 ''' 

1927 Convert positive integer to a base36 string. 

1928 ''' 

1929 

1930 if not isinstance(number, (int, long)): 

1931 raise TypeError('number must be an integer') 

1932 if number < 0: 

1933 raise ValueError('number must be positive') 

1934 

1935 # Special case for small numbers 

1936 if number < 36: 

1937 return alphabet[number] 

1938 

1939 base36 = '' 

1940 while number != 0: 

1941 number, i = divmod(number, 36) 

1942 base36 = alphabet[i] + base36 

1943 

1944 return base36 

1945 

1946 

1947def base36decode(number): 

1948 ''' 

1949 Decode base36 endcoded positive integer. 

1950 ''' 

1951 

1952 return int(number, 36) 

1953 

1954 

1955class UnpackError(Exception): 

1956 ''' 

1957 Exception raised when :py:func:`unpack_fixed` encounters an error. 

1958 ''' 

1959 

1960 pass 

1961 

1962 

1963ruler = ''.join(['%-10i' % i for i in range(8)]) \ 

1964 + '\n' + '0123456789' * 8 + '\n' 

1965 

1966 

1967def unpack_fixed(format, line, *callargs): 

1968 ''' 

1969 Unpack fixed format string, as produced by many fortran codes. 

1970 

1971 :param format: format specification 

1972 :param line: string to be processed 

1973 :param callargs: callbacks for callback fields in the format 

1974 

1975 The format is described by a string of comma-separated fields. Each field 

1976 is defined by a character for the field type followed by the field width. A 

1977 questionmark may be appended to the field description to allow the argument 

1978 to be optional (The data string is then allowed to be filled with blanks 

1979 and ``None`` is returned in this case). 

1980 

1981 The following field types are available: 

1982 

1983 ==== ================================================================ 

1984 Type Description 

1985 ==== ================================================================ 

1986 A string (full field width is extracted) 

1987 a string (whitespace at the beginning and the end is removed) 

1988 i integer value 

1989 f floating point value 

1990 @ special type, a callback must be given for the conversion 

1991 x special field type to skip parts of the string 

1992 ==== ================================================================ 

1993 ''' 

1994 

1995 ipos = 0 

1996 values = [] 

1997 icall = 0 

1998 for form in format.split(','): 

1999 form = form.strip() 

2000 optional = form[-1] == '?' 

2001 form = form.rstrip('?') 

2002 typ = form[0] 

2003 ln = int(form[1:]) 

2004 s = line[ipos:ipos+ln] 

2005 cast = { 

2006 'x': None, 

2007 'A': str, 

2008 'a': lambda x: x.strip(), 

2009 'i': int, 

2010 'f': float, 

2011 '@': 'extra'}[typ] 

2012 

2013 if cast == 'extra': 

2014 cast = callargs[icall] 

2015 icall += 1 

2016 

2017 if cast is not None: 

2018 if optional and s.strip() == '': 

2019 values.append(None) 

2020 else: 

2021 try: 

2022 values.append(cast(s)) 

2023 except Exception: 

2024 mark = [' '] * 80 

2025 mark[ipos:ipos+ln] = ['^'] * ln 

2026 mark = ''.join(mark) 

2027 raise UnpackError( 

2028 'Invalid cast to type "%s" at position [%i:%i] of ' 

2029 'line: \n%s%s\n%s' % ( 

2030 typ, ipos, ipos+ln, ruler, line.rstrip(), mark)) 

2031 

2032 ipos += ln 

2033 

2034 return values 

2035 

2036 

2037_pattern_cache = {} 

2038 

2039 

2040def _nslc_pattern(pattern): 

2041 if pattern not in _pattern_cache: 

2042 rpattern = re.compile(fnmatch.translate(pattern), re.I) 

2043 _pattern_cache[pattern] = rpattern 

2044 else: 

2045 rpattern = _pattern_cache[pattern] 

2046 

2047 return rpattern 

2048 

2049 

2050def match_nslc(patterns, nslc): 

2051 ''' 

2052 Match network-station-location-channel code against pattern or list of 

2053 patterns. 

2054 

2055 :param patterns: pattern or list of patterns 

2056 :param nslc: tuple with (network, station, location, channel) as strings 

2057 

2058 :returns: ``True`` if the pattern matches or if any of the given patterns 

2059 match; or ``False``. 

2060 

2061 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq]. 

2062 

2063 Example:: 

2064 

2065 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True 

2066 ''' 

2067 

2068 if isinstance(patterns, str): 

2069 patterns = [patterns] 

2070 

2071 if not isinstance(nslc, str): 

2072 s = '.'.join(nslc) 

2073 else: 

2074 s = nslc 

2075 

2076 for pattern in patterns: 

2077 if _nslc_pattern(pattern).match(s): 

2078 return True 

2079 

2080 return False 

2081 

2082 

2083def match_nslcs(patterns, nslcs): 

2084 ''' 

2085 Get network-station-location-channel codes that match given pattern or any 

2086 of several given patterns. 

2087 

2088 :param patterns: pattern or list of patterns 

2089 :param nslcs: list of (network, station, location, channel) tuples 

2090 

2091 See also :py:func:`match_nslc` 

2092 ''' 

2093 

2094 matching = [] 

2095 for nslc in nslcs: 

2096 if match_nslc(patterns, nslc): 

2097 matching.append(nslc) 

2098 

2099 return matching 

2100 

2101 

2102class Timeout(Exception): 

2103 pass 

2104 

2105 

2106def create_lockfile(fn, timeout=None, timewarn=10.): 

2107 t0 = time.time() 

2108 

2109 while True: 

2110 try: 

2111 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL) 

2112 os.close(f) 

2113 return 

2114 

2115 except OSError as e: 

2116 if e.errno in (errno.EEXIST, 13): # 13 occurs on windows 

2117 pass # retry 

2118 else: 

2119 raise 

2120 

2121 tnow = time.time() 

2122 

2123 if timeout is not None and tnow - t0 > timeout: 

2124 raise Timeout( 

2125 'Timeout (%gs) occured while waiting to get exclusive ' 

2126 'access to: %s' % (timeout, fn)) 

2127 

2128 if timewarn is not None and tnow - t0 > timewarn: 

2129 logger.warning( 

2130 'Waiting since %gs to get exclusive access to: %s' % ( 

2131 timewarn, fn)) 

2132 

2133 timewarn *= 2 

2134 

2135 time.sleep(0.01) 

2136 

2137 

2138def delete_lockfile(fn): 

2139 os.unlink(fn) 

2140 

2141 

2142class Lockfile(Exception): 

2143 

2144 def __init__(self, path, timeout=5, timewarn=10.): 

2145 self._path = path 

2146 self._timeout = timeout 

2147 self._timewarn = timewarn 

2148 

2149 def __enter__(self): 

2150 create_lockfile( 

2151 self._path, timeout=self._timeout, timewarn=self._timewarn) 

2152 return None 

2153 

2154 def __exit__(self, type, value, traceback): 

2155 delete_lockfile(self._path) 

2156 

2157 

2158class SoleError(Exception): 

2159 ''' 

2160 Exception raised by objects of type :py:class:`Sole`, when an concurrent 

2161 instance is running. 

2162 ''' 

2163 

2164 pass 

2165 

2166 

2167class Sole(object): 

2168 ''' 

2169 Use POSIX advisory file locking to ensure that only a single instance of a 

2170 program is running. 

2171 

2172 :param pid_path: path to lockfile to be used 

2173 

2174 Usage:: 

2175 

2176 from pyrocko.util import Sole, SoleError, setup_logging 

2177 import os 

2178 

2179 setup_logging('my_program') 

2180 

2181 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock') 

2182 try: 

2183 sole = Sole(pid_path) 

2184 

2185 except SoleError, e: 

2186 logger.fatal( str(e) ) 

2187 sys.exit(1) 

2188 ''' 

2189 

2190 def __init__(self, pid_path): 

2191 self._pid_path = pid_path 

2192 self._other_running = False 

2193 ensuredirs(self._pid_path) 

2194 self._lockfile = None 

2195 self._os = os 

2196 

2197 if not fcntl: 

2198 raise SoleError( 

2199 'Python standard library module "fcntl" not available on ' 

2200 'this platform.') 

2201 

2202 self._fcntl = fcntl 

2203 

2204 try: 

2205 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY) 

2206 except Exception: 

2207 raise SoleError( 

2208 'Cannot open lockfile (path = %s)' % self._pid_path) 

2209 

2210 try: 

2211 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB) 

2212 

2213 except IOError: 

2214 self._other_running = True 

2215 try: 

2216 f = open(self._pid_path, 'r') 

2217 pid = f.read().strip() 

2218 f.close() 

2219 except Exception: 

2220 pid = '?' 

2221 

2222 raise SoleError('Other instance is running (pid = %s)' % pid) 

2223 

2224 try: 

2225 os.ftruncate(self._lockfile, 0) 

2226 os.write(self._lockfile, '%i\n' % os.getpid()) 

2227 os.fsync(self._lockfile) 

2228 

2229 except Exception: 

2230 # the pid is only stored for user information, so this is allowed 

2231 # to fail 

2232 pass 

2233 

2234 def __del__(self): 

2235 if not self._other_running: 

2236 if self._lockfile is not None: 

2237 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN) 

2238 self._os.close(self._lockfile) 

2239 try: 

2240 self._os.unlink(self._pid_path) 

2241 except Exception: 

2242 pass 

2243 

2244 

2245re_escapequotes = re.compile(r"(['\\])") 

2246 

2247 

2248def escapequotes(s): 

2249 return re_escapequotes.sub(r'\\\1', s) 

2250 

2251 

2252class TableWriter(object): 

2253 ''' 

2254 Write table of space separated values to a file. 

2255 

2256 :param f: file like object 

2257 

2258 Strings containing spaces are quoted on output. 

2259 ''' 

2260 

2261 def __init__(self, f): 

2262 self._f = f 

2263 

2264 def writerow(self, row, minfieldwidths=None): 

2265 

2266 ''' 

2267 Write one row of values to underlying file. 

2268 

2269 :param row: iterable of values 

2270 :param minfieldwidths: minimum field widths for the values 

2271 

2272 Each value in in ``row`` is converted to a string and optionally padded 

2273 with blanks. The resulting strings are output separated with blanks. If 

2274 any values given are strings and if they contain whitespace, they are 

2275 quoted with single quotes, and any internal single quotes are 

2276 backslash-escaped. 

2277 ''' 

2278 

2279 out = [] 

2280 

2281 for i, x in enumerate(row): 

2282 w = 0 

2283 if minfieldwidths and i < len(minfieldwidths): 

2284 w = minfieldwidths[i] 

2285 

2286 if isinstance(x, str): 

2287 if re.search(r"\s|'", x): 

2288 x = "'%s'" % escapequotes(x) 

2289 

2290 x = x.ljust(w) 

2291 else: 

2292 x = str(x).rjust(w) 

2293 

2294 out.append(x) 

2295 

2296 self._f.write(' '.join(out).rstrip() + '\n') 

2297 

2298 

2299class TableReader(object): 

2300 

2301 ''' 

2302 Read table of space separated values from a file. 

2303 

2304 :param f: file-like object 

2305 

2306 This uses Pythons shlex module to tokenize lines. Should deal correctly 

2307 with quoted strings. 

2308 ''' 

2309 

2310 def __init__(self, f): 

2311 self._f = f 

2312 self.eof = False 

2313 

2314 def readrow(self): 

2315 ''' 

2316 Read one row from the underlying file, tokenize it with shlex. 

2317 

2318 :returns: tokenized line as a list of strings. 

2319 ''' 

2320 

2321 line = self._f.readline() 

2322 if not line: 

2323 self.eof = True 

2324 return [] 

2325 line.strip() 

2326 if line.startswith('#'): 

2327 return [] 

2328 

2329 return qsplit(line) 

2330 

2331 

2332def gform(number, significant_digits=3): 

2333 ''' 

2334 Pretty print floating point numbers. 

2335 

2336 Align floating point numbers at the decimal dot. 

2337 

2338 :: 

2339 

2340 | -d.dde+xxx| 

2341 | -d.dde+xx | 

2342 |-ddd. | 

2343 | -dd.d | 

2344 | -d.dd | 

2345 | -0.ddd | 

2346 | -0.0ddd | 

2347 | -0.00ddd | 

2348 | -d.dde-xx | 

2349 | -d.dde-xxx| 

2350 | nan| 

2351 

2352 

2353 The formatted string has length ``significant_digits * 2 + 6``. 

2354 ''' 

2355 

2356 no_exp_range = (pow(10., -1), 

2357 pow(10., significant_digits)) 

2358 width = significant_digits+significant_digits-1+1+1+5 

2359 

2360 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.: 

2361 s = ('%#.*g' % (significant_digits, number)).rstrip('0') 

2362 else: 

2363 s = '%.*E' % (significant_digits-1, number) 

2364 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width) 

2365 if s.strip().lower() == 'nan': 

2366 s = 'nan'.rjust(width) 

2367 return s 

2368 

2369 

2370def human_bytesize(value): 

2371 

2372 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split() 

2373 

2374 if value == 1: 

2375 return '1 Byte' 

2376 

2377 for i, ext in enumerate(exts): 

2378 x = float(value) / 1000**i 

2379 if round(x) < 10. and not value < 1000: 

2380 return '%.1f %s' % (x, ext) 

2381 if round(x) < 1000.: 

2382 return '%.0f %s' % (x, ext) 

2383 

2384 return '%i Bytes' % value 

2385 

2386 

2387re_compatibility = re.compile( 

2388 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' + 

2389 r'(dummy|poel|qseis|qssp))\.' 

2390) 

2391 

2392 

2393def pf_is_old(fn): 

2394 oldstyle = False 

2395 with open(fn, 'r') as f: 

2396 for line in f: 

2397 if re_compatibility.search(line): 

2398 oldstyle = True 

2399 

2400 return oldstyle 

2401 

2402 

2403def pf_upgrade(fn): 

2404 need = pf_is_old(fn) 

2405 if need: 

2406 fn_temp = fn + '.temp' 

2407 

2408 with open(fn, 'r') as fin: 

2409 with open(fn_temp, 'w') as fout: 

2410 for line in fin: 

2411 line = re_compatibility.sub('!pf.', line) 

2412 fout.write(line) 

2413 

2414 os.rename(fn_temp, fn) 

2415 

2416 return need 

2417 

2418 

2419def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'): 

2420 ''' 

2421 Extract leap second information from tzdata. 

2422 

2423 Based on example at http://stackoverflow.com/questions/19332902/\ 

2424 extract-historic-leap-seconds-from-tzdata 

2425 

2426 See also 'man 5 tzfile'. 

2427 ''' 

2428 

2429 from struct import unpack, calcsize 

2430 out = [] 

2431 with open(tzfile, 'rb') as f: 

2432 # read header 

2433 fmt = '>4s c 15x 6l' 

2434 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt, 

2435 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt))) 

2436 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file' 

2437 

2438 # skip over some uninteresting data 

2439 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict( 

2440 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt) 

2441 f.read(calcsize(fmt)) 

2442 

2443 # read leap-seconds 

2444 fmt = '>2l' 

2445 for i in range(leapcnt): 

2446 tleap, nleap = unpack(fmt, f.read(calcsize(fmt))) 

2447 out.append((tleap-nleap+1, nleap)) 

2448 

2449 return out 

2450 

2451 

2452class LeapSecondsError(Exception): 

2453 pass 

2454 

2455 

2456class LeapSecondsOutdated(LeapSecondsError): 

2457 pass 

2458 

2459 

2460class InvalidLeapSecondsFile(LeapSecondsOutdated): 

2461 pass 

2462 

2463 

2464def parse_leap_seconds_list(fn): 

2465 data = [] 

2466 texpires = None 

2467 try: 

2468 t0 = int(round(str_to_time('1900-01-01 00:00:00'))) 

2469 except TimeStrError: 

2470 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800 

2471 

2472 tnow = int(round(time.time())) 

2473 

2474 if not op.exists(fn): 

2475 raise LeapSecondsOutdated('no leap seconds file found') 

2476 

2477 try: 

2478 with open(fn, 'rb') as f: 

2479 for line in f: 

2480 if line.strip().startswith(b'<!DOCTYPE'): 

2481 raise InvalidLeapSecondsFile('invalid leap seconds file') 

2482 

2483 if line.startswith(b'#@'): 

2484 texpires = int(line.split()[1]) + t0 

2485 elif line.startswith(b'#') or len(line) < 5: 

2486 pass 

2487 else: 

2488 toks = line.split() 

2489 t = int(toks[0]) + t0 

2490 nleap = int(toks[1]) - 10 

2491 data.append((t, nleap)) 

2492 

2493 except IOError: 

2494 raise LeapSecondsError('cannot read leap seconds file %s' % fn) 

2495 

2496 if texpires is None or tnow > texpires: 

2497 raise LeapSecondsOutdated('leap seconds list is outdated') 

2498 

2499 return data 

2500 

2501 

2502def read_leap_seconds2(): 

2503 from pyrocko import config 

2504 conf = config.config() 

2505 fn = conf.leapseconds_path 

2506 url = conf.leapseconds_url 

2507 # check for outdated default URL 

2508 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list': 

2509 url = 'https://www.ietf.org/timezones/data/leap-seconds.list' 

2510 logger.info( 

2511 'Leap seconds default URL is now: %s\nUsing new default.' % url) 

2512 

2513 if url == 'https://www.ietf.org/timezones/data/leap-seconds.list': 

2514 url = 'https://hpiers.obspm.fr/iers/bul/bulc/ntp/leap-seconds.list' 

2515 logger.info( 

2516 'Leap seconds default URL is now: %s\nUsing new default.' % url) 

2517 

2518 for i in range(3): 

2519 try: 

2520 return parse_leap_seconds_list(fn) 

2521 

2522 except LeapSecondsOutdated: 

2523 try: 

2524 logger.info('updating leap seconds list...') 

2525 download_file(url, fn) 

2526 

2527 except Exception as e: 

2528 raise LeapSecondsError( 

2529 'cannot download leap seconds list from %s to %s (%s)' 

2530 % (url, fn, e)) 

2531 

2532 raise LeapSecondsError('Could not retrieve/read leap seconds file.') 

2533 

2534 

2535def gps_utc_offset(t_utc): 

2536 ''' 

2537 Time offset t_gps - t_utc for a given t_utc. 

2538 ''' 

2539 ls = read_leap_seconds2() 

2540 i = 0 

2541 if t_utc < ls[0][0]: 

2542 return ls[0][1] - 1 - 9 

2543 

2544 while i < len(ls) - 1: 

2545 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]: 

2546 return ls[i][1] - 9 

2547 i += 1 

2548 

2549 return ls[-1][1] - 9 

2550 

2551 

2552def utc_gps_offset(t_gps): 

2553 ''' 

2554 Time offset t_utc - t_gps for a given t_gps. 

2555 ''' 

2556 ls = read_leap_seconds2() 

2557 

2558 if t_gps < ls[0][0] + ls[0][1] - 9: 

2559 return - (ls[0][1] - 1 - 9) 

2560 

2561 i = 0 

2562 while i < len(ls) - 1: 

2563 if ls[i][0] + ls[i][1] - 9 <= t_gps \ 

2564 and t_gps < ls[i+1][0] + ls[i+1][1] - 9: 

2565 return - (ls[i][1] - 9) 

2566 i += 1 

2567 

2568 return - (ls[-1][1] - 9) 

2569 

2570 

2571def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'): 

2572 import itertools 

2573 import glob 

2574 from pyrocko.io.io_common import FileLoadError 

2575 

2576 def iload_filename(filename, **kwargs): 

2577 try: 

2578 with open(filename, 'rb') as f: 

2579 for cr in iload_fh(f, **kwargs): 

2580 yield cr 

2581 

2582 except FileLoadError as e: 

2583 e.set_context('filename', filename) 

2584 raise 

2585 

2586 def iload_dirname(dirname, **kwargs): 

2587 for entry in os.listdir(dirname): 

2588 fpath = op.join(dirname, entry) 

2589 if op.isfile(fpath): 

2590 for cr in iload_filename(fpath, **kwargs): 

2591 yield cr 

2592 

2593 def iload_glob(pattern, **kwargs): 

2594 

2595 for fn in glob.iglob(pattern): 

2596 for cr in iload_filename(fn, **kwargs): 

2597 yield cr 

2598 

2599 def iload(source, **kwargs): 

2600 if isinstance(source, str): 

2601 if op.isdir(source): 

2602 return iload_dirname(source, **kwargs) 

2603 elif op.isfile(source): 

2604 return iload_filename(source, **kwargs) 

2605 else: 

2606 return iload_glob(source, **kwargs) 

2607 

2608 elif hasattr(source, 'read'): 

2609 return iload_fh(source, **kwargs) 

2610 else: 

2611 return itertools.chain.from_iterable( 

2612 iload(subsource, **kwargs) for subsource in source) 

2613 

2614 iload_filename.__doc__ = ''' 

2615 Read %s information from named file. 

2616 ''' % doc_fmt 

2617 

2618 iload_dirname.__doc__ = ''' 

2619 Read %s information from directory of %s files. 

2620 ''' % (doc_fmt, doc_fmt) 

2621 

2622 iload_glob.__doc__ = ''' 

2623 Read %s information from files matching a glob pattern. 

2624 ''' % doc_fmt 

2625 

2626 iload.__doc__ = ''' 

2627 Load %s information from given source(s) 

2628 

2629 The ``source`` can be specified as the name of a %s file, the name of a 

2630 directory containing %s files, a glob pattern of %s files, an open 

2631 filehandle or an iterator yielding any of the forementioned sources. 

2632 

2633 This function behaves as a generator yielding %s objects. 

2634 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects) 

2635 

2636 for f in iload_filename, iload_dirname, iload_glob, iload: 

2637 f.__module__ = iload_fh.__module__ 

2638 

2639 return iload_filename, iload_dirname, iload_glob, iload 

2640 

2641 

2642class Inconsistency(Exception): 

2643 pass 

2644 

2645 

2646def consistency_check(list_of_tuples, message='values differ:'): 

2647 ''' 

2648 Check for inconsistencies. 

2649 

2650 Given a list of tuples, check that all tuple elements except for first one 

2651 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be 

2652 valid because the coordinates at the two channels are the same. 

2653 ''' 

2654 

2655 if len(list_of_tuples) >= 2: 

2656 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]): 

2657 raise Inconsistency('%s\n' % message + '\n'.join( 

2658 ' %s: %s' % (t[0], ', '.join(str(x) for x in t[1:])) 

2659 for t in list_of_tuples)) 

2660 

2661 

2662class defaultzerodict(dict): 

2663 def __missing__(self, k): 

2664 return 0 

2665 

2666 

2667def mostfrequent(x): 

2668 c = defaultzerodict() 

2669 for e in x: 

2670 c[e] += 1 

2671 

2672 return sorted(list(c.keys()), key=lambda k: c[k])[-1] 

2673 

2674 

2675def consistency_merge(list_of_tuples, 

2676 message='values differ:', 

2677 error='raise', 

2678 merge=mostfrequent): 

2679 

2680 assert error in ('raise', 'warn', 'ignore') 

2681 

2682 if len(list_of_tuples) == 0: 

2683 raise Exception('cannot merge empty sequence') 

2684 

2685 try: 

2686 consistency_check(list_of_tuples, message) 

2687 return list_of_tuples[0][1:] 

2688 except Inconsistency as e: 

2689 if error == 'raise': 

2690 raise 

2691 

2692 elif error == 'warn': 

2693 logger.warning(str(e)) 

2694 

2695 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]]) 

2696 

2697 

2698def short_to_list(nmax, it): 

2699 import itertools 

2700 

2701 if isinstance(it, list): 

2702 return it 

2703 

2704 li = [] 

2705 for i in range(nmax+1): 

2706 try: 

2707 li.append(next(it)) 

2708 except StopIteration: 

2709 return li 

2710 

2711 return itertools.chain(li, it) 

2712 

2713 

2714def parse_md(f): 

2715 try: 

2716 with open(op.join( 

2717 op.dirname(op.abspath(f)), 

2718 'README.md'), 'r') as readme: 

2719 mdstr = readme.read() 

2720 except IOError as e: 

2721 return 'Failed to get README.md: %s' % e 

2722 

2723 # Remve the title 

2724 mdstr = re.sub(r'^# .*\n?', '', mdstr) 

2725 # Append sphinx reference to `pyrocko.` modules 

2726 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr) 

2727 # Convert Subsections to toc-less rubrics 

2728 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr) 

2729 return mdstr 

2730 

2731 

2732def mpl_show(plt): 

2733 import matplotlib 

2734 if matplotlib.get_backend().lower() == 'agg': 

2735 logger.warning('Cannot show() when using matplotlib "agg" backend') 

2736 else: 

2737 plt.show() 

2738 

2739 

2740g_re_qsplit = re.compile( 

2741 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)') 

2742g_re_qsplit_sep = {} 

2743 

2744 

2745def get_re_qsplit(sep): 

2746 if sep is None: 

2747 return g_re_qsplit 

2748 else: 

2749 if sep not in g_re_qsplit_sep: 

2750 assert len(sep) == 1 

2751 assert sep not in '\'"' 

2752 esep = re.escape(sep) 

2753 g_re_qsplit_sep[sep] = re.compile( 

2754 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|' 

2755 + r'([^' + esep + r']+|(?<=' + esep + r')(?=' + esep + r')|^(?=' + esep + r')|(?<=' + esep + r')$)') # noqa 

2756 return g_re_qsplit_sep[sep] 

2757 

2758 

2759g_re_trivial = re.compile(r'\A[^\'"\s]+\Z') 

2760g_re_trivial_sep = {} 

2761 

2762 

2763def get_re_trivial(sep): 

2764 if sep is None: 

2765 return g_re_trivial 

2766 else: 

2767 if sep not in g_re_qsplit_sep: 

2768 assert len(sep) == 1 

2769 assert sep not in '\'"' 

2770 esep = re.escape(sep) 

2771 g_re_trivial_sep[sep] = re.compile(r'\A[^\'"' + esep + r']+\Z') 

2772 

2773 return g_re_trivial_sep[sep] 

2774 

2775 

2776g_re_escape_s = re.compile(r'([\\\'])') 

2777g_re_unescape_s = re.compile(r'\\([\\\'])') 

2778g_re_escape_d = re.compile(r'([\\"])') 

2779g_re_unescape_d = re.compile(r'\\([\\"])') 

2780 

2781 

2782def escape_s(s): 

2783 ''' 

2784 Backslash-escape single-quotes and backslashes. 

2785 

2786 Example: ``Jack's`` => ``Jack\\'s`` 

2787 

2788 ''' 

2789 return g_re_escape_s.sub(r'\\\1', s) 

2790 

2791 

2792def unescape_s(s): 

2793 ''' 

2794 Unescape backslash-escaped single-quotes and backslashes. 

2795 

2796 Example: ``Jack\\'s`` => ``Jack's`` 

2797 ''' 

2798 return g_re_unescape_s.sub(r'\1', s) 

2799 

2800 

2801def escape_d(s): 

2802 ''' 

2803 Backslash-escape double-quotes and backslashes. 

2804 

2805 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"`` 

2806 ''' 

2807 return g_re_escape_d.sub(r'\\\1', s) 

2808 

2809 

2810def unescape_d(s): 

2811 ''' 

2812 Unescape backslash-escaped double-quotes and backslashes. 

2813 

2814 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"`` 

2815 ''' 

2816 return g_re_unescape_d.sub(r'\1', s) 

2817 

2818 

2819def qjoin_s(it, sep=None): 

2820 ''' 

2821 Join sequence of strings into a line, single-quoting non-trivial strings. 

2822 

2823 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"`` 

2824 ''' 

2825 re_trivial = get_re_trivial(sep) 

2826 

2827 if sep is None: 

2828 sep = ' ' 

2829 

2830 return sep.join( 

2831 w if re_trivial.search(w) else "'%s'" % escape_s(w) for w in it) 

2832 

2833 

2834def qjoin_d(it, sep=None): 

2835 ''' 

2836 Join sequence of strings into a line, double-quoting non-trivial strings. 

2837 

2838 Example: ``['55', 'Pete "The Robot" Smith']`` => 

2839 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` 

2840 ''' 

2841 re_trivial = get_re_trivial(sep) 

2842 if sep is None: 

2843 sep = ' ' 

2844 

2845 return sep.join( 

2846 w if re_trivial.search(w) else '"%s"' % escape_d(w) for w in it) 

2847 

2848 

2849def qsplit(s, sep=None): 

2850 ''' 

2851 Split line into list of strings, allowing for quoted strings. 

2852 

2853 Example: ``"55 'Sparrow\\\\'s Island'"`` => 

2854 ``["55", "Sparrow's Island"]``, 

2855 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` => 

2856 ``['55', 'Pete "The Robot" Smith']`` 

2857 ''' 

2858 re_qsplit = get_re_qsplit(sep) 

2859 return [ 

2860 (unescape_d(x[0]) or unescape_s(x[1]) or x[2]) 

2861 for x in re_qsplit.findall(s)] 

2862 

2863 

2864g_have_warned_threadpoolctl = False 

2865 

2866 

2867class threadpool_limits_dummy(object): 

2868 

2869 def __init__(self, *args, **kwargs): 

2870 pass 

2871 

2872 def __enter__(self): 

2873 global g_have_warned_threadpoolctl 

2874 

2875 if not g_have_warned_threadpoolctl: 

2876 logger.warning( 

2877 'Cannot control number of BLAS threads because ' 

2878 '`threadpoolctl` module is not available. You may want to ' 

2879 'install `threadpoolctl`.') 

2880 

2881 g_have_warned_threadpoolctl = True 

2882 

2883 return self 

2884 

2885 def __exit__(self, type, value, traceback): 

2886 pass 

2887 

2888 

2889def get_threadpool_limits(): 

2890 ''' 

2891 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail. 

2892 ''' 

2893 

2894 try: 

2895 from threadpoolctl import threadpool_limits 

2896 return threadpool_limits 

2897 

2898 except ImportError: 

2899 return threadpool_limits_dummy 

2900 

2901 

2902def fmt_summary(entries, widths): 

2903 return ' | '.join( 

2904 entry.ljust(width) for (entry, width) in zip(entries, widths)) 

2905 

2906 

2907def smart_weakref(obj, callback=None): 

2908 if inspect.ismethod(obj): 

2909 return weakref.WeakMethod(obj, callback) 

2910 else: 

2911 return weakref.ref(obj, callback)