1# http://pyrocko.org - GPLv3 

2 

3# 

4# The Pyrocko Developers, 21st Century 

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

6''' 

7Utility functions for Pyrocko. 

8 

9.. _time-handling-mode: 

10 

11High precision time handling mode 

12................................. 

13 

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

15floating point values, or as high precision float (:py:class:`numpy.float128` 

16or :py:class:`numpy.float96`, whichever is available), aliased here as 

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

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

19event 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 config variable 

28:py:gattr:`~pyrocko.config.Config.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` or :py:class:`~pyrocko.trace.Trace` objects), 

35use: 

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 

64from __future__ import division, print_function 

65 

66import time 

67import logging 

68import os 

69import sys 

70import re 

71import calendar 

72import math 

73import fnmatch 

74try: 

75 import fcntl 

76except ImportError: 

77 fcntl = None 

78import optparse 

79import os.path as op 

80import errno 

81 

82import numpy as num 

83from scipy import signal 

84import pyrocko 

85from pyrocko import dummy_progressbar 

86 

87 

88try: 

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

90 from urllib.request import ( 

91 Request, build_opener, HTTPDigestAuthHandler, urlopen as _urlopen) # noqa 

92 from urllib.error import HTTPError, URLError # noqa 

93 

94except ImportError: 

95 from urllib import urlencode, quote, unquote # noqa 

96 from urllib2 import (Request, build_opener, HTTPDigestAuthHandler, # noqa 

97 HTTPError, URLError, urlopen as _urlopen) # noqa 

98 

99try: 

100 import certifi 

101 import ssl 

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

103except ImportError: 

104 g_ssl_context = None 

105 

106 

107class URLErrorSSL(URLError): 

108 

109 def __init__(self, urlerror): 

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

111 

112 def __str__(self): 

113 return ( 

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

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

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

117 'with some SSL setups of today\'s servers: %s' 

118 % URLError.__str__(self)) 

119 

120 

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

122 try: 

123 return urlopen(*args, **kwargs) 

124 except URLError as e: 

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

126 raise URLErrorSSL(e) 

127 else: 

128 raise 

129 

130 

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

132 

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

134 kwargs['context'] = g_ssl_context 

135 

136 return _urlopen(*args, **kwargs) 

137 

138 

139try: 

140 long 

141except NameError: 

142 long = int 

143 

144 

145force_dummy_progressbar = False 

146 

147 

148try: 

149 from pyrocko import util_ext 

150except ImportError: 

151 util_ext = None 

152 

153 

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

155 

156 

157try: 

158 num_full = num.full 

159except AttributeError: 

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

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

162 a.fill(fill_value) 

163 return a 

164 

165try: 

166 num_full_like = num.full_like 

167except AttributeError: 

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

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

170 a.fill(fill_value) 

171 return a 

172 

173 

174g_setup_logging_args = 'pyrocko', 'warning' 

175 

176 

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

178 ''' 

179 Initialize logging. 

180 

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

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

183 'warning', 'error', 'critical') 

184 

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

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

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

188 ''' 

189 

190 global g_setup_logging_args 

191 g_setup_logging_args = (programname, levelname) 

192 

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

194 'info': logging.INFO, 

195 'warning': logging.WARNING, 

196 'error': logging.ERROR, 

197 'critical': logging.CRITICAL} 

198 

199 logging.basicConfig( 

200 level=levels[levelname], 

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

202 

203 

204def subprocess_setup_logging_args(): 

205 ''' 

206 Get arguments from previous call to setup_logging. 

207 

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

209 in the same way as the main process. 

210 ''' 

211 return g_setup_logging_args 

212 

213 

214def data_file(fn): 

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

216 

217 

218class DownloadError(Exception): 

219 pass 

220 

221 

222class PathExists(DownloadError): 

223 pass 

224 

225 

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

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

228 status_callback=None, entries_wanted=None, 

229 recursive=False, header=None): 

230 

231 import requests 

232 from requests.auth import HTTPBasicAuth 

233 from requests.exceptions import HTTPError as req_HTTPError 

234 

235 requests.adapters.DEFAULT_RETRIES = 5 

236 urljoin = requests.compat.urljoin 

237 

238 session = requests.Session() 

239 session.header = header 

240 session.auth = None if username is None\ 

241 else HTTPBasicAuth(username, password) 

242 

243 status = { 

244 'ntotal_files': 0, 

245 'nread_files': 0, 

246 'ntotal_bytes_all_files': 0, 

247 'nread_bytes_all_files': 0, 

248 'ntotal_bytes_current_file': 0, 

249 'nread_bytes_current_file': 0, 

250 'finished': False 

251 } 

252 

253 try: 

254 url_to_size = {} 

255 

256 if callable(status_callback): 

257 status_callback(status) 

258 

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

260 raise DownloadError( 

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

262 ' but recurvise download is False' % url) 

263 

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

265 url += '/' 

266 

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

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

269 r.raise_for_status() 

270 

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

272 

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

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

275 

276 if entries_wanted is not None: 

277 files = [fn for fn in files 

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

279 

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

281 

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

283 if dn.endswith('/') 

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

285 

286 for dn in dirs: 

287 files.extend(parse_directory_tree( 

288 url, subdir=dn)) 

289 

290 return files 

291 

292 def get_content_length(url): 

293 if url not in url_to_size: 

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

295 

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

297 if content_length is None: 

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

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

300 

301 content_length = None 

302 

303 else: 

304 content_length = int(content_length) 

305 status['ntotal_bytes_all_files'] += content_length 

306 if callable(status_callback): 

307 status_callback(status) 

308 

309 url_to_size[url] = content_length 

310 

311 return url_to_size[url] 

312 

313 def download_file(url, fn): 

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

315 ensuredirs(fn) 

316 

317 fsize = get_content_length(url) 

318 status['ntotal_bytes_current_file'] = fsize 

319 status['nread_bytes_current_file'] = 0 

320 if callable(status_callback): 

321 status_callback(status) 

322 

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

324 r.raise_for_status() 

325 

326 frx = 0 

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

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

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

330 f.write(d) 

331 frx += len(d) 

332 

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

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

335 if callable(status_callback): 

336 status_callback(status) 

337 

338 os.rename(fn_tmp, fn) 

339 

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

341 logger.warning( 

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

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

344 

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

346 

347 status['nread_files'] += 1 

348 if callable(status_callback): 

349 status_callback(status) 

350 

351 if recursive: 

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

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

354 

355 files = parse_directory_tree(url) 

356 

357 dsize = 0 

358 for fn in files: 

359 file_url = urljoin(url, fn) 

360 dsize += get_content_length(file_url) or 0 

361 

362 if method == 'calcsize': 

363 return dsize 

364 

365 else: 

366 for fn in files: 

367 file_url = urljoin(url, fn) 

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

369 

370 else: 

371 status['ntotal_files'] += 1 

372 if callable(status_callback): 

373 status_callback(status) 

374 

375 fsize = get_content_length(url) 

376 if method == 'calcsize': 

377 return fsize 

378 else: 

379 download_file(url, fpath) 

380 

381 except req_HTTPError as e: 

382 logging.warning("http error: %s" % e) 

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

384 

385 finally: 

386 status['finished'] = True 

387 if callable(status_callback): 

388 status_callback(status) 

389 session.close() 

390 

391 

392def download_file( 

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

394 **kwargs): 

395 return _download( 

396 url, fpath, username, password, 

397 recursive=False, 

398 status_callback=status_callback, 

399 **kwargs) 

400 

401 

402def download_dir( 

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

404 **kwargs): 

405 

406 return _download( 

407 url, fpath, username, password, 

408 recursive=True, 

409 status_callback=status_callback, 

410 **kwargs) 

411 

412 

413class HPFloatUnavailable(Exception): 

414 pass 

415 

416 

417class dummy_hpfloat(object): 

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

419 raise HPFloatUnavailable( 

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

421 'platform.') 

422 

423 

424if hasattr(num, 'float128'): 

425 hpfloat = num.float128 

426 

427elif hasattr(num, 'float96'): 

428 hpfloat = num.float96 

429 

430else: 

431 hpfloat = dummy_hpfloat 

432 

433 

434g_time_float = None 

435g_time_dtype = None 

436 

437 

438class TimeFloatSettingError(Exception): 

439 pass 

440 

441 

442def use_high_precision_time(enabled): 

443 ''' 

