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-06 06:59 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5'''
6Utility functions for Pyrocko.
8.. _time-handling-mode:
10High precision time handling mode
11.................................
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.
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.
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.
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:
37.. code-block :: python
39 import time
40 from pyrocko import util
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
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())
50 # To get the appropriate float class, use:
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]
57 # Shortcut:
58 util.check_time_class(t1)
60Module content
61..............
63.. py:class:: hpfloat
65 Alias for NumPy's high precision float data type ``float128`` or
66 ``float96``, if available.
68 On platforms lacking support for high precision floats, an attempt to
69 create a ``hpfloat`` instance, raises :py:exc:`HPFloatUnavailable`.
71'''
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
91import numpy as num
92from scipy import signal
93import pyrocko
94from pyrocko import dummy_progressbar
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
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
111class URLErrorSSL(URLError):
113 def __init__(self, urlerror):
114 self.__dict__ = urlerror.__dict__.copy()
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))
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
135def urlopen(*args, **kwargs):
137 if 'context' not in kwargs and g_ssl_context is not None:
138 kwargs['context'] = g_ssl_context
140 return _urlopen(*args, **kwargs)
143try:
144 long
145except NameError:
146 long = int
149force_dummy_progressbar = False
152try:
153 from pyrocko import util_ext
154except ImportError:
155 util_ext = None
158logger = logging.getLogger('pyrocko.util')
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
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
178g_setup_logging_args = 'pyrocko', 'warning'
181def setup_logging(programname='pyrocko', levelname='warning'):
182 '''
183 Initialize logging.
185 :param programname: program name to be written in log
186 :param levelname: string indicating the logging level ('debug', 'info',
187 'warning', 'error', 'critical')
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 '''
194 global g_setup_logging_args
195 g_setup_logging_args = (programname, levelname)
197 levels = {'debug': logging.DEBUG,
198 'info': logging.INFO,
199 'warning': logging.WARNING,
200 'error': logging.ERROR,
201 'critical': logging.CRITICAL}
203 logging.basicConfig(
204 level=levels[levelname],
205 format=programname+':%(name)-25s - %(levelname)-8s - %(message)s')
208def subprocess_setup_logging_args():
209 '''
210 Get arguments from previous call to setup_logging.
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
218def data_file(fn):
219 return os.path.join(os.path.split(__file__)[0], 'data', fn)
222class DownloadError(Exception):
223 '''
224 Raised when a download failed.
225 '''
226 pass
229class PathExists(DownloadError):
230 '''
231 Raised when the download target file already exists.
232 '''
233 pass
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):
241 import requests
242 from requests.auth import HTTPBasicAuth
243 from requests.exceptions import HTTPError as req_HTTPError
245 requests.adapters.DEFAULT_RETRIES = 5
246 urljoin = requests.compat.urljoin
248 session = requests.Session()
249 session.header = header
250 session.auth = None if username is None\
251 else HTTPBasicAuth(username, password)
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 }
263 try:
264 url_to_size = {}
266 if callable(status_callback):
267 status_callback(status)
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)
274 if recursive and not url.endswith('/'):
275 url += '/'
277 def parse_directory_tree(url, subdir=''):
278 r = session.get(urljoin(url, subdir))
279 r.raise_for_status()
281 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text)
283 files = sorted(set(subdir + fn for fn in entries
284 if not fn.endswith('/')))
286 if entries_wanted is not None:
287 files = [fn for fn in files
288 if (fn in wanted for wanted in entries_wanted)]
290 status['ntotal_files'] += len(files)
292 dirs = sorted(set(subdir + dn for dn in entries
293 if dn.endswith('/')
294 and dn not in ('./', '../')))
296 for dn in dirs:
297 files.extend(parse_directory_tree(
298 url, subdir=dn))
300 return files
302 def get_content_length(url):
303 if url not in url_to_size:
304 r = session.head(url, headers={'Accept-Encoding': ''})
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)
311 content_length = None
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)
319 url_to_size[url] = content_length
321 return url_to_size[url]
323 def download_file(url, fn):
324 logger.info('starting download of %s...' % url)
325 ensuredirs(fn)
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)
333 r = session.get(url, stream=True, timeout=5)
334 r.raise_for_status()
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)
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)
348 os.rename(fn_tmp, fn)
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))
355 logger.info('finished download of %s' % url)
357 status['nread_files'] += 1
358 if callable(status_callback):
359 status_callback(status)
361 if recursive:
362 if op.exists(fpath) and not force:
363 raise PathExists('path %s already exists' % fpath)
365 files = parse_directory_tree(url)
367 dsize = 0
368 for fn in files:
369 file_url = urljoin(url, fn)
370 dsize += get_content_length(file_url) or 0
372 if method == 'calcsize':
373 return dsize
375 else:
376 for fn in files:
377 file_url = urljoin(url, fn)
378 download_file(file_url, op.join(fpath, fn))
380 else:
381 status['ntotal_files'] += 1
382 if callable(status_callback):
383 status_callback(status)
385 fsize = get_content_length(url)
386 if method == 'calcsize':
387 return fsize
388 else:
389 download_file(url, fpath)
391 except req_HTTPError as e:
392 logging.warning('http error: %s' % e)
393 raise DownloadError('could not download file(s) from: %s' % url)
395 finally:
396 status['finished'] = True
397 if callable(status_callback):
398 status_callback(status)
399 session.close()
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)
412def download_dir(
413 url, fpath, username=None, password=None, status_callback=None,
414 **kwargs):
416 return _download(
417 url, fpath, username, password,
418 recursive=True,
419 status_callback=status_callback,
420 **kwargs)
423class HPFloatUnavailable(Exception):
424 '''
425 Raised when a high precision float type would be required but is not
426 available.
427 '''
428 pass
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.')
438if hasattr(num, 'float128'):
439 hpfloat = num.float128
441elif hasattr(num, 'float96'):
442 hpfloat = num.float96
444else:
445 hpfloat = dummy_hpfloat
448g_time_float = None
449g_time_dtype = None
452class TimeFloatSettingError(Exception):
453 pass
456def use_high_precision_time(enabled):
457 '''
458 Globally force a specific time handling mode.
460 See :ref:`High precision time handling mode <time-handling-mode>`.
462 :param enabled: enable/disable use of high precision time type
463 :type enabled: bool
465 This function should be called before handling/reading any time data.
466 It can only be called once.
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)
475def _setup_high_precision_time_mode(enabled_app=False):
476 global g_time_float
477 global g_time_dtype
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.')
484 from pyrocko import config
486 conf = config.config()
487 enabled_config = conf.use_high_precision_time
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.')
498 enabled = enabled_config
499 mode_from = 'config variable `use_high_precision_time`'
500 notify = enabled
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`'
508 if enabled_app is not None:
509 if enabled_app != enabled:
510 notify = True
511 enabled = enabled_app
512 mode_from = 'application override'
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))
522 if notify:
523 logger.info('Pyrocko high precision time mode %s by %s.' % (
524 'activated' if enabled else 'deactivated',
525 mode_from))
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
535def get_time_float():
536 '''
537 Get the effective float class for timestamps.
539 See :ref:`High precision time handling mode <time-handling-mode>`.
541 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the
542 current time handling mode
543 '''
544 global g_time_float
546 if g_time_float is None:
547 _setup_high_precision_time_mode()
549 return g_time_float
552def get_time_dtype():
553 '''
554 Get effective NumPy float class to handle timestamps.
556 See :ref:`High precision time handling mode <time-handling-mode>`.
557 '''
559 global g_time_dtype
561 if g_time_dtype is None:
562 _setup_high_precision_time_mode()
564 return g_time_dtype
567def to_time_float(t):
568 '''
569 Convert float to valid time stamp in the current time handling mode.
571 See :ref:`High precision time handling mode <time-handling-mode>`.
572 '''
573 return get_time_float()(t)
576class TimestampTypeError(ValueError):
577 pass
580def check_time_class(t, error='raise'):
581 '''
582 Type-check variable against current time handling mode.
584 See :ref:`High precision time handling mode <time-handling-mode>`.
585 '''
587 if t == 0.0:
588 return
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()))
598 if error == 'raise':
599 raise TimestampTypeError(message)
600 elif error == 'warn':
601 logger.warning(message)
602 else:
603 assert False
606class Stopwatch(object):
607 '''
608 Simple stopwatch to measure elapsed wall clock time.
610 Usage::
612 s = Stopwatch()
613 time.sleep(1)
614 print s()
615 time.sleep(1)
616 print s()
617 '''
619 def __init__(self):
620 self.start = time.time()
622 def __call__(self):
623 return time.time() - self.start
626def wrap(text, line_length=80):
627 '''
628 Paragraph and list-aware wrapping of text.
629 '''
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)')
635 paragraphs = at_para.split(text)[::5]
636 listindents = at_para.split(text)[4::5]
637 newlist = at_para.split(text)[3::5]
639 listindents[0:0] = [None]
640 listindents.append(True)
641 newlist.append(None)
643 det_indent = re.compile(r'^ *')
645 outlines = []
646 for ip, p in enumerate(paragraphs):
647 if not p:
648 continue
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)
657 ll = line_length - len(_indent)
658 llf = ll
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)
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):
675 outlines.append('')
677 return outlines
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
692 return '\n'.join(rows)
695class BetterHelpFormatter(optparse.IndentedHelpFormatter):
697 def __init__(self, *args, **kwargs):
698 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs)
700 def format_option(self, option):
701 '''
702 From IndentedHelpFormatter but using a different wrap method.
703 '''
705 help_text_position = 4 + self.current_indent
706 help_text_width = self.width - help_text_position
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)
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('')
728 return '\n'.join(lines)
730 def format_description(self, description):
731 if not description:
732 return ''
734 if self.current_indent == 0:
735 lines = []
736 else:
737 lines = ['']
739 lines.extend(wrap(description, self.width - self.current_indent))
740 if self.current_indent == 0:
741 lines.append('\n')
743 return '\n'.join(
744 ['%*s%s' % (self.current_indent, '', line) for line in lines])
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)
754 def update(self, i):
755 self._task.update(i)
757 def finish(self):
758 self._task.done()
759 if self._context:
760 self._context.__exit__()
761 self._context = None
764def progressbar(label, maxval):
765 if force_dummy_progressbar:
766 return dummy_progressbar.ProgressBar(maxval=maxval).start()
768 return ProgressBar(label, maxval)
771def progress_beg(label):
772 '''
773 Notify user that an operation has started.
775 :param label: name of the operation
777 To be used in conjuction with :py:func:`progress_end`.
778 '''
780 sys.stderr.write(label)
781 sys.stderr.flush()
784def progress_end(label=''):
785 '''
786 Notify user that an operation has ended.
788 :param label: name of the operation
790 To be used in conjuction with :py:func:`progress_beg`.
791 '''
793 sys.stderr.write(' done. %s\n' % label)
794 sys.stderr.flush()
797class ArangeError(ValueError):
798 '''
799 Raised by :py:func:`arange2` for inconsistent range specifications.
800 '''
801 pass
804def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'):
805 '''
806 Return evenly spaced numbers over a specified interval.
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 '''
819 assert error in ('raise', 'round', 'floor', 'ceil')
821 start = dtype(start)
822 stop = dtype(stop)
823 step = dtype(step)
825 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round)
827 n = int(rnd((stop - start) / step)) + 1
828 stop_check = start + (n-1) * step
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))
835 x = num.arange(n, dtype=dtype)
836 x *= step
837 x += start
838 return x
841def polylinefit(x, y, n_or_xnodes):
842 '''
843 Fit piece-wise linear function to data.
845 :param x,y: arrays with coordinates of data
846 :param n_or_xnodes: int, number of segments or x coordinates of polyline
848 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of
849 polyline, root-mean-square error
850 '''
852 x = num.asarray(x)
853 y = num.asarray(y)
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
864 assert len(x) == len(y)
865 assert n > 0
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]
879 a[indices, i*2] = x[indices]
880 a[indices, i*2+1] = 1.0
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
889 d = num.concatenate((y, num.zeros(n-1)))
890 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2))
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
897 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2))
899 return xnodes, ynodes, rms_error
902def plf_integrate_piecewise(x_edges, x, y):
903 '''
904 Calculate definite integral of piece-wise linear function on intervals.
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.
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 '''
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):])
925def diff_fd_1d_4o(dt, data):
926 '''
927 Approximate first derivative of an array (forth order, central FD).
929 :param dt: sampling interval
930 :param data: NumPy array with data samples
932 :returns: NumPy array with same shape as input
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
940 ddata = num.empty_like(data, dtype=float)
942 if data.size >= 5:
943 ddata[2:-2] = scipy.signal.lfilter(
944 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt)
946 if data.size >= 3:
947 ddata[1] = (data[2] - data[0]) / (2. * dt)
948 ddata[-2] = (data[-1] - data[-3]) / (2. * dt)
950 if data.size >= 2:
951 ddata[0] = (data[1] - data[0]) / dt
952 ddata[-1] = (data[-1] - data[-2]) / dt
954 if data.size == 1:
955 ddata[0] = 0.0
957 return ddata
960def diff_fd_1d_2o(dt, data):
961 '''
962 Approximate first derivative of an array (second order, central FD).
964 :param dt: sampling interval
965 :param data: NumPy array with data samples
967 :returns: NumPy array with same shape as input
969 Interior points are approximated to second order, edge points to first
970 order right- or left-sided respectively.
972 Uses :py:func:`numpy.gradient`.
973 '''
975 return num.gradient(data, dt)
978def diff_fd_2d_4o(dt, data):
979 '''
980 Approximate second derivative of an array (forth order, central FD).
982 :param dt: sampling interval
983 :param data: NumPy array with data samples
985 :returns: NumPy array with same shape as input
987 Interior points are approximated to fourth order, next-to-edge points to
988 second order, edge points repeated.
989 '''
990 import scipy.signal
992 ddata = num.empty_like(data, dtype=float)
994 if data.size >= 5:
995 ddata[2:-2] = scipy.signal.lfilter(
996 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2)
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
1002 if data.size < 3:
1003 ddata[:] = 0.0
1005 return ddata
1008def diff_fd_2d_2o(dt, data):
1009 '''
1010 Approximate second derivative of an array (second order, central FD).
1012 :param dt: sampling interval
1013 :param data: NumPy array with data samples
1015 :returns: NumPy array with same shape as input
1017 Interior points are approximated to second order, edge points repeated.
1018 '''
1019 import scipy.signal
1021 ddata = num.empty_like(data, dtype=float)
1023 if data.size >= 3:
1024 ddata[1:-1] = scipy.signal.lfilter(
1025 [1., -2., 1.], [1.], data)[2:] / (dt**2)
1027 ddata[0] = ddata[1]
1028 ddata[-1] = ddata[-2]
1030 if data.size < 3:
1031 ddata[:] = 0.0
1033 return ddata
1036def diff_fd(n, order, dt, data):
1037 '''
1038 Approximate 1st or 2nd derivative of an array.
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
1045 :returns: NumPy array with same shape as input
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`.
1051 Raises :py:exc:`ValueError` for unsupported `n` or `order`.
1052 '''
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}}
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.')
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()))))
1074 return func(dt, data)
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
1087def decimate_coeffs(q, n=None, ftype='iir'):
1089 q = int(q)
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
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')
1104 b = coeffs[n, 1./q]
1105 return b, [1.], n
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
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)
1121 b, a = coeffs[n, 0.05, 0.8/q]
1122 return b, a, n
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
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'.
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`
1138 :returns: the downsampled signal (1D :class:`numpy.ndarray`)
1140 '''
1142 b, a, n = decimate_coeffs(q, n, ftype)
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
1149 y, zf = signal.lfilter(b, a, x, zi=zi_)
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()
1157class UnavailableDecimation(Exception):
1158 '''
1159 Exception raised by :py:func:`decitab` for unavailable decimation factors.
1160 '''
1162 pass
1165def gcd(a, b, epsilon=1e-7):
1166 '''
1167 Greatest common divisor.
1168 '''
1170 while b > epsilon*a:
1171 a, b = b, a % b
1173 return a
1176def lcm(a, b):
1177 '''
1178 Least common multiple.
1179 '''
1181 return a*b // gcd(a, b)
1184def mk_decitab(nmax=100):
1185 '''
1186 Make table with decimation sequences.
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 '''
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
1215 GlobalVars.decitab_nmax = nmax
1218def zfmt(n):
1219 return '%%0%ii' % (int(math.log10(n - 1)) + 1)
1222def _year_to_time(year):
1223 tt = (year, 1, 1, 0, 0, 0)
1224 return to_time_float(calendar.timegm(tt))
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
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
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
1246 except Exception:
1247 return False
1249 return True
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.
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.
1260 :returns: ``(tmin, tmax, year_min, year_max)``
1261 '''
1263 year0 = 2000
1264 year_min = year0
1265 year_max = year0
1267 itests = list(range(101))
1268 for i in range(19):
1269 itests.append(200 + i*100)
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
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
1291 return (
1292 _year_to_time(year_min),
1293 _year_to_time(year_max),
1294 year_min, year_max)
1297g_working_system_time_range = None
1300def get_working_system_time_range():
1301 '''
1302 Caching variant of :py:func:`working_system_time_range`.
1303 '''
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()
1309 return g_working_system_time_range
1312def is_working_time(t):
1313 tmin, tmax, _, _ = get_working_system_time_range()
1314 return tmin <= t <= tmax
1317def julian_day_of_year(timestamp):
1318 '''
1319 Get the day number after the 1st of January of year in ``timestamp``.
1321 :returns: day number as int
1322 '''
1324 return time.gmtime(int(timestamp)).tm_yday
1327def hour_start(timestamp):
1328 '''
1329 Get beginning of hour for any point in time.
1331 :param timestamp: time instant as system timestamp (in seconds)
1333 :returns: instant of hour start as system timestamp
1334 '''
1336 tt = time.gmtime(int(timestamp))
1337 tts = tt[0:4] + (0, 0)
1338 return to_time_float(calendar.timegm(tts))
1341def day_start(timestamp):
1342 '''
1343 Get beginning of day for any point in time.
1345 :param timestamp: time instant as system timestamp (in seconds)
1347 :returns: instant of day start as system timestamp
1348 '''
1350 tt = time.gmtime(int(timestamp))
1351 tts = tt[0:3] + (0, 0, 0)
1352 return to_time_float(calendar.timegm(tts))
1355def month_start(timestamp):
1356 '''
1357 Get beginning of month for any point in time.
1359 :param timestamp: time instant as system timestamp (in seconds)
1361 :returns: instant of month start as system timestamp
1362 '''
1364 tt = time.gmtime(int(timestamp))
1365 tts = tt[0:2] + (1, 0, 0, 0)
1366 return to_time_float(calendar.timegm(tts))
1369def year_start(timestamp):
1370 '''
1371 Get beginning of year for any point in time.
1373 :param timestamp: time instant as system timestamp (in seconds)
1375 :returns: instant of year start as system timestamp
1376 '''
1378 tt = time.gmtime(int(timestamp))
1379 tts = tt[0:1] + (1, 1, 0, 0, 0)
1380 return to_time_float(calendar.timegm(tts))
1383def iter_days(tmin, tmax):
1384 '''
1385 Yields begin and end of days until given time span is covered.
1387 :param tmin,tmax: input time span
1389 :yields: tuples with (begin, end) of days as system timestamps
1390 '''
1392 t = day_start(tmin)
1393 while t < tmax:
1394 tend = day_start(t + 26*60*60)
1395 yield t, tend
1396 t = tend
1399def iter_months(tmin, tmax):
1400 '''
1401 Yields begin and end of months until given time span is covered.
1403 :param tmin,tmax: input time span
1405 :yields: tuples with (begin, end) of months as system timestamps
1406 '''
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
1415def iter_years(tmin, tmax):
1416 '''
1417 Yields begin and end of years until given time span is covered.
1419 :param tmin,tmax: input time span
1421 :yields: tuples with (begin, end) of years as system timestamps
1422 '''
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
1431def today():
1432 return day_start(time.time())
1435def tomorrow():
1436 return day_start(time.time() + 24*60*60)
1439def decitab(n):
1440 '''
1441 Get integer decimation sequence for given downampling factor.
1443 :param n: target decimation factor
1445 :returns: tuple with downsampling sequence
1446 '''
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]
1455def ctimegm(s, format='%Y-%m-%d %H:%M:%S'):
1456 '''
1457 Convert string representing UTC time to system time.
1459 :param s: string to be interpreted
1460 :param format: format string passed to :py:func:`time.strptime`
1462 :returns: system time stamp
1464 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime.
1466 .. note::
1467 This function is to be replaced by :py:func:`str_to_time`.
1468 '''
1470 return calendar.timegm(time.strptime(s, format))
1473def gmctime(t, format='%Y-%m-%d %H:%M:%S'):
1474 '''
1475 Get string representation from system time, UTC.
1477 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime.
1479 .. note::
1480 This function is to be repaced by :py:func:`time_to_str`.
1481 '''
1483 return time.strftime(format, time.gmtime(t))
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.
1491 .. note::
1492 This function is to be replaced by :py:func:`time_to_str`.
1493 '''
1495 return time.strftime(format, time.gmtime(t))
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.
1503 .. note::
1504 This function is to be replaced by :py:func:`time_to_str`.
1505 '''
1507 return time.strftime(format, time.gmtime(t))
1510class TimeStrError(Exception):
1511 '''
1512 Raised for invalid time strings.
1513 '''
1514 pass
1517class FractionalSecondsMissing(TimeStrError):
1518 '''
1519 Exception raised by :py:func:`str_to_time` when the given string lacks
1520 fractional seconds.
1521 '''
1523 pass
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 '''
1532 pass
1535def _endswith_n(s, endings):
1536 for ix, x in enumerate(endings):
1537 if s.endswith(x):
1538 return ix
1539 return -1
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.
1546 :param s: string representing UTC time
1547 :param format: time string format
1548 :returns: system time stamp as floating point value
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 '''
1560 time_float = get_time_float()
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))
1569 return time_float(t)+tfrac
1571 fracsec = 0.
1572 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC'
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))
1581 if iend > 0 and iend != (len(s)-dotpos-1):
1582 raise FractionalSecondsWrongNumberOfDigits(
1583 'string=%s, format=%s' % (s, format))
1585 format = format[:-len(fixed_endings[iend])]
1586 fracsec = float(s[dotpos:])
1587 s = s[:dotpos]
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:])
1595 if dotpos != -1:
1596 s = s[:dotpos]
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))
1605stt = str_to_time
1608def str_to_time_fillup(s):
1609 '''
1610 Default :py:func:`str_to_time` with filling in of missing values.
1612 Allows e.g. `'2010-01-01 00:00:00'` as `'2010-01-01 00:00'`,
1613 `'2010-01-01 00'`, ..., or `'2010'`.
1614 '''
1616 if s == 'now':
1617 return time.time()
1619 if len(s) in (4, 7, 10, 13, 16):
1620 s += '0000-01-01 00:00:00'[len(s):]
1622 return str_to_time(s)
1625def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'):
1626 '''
1627 Get string representation for floating point system time.
1629 :param t: floating point system time
1630 :param format: time string format
1631 :returns: string representing UTC time
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 '''
1639 if pyrocko.grumpy > 0:
1640 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise')
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'
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))
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'])
1661 ts = float(num.floor(t))
1662 tfrac = t-ts
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.
1670 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1)
1672 return time.strftime(format, time.gmtime(ts))
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()
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.
1684 if fmt is None:
1685 fmt = '%Y-%m-%d %H:%M:%S .%r'
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])
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)
1698def gmtime_x(timestamp):
1699 etimestamp = float(num.floor(timestamp))
1700 tt = time.gmtime(etimestamp)
1701 ms = (timestamp-etimestamp)*1000
1702 return tt, ms
1705def plural_s(n):
1706 if not isinstance(n, int):
1707 n = len(n)
1709 return 's' if n != 1 else ''
1712def ensuredirs(dst):
1713 '''
1714 Create all intermediate path components for a target path.
1716 :param dst: target path
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 '''
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)
1728 dirs.reverse()
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
1738def ensuredir(dst):
1739 '''
1740 Create directory and all intermediate path components to it as needed.
1742 :param dst: directory name
1744 Nothing is done if the given target already exists.
1745 '''
1747 if os.path.exists(dst):
1748 return
1750 dst.rstrip(os.sep)
1752 ensuredirs(dst)
1753 try:
1754 os.mkdir(dst)
1755 except OSError as e:
1756 if not e.errno == errno.EEXIST:
1757 raise
1760def reuse(x):
1761 '''
1762 Get unique instance of an object.
1764 :param x: hashable object
1765 :returns: reference to x or an equivalent object
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 '''
1771 grs = GlobalVars.reuse_store
1772 if x not in grs:
1773 grs[x] = x
1774 return grs[x]
1777def deuse(x):
1778 grs = GlobalVars.reuse_store
1779 if x in grs:
1780 del grs[x]
1783class Anon(object):
1784 '''
1785 Dict-to-object utility.
1787 Any given arguments are stored as attributes.
1789 Example::
1791 a = Anon(x=1, y=2)
1792 print a.x, a.y
1793 '''
1795 def __init__(self, **dict):
1796 for k in dict:
1797 self.__dict__[k] = dict[k]
1800def iter_select_files(
1801 paths,
1802 include=None,
1803 exclude=None,
1804 selector=None,
1805 show_progress=True,
1806 pass_through=None):
1808 '''
1809 Recursively select files (generator variant).
1811 See :py:func:`select_files`.
1812 '''
1814 if show_progress:
1815 progress_beg('selecting files...')
1817 ngood = 0
1818 check_include = None
1819 if include is not None:
1820 rinclude = re.compile(include)
1822 def check_include(path):
1823 m = rinclude.search(path)
1824 if not m:
1825 return False
1827 if selector is None:
1828 return True
1830 infos = Anon(**m.groupdict())
1831 return selector(infos)
1833 check_exclude = None
1834 if exclude is not None:
1835 rexclude = re.compile(exclude)
1837 def check_exclude(path):
1838 return not bool(rexclude.search(path))
1840 if check_include and check_exclude:
1842 def check(path):
1843 return check_include(path) and check_exclude(path)
1845 elif check_include:
1846 check = check_include
1848 elif check_exclude:
1849 check = check_exclude
1851 else:
1852 check = None
1854 if isinstance(paths, str):
1855 paths = [paths]
1857 for path in paths:
1858 if pass_through and pass_through(path):
1859 if check is None or check(path):
1860 yield path
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
1876 if show_progress:
1877 progress_end('%i file%s selected.' % (ngood, plural_s(ngood)))
1880def select_files(
1881 paths, include=None, exclude=None, selector=None, show_progress=True,
1882 regex=None):
1884 '''
1885 Recursively select files.
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
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``.
1903 Examples
1905 To find all files ending in ``'.mseed'`` or ``'.msd'``::
1907 select_files(paths,
1908 include=r'\\.(mseed|msd)$')
1910 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for
1911 the year::
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
1921 return list(iter_select_files(
1922 paths, include, exclude, selector, show_progress))
1925def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
1926 '''
1927 Convert positive integer to a base36 string.
1928 '''
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')
1935 # Special case for small numbers
1936 if number < 36:
1937 return alphabet[number]
1939 base36 = ''
1940 while number != 0:
1941 number, i = divmod(number, 36)
1942 base36 = alphabet[i] + base36
1944 return base36
1947def base36decode(number):
1948 '''
1949 Decode base36 endcoded positive integer.
1950 '''
1952 return int(number, 36)
1955class UnpackError(Exception):
1956 '''
1957 Exception raised when :py:func:`unpack_fixed` encounters an error.
1958 '''
1960 pass
1963ruler = ''.join(['%-10i' % i for i in range(8)]) \
1964 + '\n' + '0123456789' * 8 + '\n'
1967def unpack_fixed(format, line, *callargs):
1968 '''
1969 Unpack fixed format string, as produced by many fortran codes.
1971 :param format: format specification
1972 :param line: string to be processed
1973 :param callargs: callbacks for callback fields in the format
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).
1981 The following field types are available:
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 '''
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]
2013 if cast == 'extra':
2014 cast = callargs[icall]
2015 icall += 1
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))
2032 ipos += ln
2034 return values
2037_pattern_cache = {}
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]
2047 return rpattern
2050def match_nslc(patterns, nslc):
2051 '''
2052 Match network-station-location-channel code against pattern or list of
2053 patterns.
2055 :param patterns: pattern or list of patterns
2056 :param nslc: tuple with (network, station, location, channel) as strings
2058 :returns: ``True`` if the pattern matches or if any of the given patterns
2059 match; or ``False``.
2061 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq].
2063 Example::
2065 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True
2066 '''
2068 if isinstance(patterns, str):
2069 patterns = [patterns]
2071 if not isinstance(nslc, str):
2072 s = '.'.join(nslc)
2073 else:
2074 s = nslc
2076 for pattern in patterns:
2077 if _nslc_pattern(pattern).match(s):
2078 return True
2080 return False
2083def match_nslcs(patterns, nslcs):
2084 '''
2085 Get network-station-location-channel codes that match given pattern or any
2086 of several given patterns.
2088 :param patterns: pattern or list of patterns
2089 :param nslcs: list of (network, station, location, channel) tuples
2091 See also :py:func:`match_nslc`
2092 '''
2094 matching = []
2095 for nslc in nslcs:
2096 if match_nslc(patterns, nslc):
2097 matching.append(nslc)
2099 return matching
2102class Timeout(Exception):
2103 pass
2106def create_lockfile(fn, timeout=None, timewarn=10.):
2107 t0 = time.time()
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
2115 except OSError as e:
2116 if e.errno in (errno.EEXIST, 13): # 13 occurs on windows
2117 pass # retry
2118 else:
2119 raise
2121 tnow = time.time()
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))
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))
2133 timewarn *= 2
2135 time.sleep(0.01)
2138def delete_lockfile(fn):
2139 os.unlink(fn)
2142class Lockfile(Exception):
2144 def __init__(self, path, timeout=5, timewarn=10.):
2145 self._path = path
2146 self._timeout = timeout
2147 self._timewarn = timewarn
2149 def __enter__(self):
2150 create_lockfile(
2151 self._path, timeout=self._timeout, timewarn=self._timewarn)
2152 return None
2154 def __exit__(self, type, value, traceback):
2155 delete_lockfile(self._path)
2158class SoleError(Exception):
2159 '''
2160 Exception raised by objects of type :py:class:`Sole`, when an concurrent
2161 instance is running.
2162 '''
2164 pass
2167class Sole(object):
2168 '''
2169 Use POSIX advisory file locking to ensure that only a single instance of a
2170 program is running.
2172 :param pid_path: path to lockfile to be used
2174 Usage::
2176 from pyrocko.util import Sole, SoleError, setup_logging
2177 import os
2179 setup_logging('my_program')
2181 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock')
2182 try:
2183 sole = Sole(pid_path)
2185 except SoleError, e:
2186 logger.fatal( str(e) )
2187 sys.exit(1)
2188 '''
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
2197 if not fcntl:
2198 raise SoleError(
2199 'Python standard library module "fcntl" not available on '
2200 'this platform.')
2202 self._fcntl = fcntl
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)
2210 try:
2211 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB)
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 = '?'
2222 raise SoleError('Other instance is running (pid = %s)' % pid)
2224 try:
2225 os.ftruncate(self._lockfile, 0)
2226 os.write(self._lockfile, '%i\n' % os.getpid())
2227 os.fsync(self._lockfile)
2229 except Exception:
2230 # the pid is only stored for user information, so this is allowed
2231 # to fail
2232 pass
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
2245re_escapequotes = re.compile(r"(['\\])")
2248def escapequotes(s):
2249 return re_escapequotes.sub(r'\\\1', s)
2252class TableWriter(object):
2253 '''
2254 Write table of space separated values to a file.
2256 :param f: file like object
2258 Strings containing spaces are quoted on output.
2259 '''
2261 def __init__(self, f):
2262 self._f = f
2264 def writerow(self, row, minfieldwidths=None):
2266 '''
2267 Write one row of values to underlying file.
2269 :param row: iterable of values
2270 :param minfieldwidths: minimum field widths for the values
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 '''
2279 out = []
2281 for i, x in enumerate(row):
2282 w = 0
2283 if minfieldwidths and i < len(minfieldwidths):
2284 w = minfieldwidths[i]
2286 if isinstance(x, str):
2287 if re.search(r"\s|'", x):
2288 x = "'%s'" % escapequotes(x)
2290 x = x.ljust(w)
2291 else:
2292 x = str(x).rjust(w)
2294 out.append(x)
2296 self._f.write(' '.join(out).rstrip() + '\n')
2299class TableReader(object):
2301 '''
2302 Read table of space separated values from a file.
2304 :param f: file-like object
2306 This uses Pythons shlex module to tokenize lines. Should deal correctly
2307 with quoted strings.
2308 '''
2310 def __init__(self, f):
2311 self._f = f
2312 self.eof = False
2314 def readrow(self):
2315 '''
2316 Read one row from the underlying file, tokenize it with shlex.
2318 :returns: tokenized line as a list of strings.
2319 '''
2321 line = self._f.readline()
2322 if not line:
2323 self.eof = True
2324 return []
2325 line.strip()
2326 if line.startswith('#'):
2327 return []
2329 return qsplit(line)
2332def gform(number, significant_digits=3):
2333 '''
2334 Pretty print floating point numbers.
2336 Align floating point numbers at the decimal dot.
2338 ::
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|
2353 The formatted string has length ``significant_digits * 2 + 6``.
2354 '''
2356 no_exp_range = (pow(10., -1),
2357 pow(10., significant_digits))
2358 width = significant_digits+significant_digits-1+1+1+5
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
2370def human_bytesize(value):
2372 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split()
2374 if value == 1:
2375 return '1 Byte'
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)
2384 return '%i Bytes' % value
2387re_compatibility = re.compile(
2388 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' +
2389 r'(dummy|poel|qseis|qssp))\.'
2390)
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
2400 return oldstyle
2403def pf_upgrade(fn):
2404 need = pf_is_old(fn)
2405 if need:
2406 fn_temp = fn + '.temp'
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)
2414 os.rename(fn_temp, fn)
2416 return need
2419def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'):
2420 '''
2421 Extract leap second information from tzdata.
2423 Based on example at http://stackoverflow.com/questions/19332902/\
2424 extract-historic-leap-seconds-from-tzdata
2426 See also 'man 5 tzfile'.
2427 '''
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'
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))
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))
2449 return out
2452class LeapSecondsError(Exception):
2453 pass
2456class LeapSecondsOutdated(LeapSecondsError):
2457 pass
2460class InvalidLeapSecondsFile(LeapSecondsOutdated):
2461 pass
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
2472 tnow = int(round(time.time()))
2474 if not op.exists(fn):
2475 raise LeapSecondsOutdated('no leap seconds file found')
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')
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))
2493 except IOError:
2494 raise LeapSecondsError('cannot read leap seconds file %s' % fn)
2496 if texpires is None or tnow > texpires:
2497 raise LeapSecondsOutdated('leap seconds list is outdated')
2499 return data
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)
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)
2518 for i in range(3):
2519 try:
2520 return parse_leap_seconds_list(fn)
2522 except LeapSecondsOutdated:
2523 try:
2524 logger.info('updating leap seconds list...')
2525 download_file(url, fn)
2527 except Exception as e:
2528 raise LeapSecondsError(
2529 'cannot download leap seconds list from %s to %s (%s)'
2530 % (url, fn, e))
2532 raise LeapSecondsError('Could not retrieve/read leap seconds file.')
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
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
2549 return ls[-1][1] - 9
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()
2558 if t_gps < ls[0][0] + ls[0][1] - 9:
2559 return - (ls[0][1] - 1 - 9)
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
2568 return - (ls[-1][1] - 9)
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
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
2582 except FileLoadError as e:
2583 e.set_context('filename', filename)
2584 raise
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
2593 def iload_glob(pattern, **kwargs):
2595 for fn in glob.iglob(pattern):
2596 for cr in iload_filename(fn, **kwargs):
2597 yield cr
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)
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)
2614 iload_filename.__doc__ = '''
2615 Read %s information from named file.
2616 ''' % doc_fmt
2618 iload_dirname.__doc__ = '''
2619 Read %s information from directory of %s files.
2620 ''' % (doc_fmt, doc_fmt)
2622 iload_glob.__doc__ = '''
2623 Read %s information from files matching a glob pattern.
2624 ''' % doc_fmt
2626 iload.__doc__ = '''
2627 Load %s information from given source(s)
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.
2633 This function behaves as a generator yielding %s objects.
2634 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects)
2636 for f in iload_filename, iload_dirname, iload_glob, iload:
2637 f.__module__ = iload_fh.__module__
2639 return iload_filename, iload_dirname, iload_glob, iload
2642class Inconsistency(Exception):
2643 pass
2646def consistency_check(list_of_tuples, message='values differ:'):
2647 '''
2648 Check for inconsistencies.
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 '''
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))
2662class defaultzerodict(dict):
2663 def __missing__(self, k):
2664 return 0
2667def mostfrequent(x):
2668 c = defaultzerodict()
2669 for e in x:
2670 c[e] += 1
2672 return sorted(list(c.keys()), key=lambda k: c[k])[-1]
2675def consistency_merge(list_of_tuples,
2676 message='values differ:',
2677 error='raise',
2678 merge=mostfrequent):
2680 assert error in ('raise', 'warn', 'ignore')
2682 if len(list_of_tuples) == 0:
2683 raise Exception('cannot merge empty sequence')
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
2692 elif error == 'warn':
2693 logger.warning(str(e))
2695 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]])
2698def short_to_list(nmax, it):
2699 import itertools
2701 if isinstance(it, list):
2702 return it
2704 li = []
2705 for i in range(nmax+1):
2706 try:
2707 li.append(next(it))
2708 except StopIteration:
2709 return li
2711 return itertools.chain(li, it)
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
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
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()
2740g_re_qsplit = re.compile(
2741 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)')
2742g_re_qsplit_sep = {}
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]
2759g_re_trivial = re.compile(r'\A[^\'"\s]+\Z')
2760g_re_trivial_sep = {}
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')
2773 return g_re_trivial_sep[sep]
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'\\([\\"])')
2782def escape_s(s):
2783 '''
2784 Backslash-escape single-quotes and backslashes.
2786 Example: ``Jack's`` => ``Jack\\'s``
2788 '''
2789 return g_re_escape_s.sub(r'\\\1', s)
2792def unescape_s(s):
2793 '''
2794 Unescape backslash-escaped single-quotes and backslashes.
2796 Example: ``Jack\\'s`` => ``Jack's``
2797 '''
2798 return g_re_unescape_s.sub(r'\1', s)
2801def escape_d(s):
2802 '''
2803 Backslash-escape double-quotes and backslashes.
2805 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"``
2806 '''
2807 return g_re_escape_d.sub(r'\\\1', s)
2810def unescape_d(s):
2811 '''
2812 Unescape backslash-escaped double-quotes and backslashes.
2814 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"``
2815 '''
2816 return g_re_unescape_d.sub(r'\1', s)
2819def qjoin_s(it, sep=None):
2820 '''
2821 Join sequence of strings into a line, single-quoting non-trivial strings.
2823 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"``
2824 '''
2825 re_trivial = get_re_trivial(sep)
2827 if sep is None:
2828 sep = ' '
2830 return sep.join(
2831 w if re_trivial.search(w) else "'%s'" % escape_s(w) for w in it)
2834def qjoin_d(it, sep=None):
2835 '''
2836 Join sequence of strings into a line, double-quoting non-trivial strings.
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 = ' '
2845 return sep.join(
2846 w if re_trivial.search(w) else '"%s"' % escape_d(w) for w in it)
2849def qsplit(s, sep=None):
2850 '''
2851 Split line into list of strings, allowing for quoted strings.
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)]
2864g_have_warned_threadpoolctl = False
2867class threadpool_limits_dummy(object):
2869 def __init__(self, *args, **kwargs):
2870 pass
2872 def __enter__(self):
2873 global g_have_warned_threadpoolctl
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`.')
2881 g_have_warned_threadpoolctl = True
2883 return self
2885 def __exit__(self, type, value, traceback):
2886 pass
2889def get_threadpool_limits():
2890 '''
2891 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail.
2892 '''
2894 try:
2895 from threadpoolctl import threadpool_limits
2896 return threadpool_limits
2898 except ImportError:
2899 return threadpool_limits_dummy
2902def fmt_summary(entries, widths):
2903 return ' | '.join(
2904 entry.ljust(width) for (entry, width) in zip(entries, widths))
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)