444 Globally force a specific time handling mode. 

445 

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

447 

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

449 :type enabled: bool 

450 

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

452 It can only be called once. 

453 

454 Special attention is required when using multiprocessing on a platform 

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

456 must be set also in the subprocess. 

457 ''' 

458 _setup_high_precision_time_mode(enabled_app=enabled) 

459 

460 

461def _setup_high_precision_time_mode(enabled_app=False): 

462 global g_time_float 

463 global g_time_dtype 

464 

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

466 raise TimeFloatSettingError( 

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

468 'fixed by an earlier call.') 

469 

470 from pyrocko import config 

471 

472 conf = config.config() 

473 enabled_config = conf.use_high_precision_time 

474 

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

476 if enabled_env is not None: 

477 try: 

478 enabled_env = int(enabled_env) == 1 

479 except ValueError: 

480 raise TimeFloatSettingError( 

481 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME ' 

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

483 

484 enabled = enabled_config 

485 mode_from = 'config variable `use_high_precision_time`' 

486 notify = enabled 

487 

488 if enabled_env is not None: 

489 if enabled_env != enabled: 

490 notify = True 

491 enabled = enabled_env 

492 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`' 

493 

494 if enabled_app is not None: 

495 if enabled_app != enabled: 

496 notify = True 

497 enabled = enabled_app 

498 mode_from = 'application override' 

499 

500 logger.debug(''' 

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

502 config: %s 

503 env: %s 

504 app: %s 

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

506 enabled_config, enabled_env, enabled_app, enabled)) 

507 

508 if notify: 

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

510 'activated' if enabled else 'deactivated', 

511 mode_from)) 

512 

513 if enabled: 

514 g_time_float = hpfloat 

515 g_time_dtype = hpfloat 

516 else: 

517 g_time_float = float 

518 g_time_dtype = num.float64 

519 

520 

521def get_time_float(): 

522 ''' 

523 Get the effective float class for timestamps. 

524 

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

526 

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

528 current time handling mode 

529 ''' 

530 global g_time_float 

531 

532 if g_time_float is None: 

533 _setup_high_precision_time_mode() 

534 

535 return g_time_float 

536 

537 

538def get_time_dtype(): 

539 ''' 

540 Get effective NumPy float class to handle timestamps. 

541 

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

543 ''' 

544 

545 global g_time_dtype 

546 

547 if g_time_dtype is None: 

548 _setup_high_precision_time_mode() 

549 

550 return g_time_dtype 

551 

552 

553def to_time_float(t): 

554 ''' 

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

556 

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

558 ''' 

559 return get_time_float()(t) 

560 

561 

562class TimestampTypeError(ValueError): 

563 pass 

564 

565 

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

567 ''' 

568 Type-check variable against current time handling mode. 

569 

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

571 ''' 

572 

573 if t == 0.0: 

574 return 

575 

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

577 message = ( 

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

579 'Pyrocko\'s currently selected time handling mode.\n\n' 

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

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

582 t, type(t), get_time_float())) 

583 

584 if error == 'raise': 

585 raise TimestampTypeError(message) 

586 elif error == 'warn': 

587 logger.warning(message) 

588 else: 

589 assert False 

590 

591 

592class Stopwatch(object): 

593 ''' 

594 Simple stopwatch to measure elapsed wall clock time. 

595 

596 Usage:: 

597 

598 s = Stopwatch() 

599 time.sleep(1) 

600 print s() 

601 time.sleep(1) 

602 print s() 

603 ''' 

604 

605 def __init__(self): 

606 self.start = time.time() 

607 

608 def __call__(self): 

609 return time.time() - self.start 

610 

611 

612def wrap(text, line_length=80): 

613 ''' 

614 Paragraph and list-aware wrapping of text. 

615 ''' 

616 

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

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

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

620 

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

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

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

624 

625 listindents[0:0] = [None] 

626 listindents.append(True) 

627 newlist.append(None) 

628 

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

630 

631 outlines = [] 

632 for ip, p in enumerate(paragraphs): 

633 if not p: 

634 continue 

635 

636 if listindents[ip] is None: 

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

638 findent = _indent 

639 else: 

640 findent = listindents[ip] 

641 _indent = ' ' * len(findent) 

642 

643 ll = line_length - len(_indent) 

644 llf = ll 

645 

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

647 p1 = ' '.join(oldlines) 

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

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

650 parout = match.group(1) 

651 if imatch == 0: 

652 outlines.append(findent + parout) 

653 else: 

654 outlines.append(_indent + parout) 

655 

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

657 listindents[ip] is None 

658 or newlist[ip] is not None 

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

660 

661 outlines.append('') 

662 

663 return outlines 

664 

665 

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

667 lines = list(lines) 

668 if not lines: 

669 return '' 

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

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

672 i = 0 

673 rows = [] 

674 while i < len(lines): 

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

676 i += nx 

677 

678 return '\n'.join(rows) 

679 

680 

681class BetterHelpFormatter(optparse.IndentedHelpFormatter): 

682 

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

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

685 

686 def format_option(self, option): 

687 ''' 

688 From IndentedHelpFormatter but using a different wrap method. 

689 ''' 

690 

691 help_text_position = 4 + self.current_indent 

692 help_text_width = self.width - help_text_position 

693 

694 opts = self.option_strings[option] 

695 opts = "%*s%s" % (self.current_indent, "", opts) 

696 if option.help: 

697 help_text = self.expand_default(option) 

698 

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

700 lines = [ 

701 '', 

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

703 ''] 

704 else: 

705 lines = [''] 

706 lines.append(opts) 

707 lines.append('') 

708 if option.help: 

709 help_lines = wrap(help_text, help_text_width) 

710 lines.extend(["%*s%s" % (help_text_position, "", line) 

711 for line in help_lines]) 

712 lines.append('') 

713 

714 return "\n".join(lines) 

715 

716 def format_description(self, description): 

717 if not description: 

718 return '' 

719 

720 if self.current_indent == 0: 

721 lines = [] 

722 else: 

723 lines = [''] 

724 

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

726 if self.current_indent == 0: 

727 lines.append('\n') 

728 

729 return '\n'.join( 

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

731 

732 

733class ProgressBar: 

734 def __init__(self, label, n): 

735 from pyrocko.progress import progress 

736 self._context = progress.view() 

737 self._context.__enter__() 

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

739 

740 def update(self, i): 

741 self._task.update(i) 

742 

743 def finish(self): 

744 self._task.done() 

745 if self._context: 

746 self._context.__exit__() 

747 self._context = None 

748 

749 

750def progressbar(label, maxval): 

751 if force_dummy_progressbar: 

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

753 

754 return ProgressBar(label, maxval) 

755 

756 

757def progress_beg(label): 

758 ''' 

759 Notify user that an operation has started. 

760 

761 :param label: name of the operation 

762 

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

764 ''' 

765 

766 sys.stderr.write(label) 

767 sys.stderr.flush() 

768 

769 

770def progress_end(label=''): 

771 ''' 

772 Notify user that an operation has ended. 

773 

774 :param label: name of the operation 

775 

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

777 ''' 

778 

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

780 sys.stderr.flush() 

781 

782 

783class ArangeError(Exception): 

784 pass 

785 

786 

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

788 ''' 

789 Return evenly spaced numbers over a specified interval. 

790 

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

792 default and with defined behaviour when stepsize is inconsistent with 

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

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

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

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

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

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

799 multiple of ``step``, respectively. 

800 ''' 

801 

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

803 

804 start = dtype(start) 

805 stop = dtype(stop) 

806 step = dtype(step) 

807 

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

809 

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

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

812 

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

814 raise ArangeError( 

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

816 % (start, stop, step)) 

817 

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

819 x *= step 

820 x += start 

821 return x 

822 

823 

824def polylinefit(x, y, n_or_xnodes): 

825 ''' 

826 Fit piece-wise linear function to data. 

827 

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

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

830 

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

832 polyline, root-mean-square error 

833 ''' 

834 

835 x = num.asarray(x) 

836 y = num.asarray(y) 

837 

838 if isinstance(n_or_xnodes, int): 

839 n = n_or_xnodes 

840 xmin = x.min() 

841 xmax = x.max() 

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

843 else: 

844 xnodes = num.asarray(n_or_xnodes) 

845 n = xnodes.size - 1 

846 

847 assert len(x) == len(y) 

848 assert n > 0 

849 

850 ndata = len(x) 

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

852 for i in range(n): 

853 xmin_block = xnodes[i] 

854 xmax_block = xnodes[i+1] 

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

856 indices = num.where( 

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

858 else: 

859 indices = num.where( 

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

861 

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

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

864 

865 w = float(ndata)*100. 

866 if i < n-1: 

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

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

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

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

871 

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

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

874 

875 ynodes = num.zeros(n+1) 

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

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

878 ynodes[1:n] *= 0.5 

879 

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

881 

882 return xnodes, ynodes, rms_error 

883 

884 

885def plf_integrate_piecewise(x_edges, x, y): 

886 ''' 

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

888 

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

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

891 be sorted. 

892 

893 :param x_edges: array with edges of the intervals 

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

895 control points 

896 ''' 

897 

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

899 ii = num.argsort(x_all) 

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

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

902 xs = x_all[ii] 

903 ys = y_all[ii] 

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

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

906 

907 

908def diff_fd_1d_4o(dt, data): 

909 ''' 

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

911 

912 :param dt: sampling interval 

913 :param data: NumPy array with data samples 

914 

915 :returns: NumPy array with same shape as input 

916 

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

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

919 order central. 

920 ''' 

921 import scipy.signal 

922 

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

924 

925 if data.size >= 5: 

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

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

928 

929 if data.size >= 3: 

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

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

932 

933 if data.size >= 2: 

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

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

936 

937 if data.size == 1: 

938 ddata[0] = 0.0 

939 

940 return ddata 

941 

942 

943def diff_fd_1d_2o(dt, data): 

944 ''' 

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

946 

947 :param dt: sampling interval 

948 :param data: NumPy array with data samples 

949 

950 :returns: NumPy array with same shape as input 

951 

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

953 order right- or left-sided respectively. 

954 

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

956 ''' 

957 

958 return num.gradient(data, dt) 

959 

960 

961def diff_fd_2d_4o(dt, data): 

962 ''' 

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

964 

965 :param dt: sampling interval 

966 :param data: NumPy array with data samples 

967 

968 :returns: NumPy array with same shape as input 

969 

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

971 second order, edge points repeated. 

972 ''' 

973 import scipy.signal 

974 

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

976 

977 if data.size >= 5: 

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

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

980 

981 if data.size >= 3: 

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

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

984 

985 if data.size < 3: 

986 ddata[:] = 0.0 

987 

988 return ddata 

989 

990 

991def diff_fd_2d_2o(dt, data): 

992 ''' 

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

994 

995 :param dt: sampling interval 

996 :param data: NumPy array with data samples 

997 

998 :returns: NumPy array with same shape as input 

999 

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

1001 ''' 

1002 import scipy.signal 

1003 

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

1005 

1006 if data.size >= 3: 

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

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

1009 

1010 ddata[0] = ddata[1] 

1011 ddata[-1] = ddata[-2] 

1012 

1013 if data.size < 3: 

1014 ddata[:] = 0.0 

1015 

1016 return ddata 

1017 

1018 

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

1020 ''' 

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

1022 

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

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

1025 :param dt: sampling interval 

1026 :param data: NumPy array with data samples 

1027 

1028 :returns: NumPy array with same shape as input 

1029 

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

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

1032 :py:func:`diff_fd_2d_4o`. 

1033 

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

1035 ''' 

1036 

1037 funcs = { 

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

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

1040 

1041 try: 

1042 funcs_n = funcs[n] 

1043 except KeyError: 

1044 raise ValueError( 

1045 'pyrocko.util.diff_fd: ' 

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

1047 

1048 try: 

1049 func = funcs_n[order] 

1050 except KeyError: 

1051 raise ValueError( 

1052 'pyrocko.util.diff_fd: ' 

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

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

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

1056 

1057 return func(dt, data) 

1058 

1059 

1060class GlobalVars(object): 

1061 reuse_store = dict() 

1062 decitab_nmax = 0 

1063 decitab = {} 

1064 decimate_fir_coeffs = {} 

1065 decimate_fir_remez_coeffs = {} 

1066 decimate_iir_coeffs = {} 

1067 re_frac = None 

1068 

1069 

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

1071 

1072 q = int(q) 

1073 

1074 if n is None: 

1075 if ftype == 'fir': 

1076 n = 30 

1077 elif ftype == 'fir-remez': 

1078 n = 45*q 

1079 else: 

1080 n = 8 

1081 

1082 if ftype == 'fir': 

1083 coeffs = GlobalVars.decimate_fir_coeffs 

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

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

1086 

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

1088 return b, [1.], n 

1089 

1090 elif ftype == 'fir-remez': 

1091 coeffs = GlobalVars.decimate_fir_remez_coeffs 

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

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

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

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

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

1097 return b, [1.], n 

1098 

1099 else: 

1100 coeffs = GlobalVars.decimate_iir_coeffs 

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

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

1103 

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

1105 return b, a, n 

1106 

1107 

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

1109 ''' 

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

1111 

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

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

1114 

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

1116 :param q: the downsampling factor 

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

1118 `fir` filter) 

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

1120 

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

1122 

1123 ''' 

1124 

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

1126 

1127 if zi is None or zi is True: 

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

1129 else: 

1130 zi_ = zi 

1131 

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

1133 

1134 if zi is not None: 

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

1136 else: 

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

1138 

1139 

1140class UnavailableDecimation(Exception): 

1141 ''' 

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

1143 ''' 

1144 

1145 pass 

1146 

1147 

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

1149 ''' 

1150 Greatest common divisor. 

1151 ''' 

1152 

1153 while b > epsilon*a: 

1154 a, b = b, a % b 

1155 

1156 return a 

1157 

1158 

1159def lcm(a, b): 

1160 ''' 

1161 Least common multiple. 

1162 ''' 

1163 

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

1165 

1166 

1167def mk_decitab(nmax=100): 

1168 ''' 

1169 Make table with decimation sequences. 

1170 

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

1172 successive application of :py:func:`decimation` with small integer 

1173 downsampling factors (because using large downampling factors can make the 

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

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

1176 ''' 

1177 

1178 tab = GlobalVars.decitab 

1179 for i in range(1, 10): 

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

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

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

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

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

1185 if p > nmax: 

1186 break 

1187 if p not in tab: 

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

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

1190 break 

1191 if i*j*k > nmax: 

1192 break 

1193 if i*j > nmax: 

1194 break 

1195 if i > nmax: 

1196 break 

1197 

1198 GlobalVars.decitab_nmax = nmax 

1199 

1200 

1201def zfmt(n): 

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

1203 

1204 

1205def _year_to_time(year): 

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

1207 return to_time_float(calendar.timegm(tt)) 

1208 

1209 

1210def _working_year(year): 

1211 try: 

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

1213 t = calendar.timegm(tt) 

1214 tt2_ = time.gmtime(t) 

1215 tt2 = tuple(tt2_)[:6] 

1216 if tt != tt2: 

1217 return False 

1218 

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

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

1221 if s != s2: 

1222 return False 

1223 

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

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

1226 if s3 != s2: 

1227 return False 

1228 

1229 except Exception: 

1230 return False 

1231 

1232 return True 

1233 

1234 

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

1236 ''' 

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

1238 

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

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

1241 which is fully supported. 

1242 

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

1244 ''' 

1245 

1246 year0 = 2000 

1247 year_min = year0 

1248 year_max = year0 

1249 

1250 itests = list(range(101)) 

1251 for i in range(19): 

1252 itests.append(200 + i*100) 

1253 

1254 for i in itests: 

1255 year = year0 - i 

1256 if year_min_lim is not None and year < year_min_lim: 

1257 year_min = year_min_lim 

1258 break 

1259 elif not _working_year(year): 

1260 break 

1261 else: 

1262 year_min = year 

1263 

1264 for i in itests: 

1265 year = year0 + i 

1266 if year_max_lim is not None and year > year_max_lim: 

1267 year_max = year_max_lim 

1268 break 

1269 elif not _working_year(year + 1): 

1270 break 

1271 else: 

1272 year_max = year 

1273 

1274 return ( 

1275 _year_to_time(year_min), 

1276 _year_to_time(year_max), 

1277 year_min, year_max) 

1278 

1279 

1280g_working_system_time_range = None 

1281 

1282 

1283def get_working_system_time_range(): 

1284 ''' 

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

1286 ''' 

1287 

1288 global g_working_system_time_range 

1289 if g_working_system_time_range is None: 

1290 g_working_system_time_range = working_system_time_range() 

1291 

1292 return g_working_system_time_range 

1293 

1294 

1295def is_working_time(t): 

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

1297 return tmin <= t <= tmax 

1298 

1299 

1300def julian_day_of_year(timestamp): 

1301 ''' 

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

1303 

1304 :returns: day number as int 

1305 ''' 

1306 

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

1308 

1309 

1310def hour_start(timestamp): 

1311 ''' 

1312 Get beginning of hour for any point in time. 

1313 

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

1315 

1316 :returns: instant of hour start as system timestamp 

1317 ''' 

1318 

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

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

1321 return to_time_float(calendar.timegm(tts)) 

1322 

1323 

1324def day_start(timestamp): 

1325 ''' 

1326 Get beginning of day for any point in time. 

1327 

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

1329 

1330 :returns: instant of day start as system timestamp 

1331 ''' 

1332 

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

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

1335 return to_time_float(calendar.timegm(tts)) 

1336 

1337 

1338def month_start(timestamp): 

1339 ''' 

1340 Get beginning of month for any point in time. 

1341 

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

1343 

1344 :returns: instant of month start as system timestamp 

1345 ''' 

1346 

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

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

1349 return to_time_float(calendar.timegm(tts)) 

1350 

1351 

1352def year_start(timestamp): 

1353 ''' 

1354 Get beginning of year for any point in time. 

1355 

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

1357 

1358 :returns: instant of year start as system timestamp 

1359 ''' 

1360 

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

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

1363 return to_time_float(calendar.timegm(tts)) 

1364 

1365 

1366def iter_days(tmin, tmax): 

1367 ''' 

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

1369 

1370 :param tmin,tmax: input time span 

1371 

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

1373 ''' 

1374 

1375 t = day_start(tmin) 

1376 while t < tmax: 

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

1378 yield t, tend 

1379 t = tend 

1380 

1381 

1382def iter_months(tmin, tmax): 

1383 ''' 

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

1385 

1386 :param tmin,tmax: input time span 

1387 

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

1389 ''' 

1390 

1391 t = month_start(tmin) 

1392 while t < tmax: 

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

1394 yield t, tend 

1395 t = tend 

1396 

1397 

1398def iter_years(tmin, tmax): 

1399 ''' 

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

1401 

1402 :param tmin,tmax: input time span 

1403 

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

1405 ''' 

1406 

1407 t = year_start(tmin) 

1408 while t < tmax: 

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

1410 yield t, tend 

1411 t = tend 

1412 

1413 

1414def today(): 

1415 return day_start(time.time()) 

1416 

1417 

1418def tomorrow(): 

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

1420 

1421 

1422def decitab(n): 

1423 ''' 

1424 Get integer decimation sequence for given downampling factor. 

1425 

1426 :param n: target decimation factor 

1427 

1428 :returns: tuple with downsampling sequence 

1429 ''' 

1430 

1431 if n > GlobalVars.decitab_nmax: 

1432 mk_decitab(n*2) 

1433 if n not in GlobalVars.decitab: 

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

1435 return GlobalVars.decitab[n] 

1436 

1437 

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

1439 ''' 

1440 Convert string representing UTC time to system time. 

1441 

1442 :param s: string to be interpreted 

1443 :param format: format string passed to :py:func:`strptime` 

1444 

1445 :returns: system time stamp 

1446 

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

1448 

1449 .. note:: 

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

1451 ''' 

1452 

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

1454 

1455 

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

1457 ''' 

1458 Get string representation from system time, UTC. 

1459 

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

1461 

1462 .. note:: 

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

1464 ''' 

1465 

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

1467 

1468 

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

1470 ''' 

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

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

1473 

1474 .. note:: 

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

1476 ''' 

1477 

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

1479 

1480 

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

1482 ''' 

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

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

1485 

1486 .. note:: 

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

1488 ''' 

1489 

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

1491 

1492 

1493class TimeStrError(Exception): 

1494 pass 

1495 

1496 

1497class FractionalSecondsMissing(TimeStrError): 

1498 ''' 

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

1500 fractional seconds. 

1501 ''' 

1502 

1503 pass 

1504 

1505 

1506class FractionalSecondsWrongNumberOfDigits(TimeStrError): 

1507 ''' 

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

1509 incorrect number of digits in the fractional seconds part. 

1510 ''' 

1511 

1512 pass 

1513 

1514 

1515def _endswith_n(s, endings): 

1516 for ix, x in enumerate(endings): 

1517 if s.endswith(x): 

1518 return ix 

1519 return -1 

1520 

1521 

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

1523 ''' 

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

1525 

1526 :param s: string representing UTC time 

1527 :param format: time string format 

1528 :returns: system time stamp as floating point value 

1529 

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

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

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

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

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

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

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

1537 seconds. 

1538 ''' 

1539 

1540 time_float = get_time_float() 

1541 

1542 if util_ext is not None: 

1543 try: 

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

1545 except util_ext.UtilExtError as e: 

1546 raise TimeStrError( 

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

1548 

1549 return time_float(t)+tfrac 

1550 

1551 fracsec = 0. 

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

1553 

1554 iend = _endswith_n(format, fixed_endings) 

1555 if iend != -1: 

1556 dotpos = s.rfind('.') 

1557 if dotpos == -1: 

1558 raise FractionalSecondsMissing( 

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

1560 

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

1562 raise FractionalSecondsWrongNumberOfDigits( 

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

1564 

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

1566 fracsec = float(s[dotpos:]) 

1567 s = s[:dotpos] 

1568 

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

1570 dotpos = s.rfind('.') 

1571 format = format[:-8] 

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

1573 fracsec = float(s[dotpos:]) 

1574 

1575 if dotpos != -1: 

1576 s = s[:dotpos] 

1577 

1578 try: 

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

1580 + fracsec 

1581 except ValueError as e: 

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

1583 

1584 

1585stt = str_to_time 

1586 

1587 

1588def str_to_time_fillup(s): 

1589 ''' 

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

1591 

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

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

1594 ''' 

1595 

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

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

1598 

1599 return str_to_time(s) 

1600 

1601 

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

1603 ''' 

1604 Get string representation for floating point system time. 

1605 

1606 :param t: floating point system time 

1607 :param format: time string format 

1608 :returns: string representing UTC time 

1609 

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

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

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

1613 with ``x`` digits precision. 

1614 ''' 

1615 

1616 if pyrocko.grumpy > 0: 

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

1618 

1619 if isinstance(format, int): 

1620 if format > 0: 

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

1622 else: 

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

1624 

1625 if util_ext is not None: 

1626 t0 = num.floor(t) 

1627 try: 

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

1629 except util_ext.UtilExtError as e: 

1630 raise TimeStrError( 

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

1632 

1633 if not GlobalVars.re_frac: 

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

1635 GlobalVars.frac_formats = dict( 

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

1637 

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

1639 tfrac = t-ts 

1640 

1641 m = GlobalVars.re_frac.search(format) 

1642 if m: 

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

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

1645 ts += 1. 

1646 

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

1648 

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

1650 

1651 

1652tts = time_to_str 

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

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

1655 

1656 

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

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

1659 # should be used. 

1660 

1661 if fmt is None: 

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

1663 

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

1665 # Setting the locale seems to have no effect. 

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

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

1668 

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

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

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

1672 return time.strftime(fmt, tt) 

1673 

1674 

1675def gmtime_x(timestamp): 

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

1677 tt = time.gmtime(etimestamp) 

1678 ms = (timestamp-etimestamp)*1000 

1679 return tt, ms 

1680 

1681 

1682def plural_s(n): 

1683 if not isinstance(n, int): 

1684 n = len(n) 

1685 

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

1687 

1688 

1689def ensuredirs(dst): 

1690 ''' 

1691 Create all intermediate path components for a target path. 

1692 

1693 :param dst: target path 

1694 

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

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

1697 ''' 

1698 

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

1700 dirs = [] 

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

1702 dirs.append(d) 

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

1704 

1705 dirs.reverse() 

1706 

1707 for d in dirs: 

1708 try: 

1709 os.mkdir(d) 

1710 except OSError as e: 

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

1712 raise 

1713 

1714 

1715def ensuredir(dst): 

1716 ''' 

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

1718 

1719 :param dst: directory name 

1720 

1721 Nothing is done if the given target already exists. 

1722 ''' 

1723 

1724 if os.path.exists(dst): 

1725 return 

1726 

1727 dst.rstrip(os.sep) 

1728 

1729 ensuredirs(dst) 

1730 try: 

1731 os.mkdir(dst) 

1732 except OSError as e: 

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

1734 raise 

1735 

1736 

1737def reuse(x): 

1738 ''' 

1739 Get unique instance of an object. 

1740 

1741 :param x: hashable object 

1742 :returns: reference to x or an equivalent object 

1743 

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

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

1746 ''' 

1747 

1748 grs = GlobalVars.reuse_store 

1749 if x not in grs: 

1750 grs[x] = x 

1751 return grs[x] 

1752 

1753 

1754def deuse(x): 

1755 grs = GlobalVars.reuse_store 

1756 if x in grs: 

1757 del grs[x] 

1758 

1759 

1760class Anon(object): 

1761 ''' 

1762 Dict-to-object utility. 

1763 

1764 Any given arguments are stored as attributes. 

1765 

1766 Example:: 

1767 

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

1769 print a.x, a.y 

1770 ''' 

1771 

1772 def __init__(self, **dict): 

1773 for k in dict: 

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

1775 

1776 

1777def iter_select_files( 

1778 paths, 

1779 include=None, 

1780 exclude=None, 

1781 selector=None, 

1782 show_progress=True, 

1783 pass_through=None): 

1784 

1785 ''' 

1786 Recursively select files (generator variant). 

1787 

1788 See :py:func:`select_files`. 

1789 ''' 

1790 

1791 if show_progress: 

1792 progress_beg('selecting files...') 

1793 

1794 ngood = 0 

1795 check_include = None 

1796 if include is not None: 

1797 rinclude = re.compile(include) 

1798 

1799 def check_include(path): 

1800 m = rinclude.search(path) 

1801 if not m: 

1802 return False 

1803 

1804 if selector is None: 

1805 return True 

1806 

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

1808 return selector(infos) 

1809 

1810 check_exclude = None 

1811 if exclude is not None: 

1812 rexclude = re.compile(exclude) 

1813 

1814 def check_exclude(path): 

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

1816 

1817 if check_include and check_exclude: 

1818 

1819 def check(path): 

1820 return check_include(path) and check_exclude(path) 

1821 

1822 elif check_include: 

1823 check = check_include 

1824 

1825 elif check_exclude: 

1826 check = check_exclude 

1827 

1828 else: 

1829 check = None 

1830 

1831 if isinstance(paths, str): 

1832 paths = [paths] 

1833 

1834 for path in paths: 

1835 if pass_through and pass_through(path): 

1836 if check is None or check(path): 

1837 yield path 

1838 

1839 elif os.path.isdir(path): 

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

1841 dirnames.sort() 

1842 filenames.sort() 

1843 for filename in filenames: 

1844 path = op.join(dirpath, filename) 

1845 if check is None or check(path): 

1846 yield os.path.abspath(path) 

1847 ngood += 1 

1848 else: 

1849 if check is None or check(path): 

1850 yield os.path.abspath(path) 

1851 ngood += 1 

1852 

1853 if show_progress: 

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

1855 

1856 

1857def select_files( 

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

1859 regex=None): 

1860 

1861 ''' 

1862 Recursively select files. 

1863 

1864 :param paths: entry path names 

1865 :param include: pattern for conditional inclusion 

1866 :param exclude: pattern for conditional exclusion 

1867 :param selector: callback for conditional inclusion 

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

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

1870 :returns: list of path names 

1871 

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

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

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

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

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

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

1878 groups given in ``include``. 

1879 

1880 Examples 

1881 

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

1883 

1884 select_files(paths, 

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

1886 

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

1888 the year:: 

1889 

1890 select_files(paths, 

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

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

1893 ''' 

1894 if regex is not None: 

1895 assert include is None 

1896 include = regex 

1897 

1898 return list(iter_select_files( 

1899 paths, include, exclude, selector, show_progress)) 

1900 

1901 

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

1903 ''' 

1904 Convert positive integer to a base36 string. 

1905 ''' 

1906 

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

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

1909 if number < 0: 

1910 raise ValueError('number must be positive') 

1911 

1912 # Special case for small numbers 

1913 if number < 36: 

1914 return alphabet[number] 

1915 

1916 base36 = '' 

1917 while number != 0: 

1918 number, i = divmod(number, 36) 

1919 base36 = alphabet[i] + base36 

1920 

1921 return base36 

1922 

1923 

1924def base36decode(number): 

1925 ''' 

1926 Decode base36 endcoded positive integer. 

1927 ''' 

1928 

1929 return int(number, 36) 

1930 

1931 

1932class UnpackError(Exception): 

1933 ''' 

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

1935 ''' 

1936 

1937 pass 

1938 

1939 

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

1941 + '\n' + '0123456789' * 8 + '\n' 

1942 

1943 

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

1945 ''' 

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

1947 

1948 :param format: format specification 

1949 :param line: string to be processed 

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

1951 

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

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

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

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

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

1957 

1958 The following field types are available: 

1959 

1960 ==== ================================================================ 

1961 Type Description 

1962 ==== ================================================================ 

1963 A string (full field width is extracted) 

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

1965 i integer value 

1966 f floating point value 

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

1968 x special field type to skip parts of the string 

1969 ==== ================================================================ 

1970 ''' 

1971 

1972 ipos = 0 

1973 values = [] 

1974 icall = 0 

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

1976 form = form.strip() 

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

1978 form = form.rstrip('?') 

1979 typ = form[0] 

1980 ln = int(form[1:]) 

1981 s = line[ipos:ipos+ln] 

1982 cast = { 

1983 'x': None, 

1984 'A': str, 

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

1986 'i': int, 

1987 'f': float, 

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

1989 

1990 if cast == 'extra': 

1991 cast = callargs[icall] 

1992 icall += 1 

1993 

1994 if cast is not None: 

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

1996 values.append(None) 

1997 else: 

1998 try: 

1999 values.append(cast(s)) 

2000 except Exception: 

2001 mark = [' '] * 80 

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

2003 mark = ''.join(mark) 

2004 raise UnpackError( 

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

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

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

2008 

2009 ipos += ln 

2010 

2011 return values 

2012 

2013 

2014_pattern_cache = {} 

2015 

2016 

2017def _nslc_pattern(pattern): 

2018 if pattern not in _pattern_cache: 

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

2020 _pattern_cache[pattern] = rpattern 

2021 else: 

2022 rpattern = _pattern_cache[pattern] 

2023 

2024 return rpattern 

2025 

2026 

2027def match_nslc(patterns, nslc): 

2028 ''' 

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

2030 patterns. 

2031 

2032 :param patterns: pattern or list of patterns 

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

2034 

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

2036 match; or ``False``. 

2037 

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

2039 

2040 Example:: 

2041 

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

2043 ''' 

2044 

2045 if isinstance(patterns, str): 

2046 patterns = [patterns] 

2047 

2048 if not isinstance(nslc, str): 

2049 s = '.'.join(nslc) 

2050 else: 

2051 s = nslc 

2052 

2053 for pattern in patterns: 

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

2055 return True 

2056 

2057 return False 

2058 

2059 

2060def match_nslcs(patterns, nslcs): 

2061 ''' 

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

2063 of several given patterns. 

2064 

2065 :param patterns: pattern or list of patterns 

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

2067 

2068 See also :py:func:`match_nslc` 

2069 ''' 

2070 

2071 matching = [] 

2072 for nslc in nslcs: 

2073 if match_nslc(patterns, nslc): 

2074 matching.append(nslc) 

2075 

2076 return matching 

2077 

2078 

2079class Timeout(Exception): 

2080 pass 

2081 

2082 

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

2084 t0 = time.time() 

2085 

2086 while True: 

2087 try: 

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

2089 os.close(f) 

2090 return 

2091 

2092 except OSError as e: 

2093 if e.errno == errno.EEXIST: 

2094 tnow = time.time() 

2095 

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

2097 raise Timeout( 

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

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

2100 

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

2102 logger.warning( 

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

2104 timewarn, fn)) 

2105 

2106 timewarn *= 2 

2107 

2108 time.sleep(0.01) 

2109 else: 

2110 raise 

2111 

2112 

2113def delete_lockfile(fn): 

2114 os.unlink(fn) 

2115 

2116 

2117class Lockfile(Exception): 

2118 

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

2120 self._path = path 

2121 self._timeout = timeout 

2122 self._timewarn = timewarn 

2123 

2124 def __enter__(self): 

2125 create_lockfile( 

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

2127 return None 

2128 

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

2130 delete_lockfile(self._path) 

2131 

2132 

2133class SoleError(Exception): 

2134 ''' 

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

2136 instance is running. 

2137 ''' 

2138 

2139 pass 

2140 

2141 

2142class Sole(object): 

2143 ''' 

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

2145 program is running. 

2146 

2147 :param pid_path: path to lockfile to be used 

2148 

2149 Usage:: 

2150 

2151 from pyrocko.util import Sole, SoleError, setup_logging 

2152 import os 

2153 

2154 setup_logging('my_program') 

2155 

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

2157 try: 

2158 sole = Sole(pid_path) 

2159 

2160 except SoleError, e: 

2161 logger.fatal( str(e) ) 

2162 sys.exit(1) 

2163 ''' 

2164 

2165 def __init__(self, pid_path): 

2166 self._pid_path = pid_path 

2167 self._other_running = False 

2168 ensuredirs(self._pid_path) 

2169 self._lockfile = None 

2170 self._os = os 

2171 

2172 if not fcntl: 

2173 raise SoleError( 

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

2175 'this platform.') 

2176 

2177 self._fcntl = fcntl 

2178 

2179 try: 

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

2181 except Exception: 

2182 raise SoleError( 

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

2184 

2185 try: 

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

2187 

2188 except IOError: 

2189 self._other_running = True 

2190 try: 

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

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

2193 f.close() 

2194 except Exception: 

2195 pid = '?' 

2196 

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

2198 

2199 try: 

2200 os.ftruncate(self._lockfile, 0) 

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

2202 os.fsync(self._lockfile) 

2203 

2204 except Exception: 

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

2206 # to fail 

2207 pass 

2208 

2209 def __del__(self): 

2210 if not self._other_running: 

2211 if self._lockfile is not None: 

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

2213 self._os.close(self._lockfile) 

2214 try: 

2215 self._os.unlink(self._pid_path) 

2216 except Exception: 

2217 pass 

2218 

2219 

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

2221 

2222 

2223def escapequotes(s): 

2224 return re_escapequotes.sub(r"\\\1", s) 

2225 

2226 

2227class TableWriter(object): 

2228 ''' 

2229 Write table of space separated values to a file. 

2230 

2231 :param f: file like object 

2232 

2233 Strings containing spaces are quoted on output. 

2234 ''' 

2235 

2236 def __init__(self, f): 

2237 self._f = f 

2238 

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

2240 

2241 ''' 

2242 Write one row of values to underlying file. 

2243 

2244 :param row: iterable of values 

2245 :param minfieldwidths: minimum field widths for the values 

2246 

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

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

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

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

2251 backslash-escaped. 

2252 ''' 

2253 

2254 out = [] 

2255 

2256 for i, x in enumerate(row): 

2257 w = 0 

2258 if minfieldwidths and i < len(minfieldwidths): 

2259 w = minfieldwidths[i] 

2260 

2261 if isinstance(x, str): 

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

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

2264 

2265 x = x.ljust(w) 

2266 else: 

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

2268 

2269 out.append(x) 

2270 

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

2272 

2273 

2274class TableReader(object): 

2275 

2276 ''' 

2277 Read table of space separated values from a file. 

2278 

2279 :param f: file-like object 

2280 

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

2282 with quoted strings. 

2283 ''' 

2284 

2285 def __init__(self, f): 

2286 self._f = f 

2287 self.eof = False 

2288 

2289 def readrow(self): 

2290 ''' 

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

2292 

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

2294 ''' 

2295 

2296 line = self._f.readline() 

2297 if not line: 

2298 self.eof = True 

2299 return [] 

2300 line.strip() 

2301 if line.startswith('#'): 

2302 return [] 

2303 

2304 return qsplit(line) 

2305 

2306 

2307def gform(number, significant_digits=3): 

2308 ''' 

2309 Pretty print floating point numbers. 

2310 

2311 Align floating point numbers at the decimal dot. 

2312 

2313 :: 

2314 

2315 | -d.dde+xxx| 

2316 | -d.dde+xx | 

2317 |-ddd. | 

2318 | -dd.d | 

2319 | -d.dd | 

2320 | -0.ddd | 

2321 | -0.0ddd | 

2322 | -0.00ddd | 

2323 | -d.dde-xx | 

2324 | -d.dde-xxx| 

2325 | nan| 

2326 

2327 

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

2329 ''' 

2330 

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

2332 pow(10., significant_digits)) 

2333 width = significant_digits+significant_digits-1+1+1+5 

2334 

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

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

2337 else: 

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

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

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

2341 s = 'nan'.rjust(width) 

2342 return s 

2343 

2344 

2345def human_bytesize(value): 

2346 

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

2348 

2349 if value == 1: 

2350 return '1 Byte' 

2351 

2352 for i, ext in enumerate(exts): 

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

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

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

2356 if round(x) < 1000.: 

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

2358 

2359 return '%i Bytes' % value 

2360 

2361 

2362re_compatibility = re.compile( 

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

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

2365) 

2366 

2367 

2368def pf_is_old(fn): 

2369 oldstyle = False 

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

2371 for line in f: 

2372 if re_compatibility.search(line): 

2373 oldstyle = True 

2374 

2375 return oldstyle 

2376 

2377 

2378def pf_upgrade(fn): 

2379 need = pf_is_old(fn) 

2380 if need: 

2381 fn_temp = fn + '.temp' 

2382 

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

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

2385 for line in fin: 

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

2387 fout.write(line) 

2388 

2389 os.rename(fn_temp, fn) 

2390 

2391 return need 

2392 

2393 

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

2395 ''' 

2396 Extract leap second information from tzdata. 

2397 

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

2399 extract-historic-leap-seconds-from-tzdata 

2400 

2401 See also 'man 5 tzfile'. 

2402 ''' 

2403 

2404 from struct import unpack, calcsize 

2405 out = [] 

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

2407 # read header 

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

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

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

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

2412 

2413 # skip over some uninteresting data 

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

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

2416 f.read(calcsize(fmt)) 

2417 

2418 # read leap-seconds 

2419 fmt = '>2l' 

2420 for i in range(leapcnt): 

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

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

2423 

2424 return out 

2425 

2426 

2427class LeapSecondsError(Exception): 

2428 pass 

2429 

2430 

2431class LeapSecondsOutdated(LeapSecondsError): 

2432 pass 

2433 

2434 

2435class InvalidLeapSecondsFile(LeapSecondsOutdated): 

2436 pass 

2437 

2438 

2439def parse_leap_seconds_list(fn): 

2440 data = [] 

2441 texpires = None 

2442 try: 

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

2444 except TimeStrError: 

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

2446 

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

2448 

2449 if not op.exists(fn): 

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

2451 

2452 try: 

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

2454 for line in f: 

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

2456 raise InvalidLeapSecondsFile('invalid leap seconds file') 

2457 

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

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

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

2461 pass 

2462 else: 

2463 toks = line.split() 

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

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

2466 data.append((t, nleap)) 

2467 

2468 except IOError: 

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

2470 

2471 if texpires is None or tnow > texpires: 

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

2473 

2474 return data 

2475 

2476 

2477def read_leap_seconds2(): 

2478 from pyrocko import config 

2479 conf = config.config() 

2480 fn = conf.leapseconds_path 

2481 url = conf.leapseconds_url 

2482 # check for outdated default URL 

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

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

2485 logger.info( 

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

2487 

2488 for i in range(3): 

2489 try: 

2490 return parse_leap_seconds_list(fn) 

2491 

2492 except LeapSecondsOutdated: 

2493 try: 

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

2495 download_file(url, fn) 

2496 

2497 except Exception as e: 

2498 raise LeapSecondsError( 

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

2500 % (url, fn, e)) 

2501 

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

2503 

2504 

2505def gps_utc_offset(t_utc): 

2506 ''' 

2507 Time offset t_gps - t_utc for a given t_utc. 

2508 ''' 

2509 ls = read_leap_seconds2() 

2510 i = 0 

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

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

2513 

2514 while i < len(ls) - 1: 

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

2516 return ls[i][1] - 9 

2517 i += 1 

2518 

2519 return ls[-1][1] - 9 

2520 

2521 

2522def utc_gps_offset(t_gps): 

2523 ''' 

2524 Time offset t_utc - t_gps for a given t_gps. 

2525 ''' 

2526 ls = read_leap_seconds2() 

2527 

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

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

2530 

2531 i = 0 

2532 while i < len(ls) - 1: 

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

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

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

2536 i += 1 

2537 

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

2539 

2540 

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

2542 import itertools 

2543 import glob 

2544 from pyrocko.io.io_common import FileLoadError 

2545 

2546 def iload_filename(filename, **kwargs): 

2547 try: 

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

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

2550 yield cr 

2551 

2552 except FileLoadError as e: 

2553 e.set_context('filename', filename) 

2554 raise 

2555 

2556 def iload_dirname(dirname, **kwargs): 

2557 for entry in os.listdir(dirname): 

2558 fpath = op.join(dirname, entry) 

2559 if op.isfile(fpath): 

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

2561 yield cr 

2562 

2563 def iload_glob(pattern, **kwargs): 

2564 

2565 for fn in glob.iglob(pattern): 

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

2567 yield cr 

2568 

2569 def iload(source, **kwargs): 

2570 if isinstance(source, str): 

2571 if op.isdir(source): 

2572 return iload_dirname(source, **kwargs) 

2573 elif op.isfile(source): 

2574 return iload_filename(source, **kwargs) 

2575 else: 

2576 return iload_glob(source, **kwargs) 

2577 

2578 elif hasattr(source, 'read'): 

2579 return iload_fh(source, **kwargs) 

2580 else: 

2581 return itertools.chain.from_iterable( 

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

2583 

2584 iload_filename.__doc__ = ''' 

2585 Read %s information from named file. 

2586 ''' % doc_fmt 

2587 

2588 iload_dirname.__doc__ = ''' 

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

2590 ''' % (doc_fmt, doc_fmt) 

2591 

2592 iload_glob.__doc__ = ''' 

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

2594 ''' % doc_fmt 

2595 

2596 iload.__doc__ = ''' 

2597 Load %s information from given source(s) 

2598 

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

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

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

2602 

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

2604 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects) 

2605 

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

2607 f.__module__ = iload_fh.__module__ 

2608 

2609 return iload_filename, iload_dirname, iload_glob, iload 

2610 

2611 

2612class Inconsistency(Exception): 

2613 pass 

2614 

2615 

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

2617 ''' 

2618 Check for inconsistencies. 

2619 

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

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

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

2623 ''' 

2624 

2625 if len(list_of_tuples) >= 2: 

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

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

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

2629 for t in list_of_tuples)) 

2630 

2631 

2632class defaultzerodict(dict): 

2633 def __missing__(self, k): 

2634 return 0 

2635 

2636 

2637def mostfrequent(x): 

2638 c = defaultzerodict() 

2639 for e in x: 

2640 c[e] += 1 

2641 

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

2643 

2644 

2645def consistency_merge(list_of_tuples, 

2646 message='values differ:', 

2647 error='raise', 

2648 merge=mostfrequent): 

2649 

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

2651 

2652 if len(list_of_tuples) == 0: 

2653 raise Exception('cannot merge empty sequence') 

2654 

2655 try: 

2656 consistency_check(list_of_tuples, message) 

2657 return list_of_tuples[0][1:] 

2658 except Inconsistency as e: 

2659 if error == 'raise': 

2660 raise 

2661 

2662 elif error == 'warn': 

2663 logger.warning(str(e)) 

2664 

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

2666 

2667 

2668def short_to_list(nmax, it): 

2669 import itertools 

2670 

2671 if isinstance(it, list): 

2672 return it 

2673 

2674 li = [] 

2675 for i in range(nmax+1): 

2676 try: 

2677 li.append(next(it)) 

2678 except StopIteration: 

2679 return li 

2680 

2681 return itertools.chain(li, it) 

2682 

2683 

2684def parse_md(f): 

2685 try: 

2686 with open(op.join( 

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

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

2689 mdstr = readme.read() 

2690 except IOError as e: 

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

2692 

2693 # Remve the title 

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

2695 # Append sphinx reference to `pyrocko.` modules 

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

2697 # Convert Subsections to toc-less rubrics 

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

2699 return mdstr 

2700 

2701 

2702def mpl_show(plt): 

2703 import matplotlib 

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

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

2706 else: 

2707 plt.show() 

2708 

2709 

2710g_re_qsplit = re.compile( 

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

2712g_re_qsplit_sep = {} 

2713 

2714 

2715def get_re_qsplit(sep): 

2716 if sep is None: 

2717 return g_re_qsplit 

2718 else: 

2719 if sep not in g_re_qsplit_sep: 

2720 assert len(sep) == 1 

2721 assert sep not in '\'"' 

2722 esep = re.escape(sep) 

2723 g_re_qsplit_sep[sep] = re.compile( 

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

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

2726 return g_re_qsplit_sep[sep] 

2727 

2728 

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

2730g_re_trivial_sep = {} 

2731 

2732 

2733def get_re_trivial(sep): 

2734 if sep is None: 

2735 return g_re_trivial 

2736 else: 

2737 if sep not in g_re_qsplit_sep: 

2738 assert len(sep) == 1 

2739 assert sep not in '\'"' 

2740 esep = re.escape(sep) 

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

2742 

2743 return g_re_trivial_sep[sep] 

2744 

2745 

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

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

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

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

2750 

2751 

2752def escape_s(s): 

2753 ''' 

2754 Backslash-escape single-quotes and backslashes. 

2755 

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

2757 

2758 ''' 

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

2760 

2761 

2762def unescape_s(s): 

2763 ''' 

2764 Unescape backslash-escaped single-quotes and backslashes. 

2765 

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

2767 ''' 

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

2769 

2770 

2771def escape_d(s): 

2772 ''' 

2773 Backslash-escape double-quotes and backslashes. 

2774 

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

2776 ''' 

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

2778 

2779 

2780def unescape_d(s): 

2781 ''' 

2782 Unescape backslash-escaped double-quotes and backslashes. 

2783 

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

2785 ''' 

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

2787 

2788 

2789def qjoin_s(it, sep=None): 

2790 ''' 

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

2792 

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

2794 ''' 

2795 re_trivial = get_re_trivial(sep) 

2796 

2797 if sep is None: 

2798 sep = ' ' 

2799 

2800 return sep.join( 

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

2802 

2803 

2804def qjoin_d(it, sep=None): 

2805 ''' 

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

2807 

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

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

2810 ''' 

2811 re_trivial = get_re_trivial(sep) 

2812 if sep is None: 

2813 sep = ' ' 

2814 

2815 return sep.join( 

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

2817 

2818 

2819def qsplit(s, sep=None): 

2820 ''' 

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

2822 

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

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

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

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

2827 ''' 

2828 re_qsplit = get_re_qsplit(sep) 

2829 return [ 

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

2831 for x in re_qsplit.findall(s)] 

2832 

2833 

2834g_have_warned_threadpoolctl = False 

2835 

2836 

2837class threadpool_limits_dummy(object): 

2838 

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

2840 pass 

2841 

2842 def __enter__(self): 

2843 global g_have_warned_threadpoolctl 

2844 

2845 if not g_have_warned_threadpoolctl: 

2846 logger.warning( 

2847 'Cannot control number of BLAS threads because ' 

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

2849 'install `threadpoolctl`.') 

2850 

2851 g_have_warned_threadpoolctl = True 

2852 

2853 return self 

2854 

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

2856 pass 

2857 

2858 

2859def get_threadpool_limits(): 

2860 ''' 

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

2862 ''' 

2863 

2864 try: 

2865 from threadpoolctl import threadpool_limits 

2866 return threadpool_limits 

2867 

2868 except ImportError: 

2869 return threadpool_limits_dummy