Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/util.py: 78%
1238 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-03-07 11:54 +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')
161# fallbacks num_full and num_full_like are not needed anymore but
162# kept here because downstream code may still use these.
163try:
164 num_full = num.full
165except AttributeError:
166 def num_full(shape, fill_value, dtype=None, order='C'):
167 a = num.empty(shape, dtype=dtype, order=order)
168 a.fill(fill_value)
169 return a
171try:
172 num_full_like = num.full_like
173except AttributeError:
174 def num_full_like(arr, fill_value, dtype=None, order='K', subok=True):
175 a = num.empty_like(arr, dtype=dtype, order=order, subok=subok)
176 a.fill(fill_value)
177 return a
180g_setup_logging_args = 'pyrocko', 'warning'
183def setup_logging(programname='pyrocko', levelname='warning'):
184 '''
185 Initialize logging.
187 :param programname: program name to be written in log
188 :param levelname: string indicating the logging level ('debug', 'info',
189 'warning', 'error', 'critical')
191 This function is called at startup by most pyrocko programs to set up a
192 consistent logging format. This is simply a shortcut to a call to
193 :py:func:`logging.basicConfig()`.
194 '''
196 global g_setup_logging_args
197 g_setup_logging_args = (programname, levelname)
199 levels = {'debug': logging.DEBUG,
200 'info': logging.INFO,
201 'warning': logging.WARNING,
202 'error': logging.ERROR,
203 'critical': logging.CRITICAL}
205 logging.basicConfig(
206 level=levels[levelname],
207 format=programname+':%(name)-25s - %(levelname)-8s - %(message)s')
210def subprocess_setup_logging_args():
211 '''
212 Get arguments from previous call to setup_logging.
214 These can be sent down to a worker process so it can setup its logging
215 in the same way as the main process.
216 '''
217 return g_setup_logging_args
220def data_file(fn):
221 return os.path.join(os.path.split(__file__)[0], 'data', fn)
224class DownloadError(Exception):
225 '''
226 Raised when a download failed.
227 '''
228 pass
231class PathExists(DownloadError):
232 '''
233 Raised when the download target file already exists.
234 '''
235 pass
238def _download(url, fpath, username=None, password=None,
239 force=False, method='download', stats=None,
240 status_callback=None, entries_wanted=None,
241 recursive=False, header=None):
243 import requests
244 from requests.auth import HTTPBasicAuth
245 from requests.exceptions import HTTPError as req_HTTPError
247 requests.adapters.DEFAULT_RETRIES = 5
248 urljoin = requests.compat.urljoin
250 session = requests.Session()
251 session.header = header
252 session.auth = None if username is None\
253 else HTTPBasicAuth(username, password)
255 status = {
256 'ntotal_files': 0,
257 'nread_files': 0,
258 'ntotal_bytes_all_files': 0,
259 'nread_bytes_all_files': 0,
260 'ntotal_bytes_current_file': 0,
261 'nread_bytes_current_file': 0,
262 'finished': False
263 }
265 try:
266 url_to_size = {}
268 if callable(status_callback):
269 status_callback(status)
271 if not recursive and url.endswith('/'):
272 raise DownloadError(
273 'URL: %s appears to be a directory'
274 ' but recurvise download is False' % url)
276 if recursive and not url.endswith('/'):
277 url += '/'
279 def parse_directory_tree(url, subdir=''):
280 r = session.get(urljoin(url, subdir))
281 r.raise_for_status()
283 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text)
285 files = sorted(set(subdir + fn for fn in entries
286 if not fn.endswith('/')))
288 if entries_wanted is not None:
289 files = [fn for fn in files
290 if (fn in wanted for wanted in entries_wanted)]
292 status['ntotal_files'] += len(files)
294 dirs = sorted(set(subdir + dn for dn in entries
295 if dn.endswith('/')
296 and dn not in ('./', '../')))
298 for dn in dirs:
299 files.extend(parse_directory_tree(
300 url, subdir=dn))
302 return files
304 def get_content_length(url):
305 if url not in url_to_size:
306 r = session.head(url, headers={'Accept-Encoding': ''})
308 content_length = r.headers.get('content-length', None)
309 if content_length is None:
310 logger.debug('Could not get HTTP header '
311 'Content-Length for %s' % url)
313 content_length = None
315 else:
316 content_length = int(content_length)
317 status['ntotal_bytes_all_files'] += content_length
318 if callable(status_callback):
319 status_callback(status)
321 url_to_size[url] = content_length
323 return url_to_size[url]
325 def download_file(url, fn):
326 logger.info('starting download of %s...' % url)
327 ensuredirs(fn)
329 fsize = get_content_length(url)
330 status['ntotal_bytes_current_file'] = fsize
331 status['nread_bytes_current_file'] = 0
332 if callable(status_callback):
333 status_callback(status)
335 r = session.get(url, stream=True, timeout=5)
336 r.raise_for_status()
338 frx = 0
339 fn_tmp = fn + '.%i.temp' % os.getpid()
340 with open(fn_tmp, 'wb') as f:
341 for d in r.iter_content(chunk_size=1024):
342 f.write(d)
343 frx += len(d)
345 status['nread_bytes_all_files'] += len(d)
346 status['nread_bytes_current_file'] += len(d)
347 if callable(status_callback):
348 status_callback(status)
350 os.rename(fn_tmp, fn)
352 if fsize is not None and frx != fsize:
353 logger.warning(
354 'HTTP header Content-Length: %i bytes does not match '
355 'download size %i bytes' % (fsize, frx))
357 logger.info('finished download of %s' % url)
359 status['nread_files'] += 1
360 if callable(status_callback):
361 status_callback(status)
363 if recursive:
364 if op.exists(fpath) and not force:
365 raise PathExists('path %s already exists' % fpath)
367 files = parse_directory_tree(url)
369 dsize = 0
370 for fn in files:
371 file_url = urljoin(url, fn)
372 dsize += get_content_length(file_url) or 0
374 if method == 'calcsize':
375 return dsize
377 else:
378 for fn in files:
379 file_url = urljoin(url, fn)
380 download_file(file_url, op.join(fpath, fn))
382 else:
383 status['ntotal_files'] += 1
384 if callable(status_callback):
385 status_callback(status)
387 fsize = get_content_length(url)
388 if method == 'calcsize':
389 return fsize
390 else:
391 download_file(url, fpath)
393 except req_HTTPError as e:
394 logging.warning('http error: %s' % e)
395 raise DownloadError('could not download file(s) from: %s' % url)
397 finally:
398 status['finished'] = True
399 if callable(status_callback):
400 status_callback(status)
401 session.close()
404def download_file(
405 url, fpath, username=None, password=None, status_callback=None,
406 **kwargs):
407 return _download(
408 url, fpath, username, password,
409 recursive=False,
410 status_callback=status_callback,
411 **kwargs)
414def download_dir(
415 url, fpath, username=None, password=None, status_callback=None,
416 **kwargs):
418 return _download(
419 url, fpath, username, password,
420 recursive=True,
421 status_callback=status_callback,
422 **kwargs)
425class HPFloatUnavailable(Exception):
426 '''
427 Raised when a high precision float type would be required but is not
428 available.
429 '''
430 pass
433class dummy_hpfloat(object):
434 def __init__(self, *args, **kwargs):
435 raise HPFloatUnavailable(
436 'NumPy lacks support for float128 or float96 data type on this '
437 'platform.')
440if hasattr(num, 'float128'):
441 hpfloat = num.float128
443elif hasattr(num, 'float96'):
444 hpfloat = num.float96
446else:
447 hpfloat = dummy_hpfloat
450g_time_float = None
451g_time_dtype = None
454class TimeFloatSettingError(Exception):
455 pass
458def use_high_precision_time(enabled):
459 '''
460 Globally force a specific time handling mode.
462 See :ref:`High precision time handling mode <time-handling-mode>`.
464 :param enabled: enable/disable use of high precision time type
465 :type enabled: bool
467 This function should be called before handling/reading any time data.
468 It can only be called once.
470 Special attention is required when using multiprocessing on a platform
471 which does not use fork under the hood. In such cases, the desired setting
472 must be set also in the subprocess.
473 '''
474 _setup_high_precision_time_mode(enabled_app=enabled)
477def _setup_high_precision_time_mode(enabled_app=False):
478 global g_time_float
479 global g_time_dtype
481 if not (g_time_float is None and g_time_dtype is None):
482 raise TimeFloatSettingError(
483 'Cannot set time handling mode: too late, it has already been '
484 'fixed by an earlier call.')
486 from pyrocko import config
488 conf = config.config()
489 enabled_config = conf.use_high_precision_time
491 enabled_env = os.environ.get('PYROCKO_USE_HIGH_PRECISION_TIME', None)
492 if enabled_env is not None:
493 try:
494 enabled_env = int(enabled_env) == 1
495 except ValueError:
496 raise TimeFloatSettingError(
497 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME '
498 'should be set to 0 or 1.')
500 enabled = enabled_config
501 mode_from = 'config variable `use_high_precision_time`'
502 notify = enabled
504 if enabled_env is not None:
505 if enabled_env != enabled:
506 notify = True
507 enabled = enabled_env
508 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`'
510 if enabled_app is not None:
511 if enabled_app != enabled:
512 notify = True
513 enabled = enabled_app
514 mode_from = 'application override'
516 logger.debug('''
517Pyrocko high precision time mode selection (latter override earlier):
518 config: %s
519 env: %s
520 app: %s
521 -> enabled: %s'''.lstrip() % (
522 enabled_config, enabled_env, enabled_app, enabled))
524 if notify:
525 logger.info('Pyrocko high precision time mode %s by %s.' % (
526 'activated' if enabled else 'deactivated',
527 mode_from))
529 if enabled:
530 g_time_float = hpfloat
531 g_time_dtype = hpfloat
532 else:
533 g_time_float = float
534 g_time_dtype = num.float64
537def get_time_float():
538 '''
539 Get the effective float class for timestamps.
541 See :ref:`High precision time handling mode <time-handling-mode>`.
543 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the
544 current time handling mode
545 '''
546 global g_time_float
548 if g_time_float is None:
549 _setup_high_precision_time_mode()
551 return g_time_float
554def get_time_dtype():
555 '''
556 Get effective NumPy float class to handle timestamps.
558 See :ref:`High precision time handling mode <time-handling-mode>`.
559 '''
561 global g_time_dtype
563 if g_time_dtype is None:
564 _setup_high_precision_time_mode()
566 return g_time_dtype
569def to_time_float(t):
570 '''
571 Convert float to valid time stamp in the current time handling mode.
573 See :ref:`High precision time handling mode <time-handling-mode>`.
574 '''
575 return get_time_float()(t)
578class TimestampTypeError(ValueError):
579 pass
582def check_time_class(t, error='raise'):
583 '''
584 Type-check variable against current time handling mode.
586 See :ref:`High precision time handling mode <time-handling-mode>`.
587 '''
589 if t == 0.0:
590 return
592 if not isinstance(t, get_time_float()):
593 message = (
594 'Timestamp %g is of type %s but should be of type %s with '
595 "Pyrocko's currently selected time handling mode.\n\n"
596 'See https://pyrocko.org/docs/current/library/reference/util.html'
597 '#high-precision-time-handling-mode' % (
598 t, type(t), get_time_float()))
600 if error == 'raise':
601 raise TimestampTypeError(message)
602 elif error == 'warn':
603 logger.warning(message)
604 else:
605 assert False
608class Stopwatch(object):
609 '''
610 Simple stopwatch to measure elapsed wall clock time.
612 Usage::
614 s = Stopwatch()
615 time.sleep(1)
616 print s()
617 time.sleep(1)
618 print s()
619 '''
621 def __init__(self):
622 self.start = time.time()
624 def __call__(self):
625 return time.time() - self.start
628def wrap(text, line_length=80):
629 '''
630 Paragraph and list-aware wrapping of text.
631 '''
633 text = text.strip('\n')
634 at_lineend = re.compile(r' *\n')
635 at_para = re.compile(r'((^|(\n\s*)?\n)(\s+[*] )|\n\s*\n)')
637 paragraphs = at_para.split(text)[::5]
638 listindents = at_para.split(text)[4::5]
639 newlist = at_para.split(text)[3::5]
641 listindents[0:0] = [None]
642 listindents.append(True)
643 newlist.append(None)
645 det_indent = re.compile(r'^ *')
647 outlines = []
648 for ip, p in enumerate(paragraphs):
649 if not p:
650 continue
652 if listindents[ip] is None:
653 _indent = det_indent.findall(p)[0]
654 findent = _indent
655 else:
656 findent = listindents[ip]
657 _indent = ' ' * len(findent)
659 ll = line_length - len(_indent)
660 llf = ll
662 oldlines = [s.strip() for s in at_lineend.split(p.rstrip())]
663 p1 = ' '.join(oldlines)
664 possible = re.compile(r'(^.{1,%i}|.{1,%i})( |$)' % (llf, ll))
665 for imatch, match in enumerate(possible.finditer(p1)):
666 parout = match.group(1)
667 if imatch == 0:
668 outlines.append(findent + parout)
669 else:
670 outlines.append(_indent + parout)
672 if ip != len(paragraphs)-1 and (
673 listindents[ip] is None
674 or newlist[ip] is not None
675 or listindents[ip+1] is None):
677 outlines.append('')
679 return outlines
682def ewrap(lines, width=80, indent=''):
683 lines = list(lines)
684 if not lines:
685 return ''
686 fwidth = max(len(s) for s in lines)
687 nx = max(1, (80-len(indent)) // (fwidth+1))
688 i = 0
689 rows = []
690 while i < len(lines):
691 rows.append(indent + ' '.join(x.ljust(fwidth) for x in lines[i:i+nx]))
692 i += nx
694 return '\n'.join(rows)
697class BetterHelpFormatter(optparse.IndentedHelpFormatter):
699 def __init__(self, *args, **kwargs):
700 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs)
702 def format_option(self, option):
703 '''
704 From IndentedHelpFormatter but using a different wrap method.
705 '''
707 help_text_position = 4 + self.current_indent
708 help_text_width = self.width - help_text_position
710 opts = self.option_strings[option]
711 opts = '%*s%s' % (self.current_indent, '', opts)
712 if option.help:
713 help_text = self.expand_default(option)
715 if self.help_position + len(help_text) + 1 <= self.width:
716 lines = [
717 '',
718 '%-*s %s' % (self.help_position, opts, help_text),
719 '']
720 else:
721 lines = ['']
722 lines.append(opts)
723 lines.append('')
724 if option.help:
725 help_lines = wrap(help_text, help_text_width)
726 lines.extend(['%*s%s' % (help_text_position, '', line)
727 for line in help_lines])
728 lines.append('')
730 return '\n'.join(lines)
732 def format_description(self, description):
733 if not description:
734 return ''
736 if self.current_indent == 0:
737 lines = []
738 else:
739 lines = ['']
741 lines.extend(wrap(description, self.width - self.current_indent))
742 if self.current_indent == 0:
743 lines.append('\n')
745 return '\n'.join(
746 ['%*s%s' % (self.current_indent, '', line) for line in lines])
749class ProgressBar:
750 def __init__(self, label, n):
751 from pyrocko import progress
752 self._context = progress.view()
753 self._context.__enter__()
754 self._task = progress.task(label, n)
756 def update(self, i):
757 self._task.update(i)
759 def finish(self):
760 self._task.done()
761 if self._context:
762 self._context.__exit__()
763 self._context = None
766def progressbar(label, maxval):
767 if force_dummy_progressbar:
768 return dummy_progressbar.ProgressBar(maxval=maxval).start()
770 return ProgressBar(label, maxval)
773def progress_beg(label):
774 '''
775 Notify user that an operation has started.
777 :param label: name of the operation
779 To be used in conjuction with :py:func:`progress_end`.
780 '''
782 sys.stderr.write(label)
783 sys.stderr.flush()
786def progress_end(label=''):
787 '''
788 Notify user that an operation has ended.
790 :param label: name of the operation
792 To be used in conjuction with :py:func:`progress_beg`.
793 '''
795 sys.stderr.write(' done. %s\n' % label)
796 sys.stderr.flush()
799class ArangeError(ValueError):
800 '''
801 Raised by :py:func:`arange2` for inconsistent range specifications.
802 '''
803 pass
806def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'):
807 '''
808 Return evenly spaced numbers over a specified interval.
810 Like :py:func:`numpy.arange` but returning floating point numbers by
811 default and with defined behaviour when stepsize is inconsistent with
812 interval bounds. It is considered inconsistent if the difference between
813 the closest multiple of ``step`` and ``stop`` is larger than ``epsilon *
814 step``. Inconsistencies are handled according to the ``error`` parameter.
815 If it is set to ``'raise'`` an exception of type :py:exc:`ArangeError` is
816 raised. If it is set to ``'round'``, ``'floor'``, or ``'ceil'``, ``stop``
817 is silently changed to the closest, the next smaller, or next larger
818 multiple of ``step``, respectively.
819 '''
821 assert error in ('raise', 'round', 'floor', 'ceil')
823 start = dtype(start)
824 stop = dtype(stop)
825 step = dtype(step)
827 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round)
829 n = int(rnd((stop - start) / step)) + 1
830 stop_check = start + (n-1) * step
832 if error == 'raise' and abs(stop_check - stop) > step * epsilon:
833 raise ArangeError(
834 'inconsistent range specification: start=%g, stop=%g, step=%g'
835 % (start, stop, step))
837 x = num.arange(n, dtype=dtype)
838 x *= step
839 x += start
840 return x
843def polylinefit(x, y, n_or_xnodes):
844 '''
845 Fit piece-wise linear function to data.
847 :param x,y: arrays with coordinates of data
848 :param n_or_xnodes: int, number of segments or x coordinates of polyline
850 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of
851 polyline, root-mean-square error
852 '''
854 x = num.asarray(x)
855 y = num.asarray(y)
857 if isinstance(n_or_xnodes, int):
858 n = n_or_xnodes
859 xmin = x.min()
860 xmax = x.max()
861 xnodes = num.linspace(xmin, xmax, n+1)
862 else:
863 xnodes = num.asarray(n_or_xnodes)
864 n = xnodes.size - 1
866 assert len(x) == len(y)
867 assert n > 0
869 ndata = len(x)
870 a = num.zeros((ndata+(n-1), n*2))
871 for i in range(n):
872 xmin_block = xnodes[i]
873 xmax_block = xnodes[i+1]
874 if i == n-1: # don't loose last point
875 indices = num.where(
876 num.logical_and(xmin_block <= x, x <= xmax_block))[0]
877 else:
878 indices = num.where(
879 num.logical_and(xmin_block <= x, x < xmax_block))[0]
881 a[indices, i*2] = x[indices]
882 a[indices, i*2+1] = 1.0
884 w = float(ndata)*100.
885 if i < n-1:
886 a[ndata+i, i*2] = xmax_block*w
887 a[ndata+i, i*2+1] = 1.0*w
888 a[ndata+i, i*2+2] = -xmax_block*w
889 a[ndata+i, i*2+3] = -1.0*w
891 d = num.concatenate((y, num.zeros(n-1)))
892 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2))
894 ynodes = num.zeros(n+1)
895 ynodes[:n] = model[:, 0]*xnodes[:n] + model[:, 1]
896 ynodes[1:] += model[:, 0]*xnodes[1:] + model[:, 1]
897 ynodes[1:n] *= 0.5
899 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2))
901 return xnodes, ynodes, rms_error
904def plf_integrate_piecewise(x_edges, x, y):
905 '''
906 Calculate definite integral of piece-wise linear function on intervals.
908 Use trapezoidal rule to calculate definite integral of a piece-wise linear
909 function for a series of consecutive intervals. ``x_edges`` and ``x`` must
910 be sorted.
912 :param x_edges: array with edges of the intervals
913 :param x,y: arrays with coordinates of piece-wise linear function's
914 control points
915 '''
917 x_all = num.concatenate((x, x_edges))
918 ii = num.argsort(x_all)
919 y_edges = num.interp(x_edges, x, y)
920 y_all = num.concatenate((y, y_edges))
921 xs = x_all[ii]
922 ys = y_all[ii]
923 y_all[ii[1:]] = num.cumsum((xs[1:] - xs[:-1]) * 0.5 * (ys[1:] + ys[:-1]))
924 return num.diff(y_all[-len(y_edges):])
927def diff_fd_1d_4o(dt, data):
928 '''
929 Approximate first derivative of an array (forth order, central FD).
931 :param dt: sampling interval
932 :param data: NumPy array with data samples
934 :returns: NumPy array with same shape as input
936 Interior points are approximated to fourth order, edge points to first
937 order right- or left-sided respectively, points next to edge to second
938 order central.
939 '''
940 import scipy.signal
942 ddata = num.empty_like(data, dtype=float)
944 if data.size >= 5:
945 ddata[2:-2] = scipy.signal.lfilter(
946 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt)
948 if data.size >= 3:
949 ddata[1] = (data[2] - data[0]) / (2. * dt)
950 ddata[-2] = (data[-1] - data[-3]) / (2. * dt)
952 if data.size >= 2:
953 ddata[0] = (data[1] - data[0]) / dt
954 ddata[-1] = (data[-1] - data[-2]) / dt
956 if data.size == 1:
957 ddata[0] = 0.0
959 return ddata
962def diff_fd_1d_2o(dt, data):
963 '''
964 Approximate first derivative of an array (second order, central FD).
966 :param dt: sampling interval
967 :param data: NumPy array with data samples
969 :returns: NumPy array with same shape as input
971 Interior points are approximated to second order, edge points to first
972 order right- or left-sided respectively.
974 Uses :py:func:`numpy.gradient`.
975 '''
977 return num.gradient(data, dt)
980def diff_fd_2d_4o(dt, data):
981 '''
982 Approximate second derivative of an array (forth order, central FD).
984 :param dt: sampling interval
985 :param data: NumPy array with data samples
987 :returns: NumPy array with same shape as input
989 Interior points are approximated to fourth order, next-to-edge points to
990 second order, edge points repeated.
991 '''
992 import scipy.signal
994 ddata = num.empty_like(data, dtype=float)
996 if data.size >= 5:
997 ddata[2:-2] = scipy.signal.lfilter(
998 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2)
1000 if data.size >= 3:
1001 ddata[:2] = (data[2] - 2.0 * data[1] + data[0]) / dt**2
1002 ddata[-2:] = (data[-1] - 2.0 * data[-2] + data[-3]) / dt**2
1004 if data.size < 3:
1005 ddata[:] = 0.0
1007 return ddata
1010def diff_fd_2d_2o(dt, data):
1011 '''
1012 Approximate second derivative of an array (second order, central FD).
1014 :param dt: sampling interval
1015 :param data: NumPy array with data samples
1017 :returns: NumPy array with same shape as input
1019 Interior points are approximated to second order, edge points repeated.
1020 '''
1021 import scipy.signal
1023 ddata = num.empty_like(data, dtype=float)
1025 if data.size >= 3:
1026 ddata[1:-1] = scipy.signal.lfilter(
1027 [1., -2., 1.], [1.], data)[2:] / (dt**2)
1029 ddata[0] = ddata[1]
1030 ddata[-1] = ddata[-2]
1032 if data.size < 3:
1033 ddata[:] = 0.0
1035 return ddata
1038def diff_fd(n, order, dt, data):
1039 '''
1040 Approximate 1st or 2nd derivative of an array.
1042 :param n: 1 for first derivative, 2 for second
1043 :param order: order of the approximation 2 and 4 are supported
1044 :param dt: sampling interval
1045 :param data: NumPy array with data samples
1047 :returns: NumPy array with same shape as input
1049 This is a frontend to the functions :py:func:`diff_fd_1d_2o`,
1050 :py:func:`diff_fd_1d_4o`, :py:func:`diff_fd_2d_2o`, and
1051 :py:func:`diff_fd_2d_4o`.
1053 Raises :py:exc:`ValueError` for unsupported `n` or `order`.
1054 '''
1056 funcs = {
1057 1: {2: diff_fd_1d_2o, 4: diff_fd_1d_4o},
1058 2: {2: diff_fd_2d_2o, 4: diff_fd_2d_4o}}
1060 try:
1061 funcs_n = funcs[n]
1062 except KeyError:
1063 raise ValueError(
1064 'pyrocko.util.diff_fd: '
1065 'Only 1st and 2sd derivatives are supported.')
1067 try:
1068 func = funcs_n[order]
1069 except KeyError:
1070 raise ValueError(
1071 'pyrocko.util.diff_fd: '
1072 'Order %i is not supported for %s derivative. Supported: %s' % (
1073 order, ['', '1st', '2nd'][n],
1074 ', '.join('%i' % order for order in sorted(funcs_n.keys()))))
1076 return func(dt, data)
1079class GlobalVars(object):
1080 reuse_store = dict()
1081 decitab_nmax = 0
1082 decitab = {}
1083 decimate_fir_coeffs = {}
1084 decimate_fir_remez_coeffs = {}
1085 decimate_iir_coeffs = {}
1086 re_frac = None
1089def decimate_coeffs(q, n=None, ftype='iir'):
1091 q = int(q)
1093 if n is None:
1094 if ftype == 'fir':
1095 n = 30
1096 elif ftype == 'fir-remez':
1097 n = 45*q
1098 else:
1099 n = 8
1101 if ftype == 'fir':
1102 coeffs = GlobalVars.decimate_fir_coeffs
1103 if (n, 1./q) not in coeffs:
1104 coeffs[n, 1./q] = signal.firwin(n+1, .75/q, window='hamming')
1106 b = coeffs[n, 1./q]
1107 return b, [1.], n
1109 elif ftype == 'fir-remez':
1110 coeffs = GlobalVars.decimate_fir_remez_coeffs
1111 if (n, 1./q) not in coeffs:
1112 coeffs[n, 1./q] = signal.remez(
1113 n+1, (0., .75/q, 1./q, 1.),
1114 (1., 0.), fs=2, weight=(1, 50))
1115 b = coeffs[n, 1./q]
1116 return b, [1.], n
1118 else:
1119 coeffs = GlobalVars.decimate_iir_coeffs
1120 if (n, 0.05, 0.8/q) not in coeffs:
1121 coeffs[n, 0.05, 0.8/q] = signal.cheby1(n, 0.05, 0.8/q)
1123 b, a = coeffs[n, 0.05, 0.8/q]
1124 return b, a, n
1127def decimate(x, q, n=None, ftype='iir', zi=None, ioff=0):
1128 '''
1129 Downsample the signal x by an integer factor q, using an order n filter
1131 By default, an order 8 Chebyshev type I filter is used or a 30 point FIR
1132 filter with hamming window if ftype is 'fir'.
1134 :param x: the signal to be downsampled (1D :class:`numpy.ndarray`)
1135 :param q: the downsampling factor
1136 :param n: order of the filter (1 less than the length of the filter for a
1137 `fir` filter)
1138 :param ftype: type of the filter; can be `iir`, `fir` or `fir-remez`
1140 :returns: the downsampled signal (1D :class:`numpy.ndarray`)
1142 '''
1144 b, a, n = decimate_coeffs(q, n, ftype)
1146 if zi is None or zi is True:
1147 zi_ = num.zeros(max(len(a), len(b))-1, dtype=float)
1148 else:
1149 zi_ = zi
1151 y, zf = signal.lfilter(b, a, x, zi=zi_)
1153 if zi is not None:
1154 return y[n//2+ioff::q].copy(), zf
1155 else:
1156 return y[n//2+ioff::q].copy()
1159class UnavailableDecimation(Exception):
1160 '''
1161 Exception raised by :py:func:`decitab` for unavailable decimation factors.
1162 '''
1164 pass
1167def gcd(a, b, epsilon=1e-7):
1168 '''
1169 Greatest common divisor.
1170 '''
1172 while b > epsilon*a:
1173 a, b = b, a % b
1175 return a
1178def lcm(a, b):
1179 '''
1180 Least common multiple.
1181 '''
1183 return a*b // gcd(a, b)
1186def mk_decitab(nmax=100):
1187 '''
1188 Make table with decimation sequences.
1190 Decimation from one sampling rate to a lower one is achieved by a
1191 successive application of :py:func:`decimate` with small integer
1192 downsampling factors (because using large downsampling factors can make the
1193 decimation unstable or slow). This function sets up a table with downsample
1194 sequences for factors up to ``nmax``.
1195 '''
1197 tab = GlobalVars.decitab
1198 for i in range(1, 10):
1199 for j in range(1, i+1):
1200 for k in range(1, j+1):
1201 for l_ in range(1, k+1):
1202 for m in range(1, l_+1):
1203 p = i*j*k*l_*m
1204 if p > nmax:
1205 break
1206 if p not in tab:
1207 tab[p] = (i, j, k, l_, m)
1208 if i*j*k*l_ > nmax:
1209 break
1210 if i*j*k > nmax:
1211 break
1212 if i*j > nmax:
1213 break
1214 if i > nmax:
1215 break
1217 GlobalVars.decitab_nmax = nmax
1220def zfmt(n):
1221 return '%%0%ii' % (int(math.log10(n - 1)) + 1)
1224def _year_to_time(year):
1225 tt = (year, 1, 1, 0, 0, 0)
1226 return to_time_float(calendar.timegm(tt))
1229def _working_year(year):
1230 try:
1231 tt = (year, 1, 1, 0, 0, 0)
1232 t = calendar.timegm(tt)
1233 tt2_ = time.gmtime(t)
1234 tt2 = tuple(tt2_)[:6]
1235 if tt != tt2:
1236 return False
1238 s = '%i-01-01 00:00:00' % year
1239 s2 = time.strftime('%Y-%m-%d %H:%M:%S', tt2_)
1240 if s != s2:
1241 return False
1243 t3 = str_to_time(s2, format='%Y-%m-%d %H:%M:%S')
1244 s3 = time_to_str(t3, format='%Y-%m-%d %H:%M:%S')
1245 if s3 != s2:
1246 return False
1248 except Exception:
1249 return False
1251 return True
1254def working_system_time_range(year_min_lim=None, year_max_lim=None):
1255 '''
1256 Check time range supported by the systems's time conversion functions.
1258 Returns system time stamps of start of year of first/last fully supported
1259 year span. If this is before 1900 or after 2100, return first/last century
1260 which is fully supported.
1262 :returns: ``(tmin, tmax, year_min, year_max)``
1263 '''
1265 year0 = 2000
1266 year_min = year0
1267 year_max = year0
1269 itests = list(range(101))
1270 for i in range(19):
1271 itests.append(200 + i*100)
1273 for i in itests:
1274 year = year0 - i
1275 if year_min_lim is not None and year < year_min_lim:
1276 year_min = year_min_lim
1277 break
1278 elif not _working_year(year):
1279 break
1280 else:
1281 year_min = year
1283 for i in itests:
1284 year = year0 + i
1285 if year_max_lim is not None and year > year_max_lim:
1286 year_max = year_max_lim
1287 break
1288 elif not _working_year(year + 1):
1289 break
1290 else:
1291 year_max = year
1293 return (
1294 _year_to_time(year_min),
1295 _year_to_time(year_max),
1296 year_min, year_max)
1299g_working_system_time_range = None
1302def get_working_system_time_range():
1303 '''
1304 Caching variant of :py:func:`working_system_time_range`.
1305 '''
1307 global g_working_system_time_range
1308 if g_working_system_time_range is None:
1309 g_working_system_time_range = working_system_time_range()
1311 return g_working_system_time_range
1314def is_working_time(t):
1315 tmin, tmax, _, _ = get_working_system_time_range()
1316 return tmin <= t <= tmax
1319def julian_day_of_year(timestamp):
1320 '''
1321 Get the day number after the 1st of January of year in ``timestamp``.
1323 :returns: day number as int
1324 '''
1326 return time.gmtime(int(timestamp)).tm_yday
1329def hour_start(timestamp):
1330 '''
1331 Get beginning of hour for any point in time.
1333 :param timestamp: time instant as system timestamp (in seconds)
1335 :returns: instant of hour start as system timestamp
1336 '''
1338 tt = time.gmtime(int(timestamp))
1339 tts = tt[0:4] + (0, 0)
1340 return to_time_float(calendar.timegm(tts))
1343def day_start(timestamp):
1344 '''
1345 Get beginning of day for any point in time.
1347 :param timestamp: time instant as system timestamp (in seconds)
1349 :returns: instant of day start as system timestamp
1350 '''
1352 tt = time.gmtime(int(timestamp))
1353 tts = tt[0:3] + (0, 0, 0)
1354 return to_time_float(calendar.timegm(tts))
1357def month_start(timestamp):
1358 '''
1359 Get beginning of month for any point in time.
1361 :param timestamp: time instant as system timestamp (in seconds)
1363 :returns: instant of month start as system timestamp
1364 '''
1366 tt = time.gmtime(int(timestamp))
1367 tts = tt[0:2] + (1, 0, 0, 0)
1368 return to_time_float(calendar.timegm(tts))
1371def year_start(timestamp):
1372 '''
1373 Get beginning of year for any point in time.
1375 :param timestamp: time instant as system timestamp (in seconds)
1377 :returns: instant of year start as system timestamp
1378 '''
1380 tt = time.gmtime(int(timestamp))
1381 tts = tt[0:1] + (1, 1, 0, 0, 0)
1382 return to_time_float(calendar.timegm(tts))
1385def iter_days(tmin, tmax):
1386 '''
1387 Yields begin and end of days until given time span is covered.
1389 :param tmin,tmax: input time span
1391 :yields: tuples with (begin, end) of days as system timestamps
1392 '''
1394 t = day_start(tmin)
1395 while t < tmax:
1396 tend = day_start(t + 26*60*60)
1397 yield t, tend
1398 t = tend
1401def iter_months(tmin, tmax):
1402 '''
1403 Yields begin and end of months until given time span is covered.
1405 :param tmin,tmax: input time span
1407 :yields: tuples with (begin, end) of months as system timestamps
1408 '''
1410 t = month_start(tmin)
1411 while t < tmax:
1412 tend = month_start(t + 24*60*60*33)
1413 yield t, tend
1414 t = tend
1417def iter_years(tmin, tmax):
1418 '''
1419 Yields begin and end of years until given time span is covered.
1421 :param tmin,tmax: input time span
1423 :yields: tuples with (begin, end) of years as system timestamps
1424 '''
1426 t = year_start(tmin)
1427 while t < tmax:
1428 tend = year_start(t + 24*60*60*369)
1429 yield t, tend
1430 t = tend
1433def today():
1434 return day_start(time.time())
1437def tomorrow():
1438 return day_start(time.time() + 24*60*60)
1441def decitab(n):
1442 '''
1443 Get integer decimation sequence for given downampling factor.
1445 :param n: target decimation factor
1447 :returns: tuple with downsampling sequence
1448 '''
1450 if n > GlobalVars.decitab_nmax:
1451 mk_decitab(n*2)
1452 if n not in GlobalVars.decitab:
1453 raise UnavailableDecimation('ratio = %g' % n)
1454 return GlobalVars.decitab[n]
1457def ctimegm(s, format='%Y-%m-%d %H:%M:%S'):
1458 '''
1459 Convert string representing UTC time to system time.
1461 :param s: string to be interpreted
1462 :param format: format string passed to :py:func:`time.strptime`
1464 :returns: system time stamp
1466 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime.
1468 .. note::
1469 This function is to be replaced by :py:func:`str_to_time`.
1470 '''
1472 return calendar.timegm(time.strptime(s, format))
1475def gmctime(t, format='%Y-%m-%d %H:%M:%S'):
1476 '''
1477 Get string representation from system time, UTC.
1479 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime.
1481 .. note::
1482 This function is to be repaced by :py:func:`time_to_str`.
1483 '''
1485 return time.strftime(format, time.gmtime(t))
1488def gmctime_v(t, format='%a, %d %b %Y %H:%M:%S'):
1489 '''
1490 Get string representation from system time, UTC. Same as
1491 :py:func:`gmctime` but with a more verbose default format.
1493 .. note::
1494 This function is to be replaced by :py:func:`time_to_str`.
1495 '''
1497 return time.strftime(format, time.gmtime(t))
1500def gmctime_fn(t, format='%Y-%m-%d_%H-%M-%S'):
1501 '''
1502 Get string representation from system time, UTC. Same as
1503 :py:func:`gmctime` but with a default usable in filenames.
1505 .. note::
1506 This function is to be replaced by :py:func:`time_to_str`.
1507 '''
1509 return time.strftime(format, time.gmtime(t))
1512class TimeStrError(Exception):
1513 '''
1514 Raised for invalid time strings.
1515 '''
1516 pass
1519class FractionalSecondsMissing(TimeStrError):
1520 '''
1521 Exception raised by :py:func:`str_to_time` when the given string lacks
1522 fractional seconds.
1523 '''
1525 pass
1528class FractionalSecondsWrongNumberOfDigits(TimeStrError):
1529 '''
1530 Exception raised by :py:func:`str_to_time` when the given string has an
1531 incorrect number of digits in the fractional seconds part.
1532 '''
1534 pass
1537def _endswith_n(s, endings):
1538 for ix, x in enumerate(endings):
1539 if s.endswith(x):
1540 return ix
1541 return -1
1544def str_to_time(s, format='%Y-%m-%d %H:%M:%S.OPTFRAC'):
1545 '''
1546 Convert string representing UTC time to floating point system time.
1548 :param s: string representing UTC time
1549 :param format: time string format
1550 :returns: system time stamp as floating point value
1552 Uses the semantics of :py:func:`time.strptime` but allows for fractional
1553 seconds. If the format ends with ``'.FRAC'``, anything after a dot is
1554 interpreted as fractional seconds. If the format ends with ``'.OPTFRAC'``,
1555 the fractional part, including the dot is made optional. The latter has the
1556 consequence, that the time strings and the format may not contain any other
1557 dots. If the format ends with ``'.xFRAC'`` where x is 1, 2, or 3, it is
1558 ensured, that exactly that number of digits are present in the fractional
1559 seconds.
1560 '''
1562 time_float = get_time_float()
1564 if util_ext is not None:
1565 try:
1566 t, tfrac = util_ext.stt(s, format)
1567 except util_ext.UtilExtError as e:
1568 raise TimeStrError(
1569 '%s, string=%s, format=%s' % (str(e), s, format))
1571 return time_float(t)+tfrac
1573 fracsec = 0.
1574 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC'
1576 iend = _endswith_n(format, fixed_endings)
1577 if iend != -1:
1578 dotpos = s.rfind('.')
1579 if dotpos == -1:
1580 raise FractionalSecondsMissing(
1581 'string=%s, format=%s' % (s, format))
1583 if iend > 0 and iend != (len(s)-dotpos-1):
1584 raise FractionalSecondsWrongNumberOfDigits(
1585 'string=%s, format=%s' % (s, format))
1587 format = format[:-len(fixed_endings[iend])]
1588 fracsec = float(s[dotpos:])
1589 s = s[:dotpos]
1591 elif format.endswith('.OPTFRAC'):
1592 dotpos = s.rfind('.')
1593 format = format[:-8]
1594 if dotpos != -1 and len(s[dotpos:]) > 1:
1595 fracsec = float(s[dotpos:])
1597 if dotpos != -1:
1598 s = s[:dotpos]
1600 try:
1601 return time_float(calendar.timegm(time.strptime(s, format))) \
1602 + fracsec
1603 except ValueError as e:
1604 raise TimeStrError('%s, string=%s, format=%s' % (str(e), s, format))
1607stt = str_to_time
1610def str_to_time_fillup(s):
1611 '''
1612 Default :py:func:`str_to_time` with filling in of missing values.
1614 Allows e.g. `'2010-01-01 00:00:00'` as `'2010-01-01 00:00'`,
1615 `'2010-01-01 00'`, ..., or `'2010'`.
1616 '''
1618 if s == 'now':
1619 return time.time()
1620 elif s == 'today':
1621 return today()
1622 elif s == 'tomorrow':
1623 return tomorrow()
1625 if len(s) in (4, 7, 10, 13, 16):
1626 s += '0000-01-01 00:00:00'[len(s):]
1628 return str_to_time(s)
1631def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'):
1632 '''
1633 Get string representation for floating point system time.
1635 :param t: floating point system time
1636 :param format: time string format
1637 :returns: string representing UTC time
1639 Uses the semantics of :py:func:`time.strftime` but additionally allows for
1640 fractional seconds. If ``format`` contains ``'.xFRAC'``, where ``x`` is a
1641 digit between 1 and 9, this is replaced with the fractional part of ``t``
1642 with ``x`` digits precision.
1643 '''
1645 if pyrocko.grumpy > 0:
1646 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise')
1648 if isinstance(format, int):
1649 if format > 0:
1650 format = '%Y-%m-%d %H:%M:%S.' + '%iFRAC' % format
1651 else:
1652 format = '%Y-%m-%d %H:%M:%S'
1654 if util_ext is not None:
1655 t0 = num.floor(t)
1656 try:
1657 return util_ext.tts(int(t0), float(t - t0), format)
1658 except util_ext.UtilExtError as e:
1659 raise TimeStrError(
1660 '%s, timestamp=%f, format=%s' % (str(e), t, format))
1662 if not GlobalVars.re_frac:
1663 GlobalVars.re_frac = re.compile(r'\.[1-9]FRAC')
1664 GlobalVars.frac_formats = dict(
1665 [('.%sFRAC' % x, '%.'+x+'f') for x in '123456789'])
1667 ts = float(num.floor(t))
1668 tfrac = t-ts
1670 m = GlobalVars.re_frac.search(format)
1671 if m:
1672 sfrac = (GlobalVars.frac_formats[m.group(0)] % tfrac)
1673 if sfrac[0] == '1':
1674 ts += 1.
1676 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1)
1678 return time.strftime(format, time.gmtime(ts))
1681tts = time_to_str
1682_abbr_weekday = 'Mon Tue Wed Thu Fri Sat Sun'.split()
1683_abbr_month = 'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split()
1686def mystrftime(fmt=None, tt=None, milliseconds=0):
1687 # Needed by Snuffler for the time axis. In other cases `time_to_str`
1688 # should be used.
1690 if fmt is None:
1691 fmt = '%Y-%m-%d %H:%M:%S .%r'
1693 # Get these two locale independent, needed by Snuffler.
1694 # Setting the locale seems to have no effect.
1695 fmt = fmt.replace('%a', _abbr_weekday[tt.tm_wday])
1696 fmt = fmt.replace('%b', _abbr_month[tt.tm_mon-1])
1698 fmt = fmt.replace('%r', '%03i' % int(round(milliseconds)))
1699 fmt = fmt.replace('%u', '%06i' % int(round(milliseconds*1000)))
1700 fmt = fmt.replace('%n', '%09i' % int(round(milliseconds*1000000)))
1701 return time.strftime(fmt, tt)
1704def gmtime_x(timestamp):
1705 etimestamp = float(num.floor(timestamp))
1706 tt = time.gmtime(etimestamp)
1707 ms = (timestamp-etimestamp)*1000
1708 return tt, ms
1711def plural_s(n):
1712 if not isinstance(n, int):
1713 n = len(n)
1715 return 's' if n != 1 else ''
1718def ensuredirs(dst):
1719 '''
1720 Create all intermediate path components for a target path.
1722 :param dst: target path
1724 The leaf part of the target path is not created (use :py:func:`ensuredir`
1725 if a the target path is a directory to be created).
1726 '''
1728 d, x = os.path.split(dst.rstrip(os.sep))
1729 dirs = []
1730 while d and not os.path.exists(d):
1731 dirs.append(d)
1732 d, x = os.path.split(d)
1734 dirs.reverse()
1736 for d in dirs:
1737 try:
1738 os.mkdir(d)
1739 except OSError as e:
1740 if not e.errno == errno.EEXIST:
1741 raise
1744def ensuredir(dst):
1745 '''
1746 Create directory and all intermediate path components to it as needed.
1748 :param dst: directory name
1750 Nothing is done if the given target already exists.
1751 '''
1753 if os.path.exists(dst):
1754 return
1756 dst.rstrip(os.sep)
1758 ensuredirs(dst)
1759 try:
1760 os.mkdir(dst)
1761 except OSError as e:
1762 if not e.errno == errno.EEXIST:
1763 raise
1766def reuse(x):
1767 '''
1768 Get unique instance of an object.
1770 :param x: hashable object
1771 :returns: reference to x or an equivalent object
1773 Cache object ``x`` in a global dict for reuse, or if x already
1774 is in that dict, return a reference to it.
1775 '''
1777 grs = GlobalVars.reuse_store
1778 if x not in grs:
1779 grs[x] = x
1780 return grs[x]
1783def deuse(x):
1784 grs = GlobalVars.reuse_store
1785 if x in grs:
1786 del grs[x]
1789class Anon(object):
1790 '''
1791 Dict-to-object utility.
1793 Any given arguments are stored as attributes.
1795 Example::
1797 a = Anon(x=1, y=2)
1798 print a.x, a.y
1799 '''
1801 def __init__(self, **dict):
1802 for k in dict:
1803 self.__dict__[k] = dict[k]
1806def iter_select_files(
1807 paths,
1808 include=None,
1809 exclude=None,
1810 selector=None,
1811 show_progress=True,
1812 pass_through=None):
1814 '''
1815 Recursively select files (generator variant).
1817 See :py:func:`select_files`.
1818 '''
1820 if show_progress:
1821 progress_beg('selecting files...')
1823 ngood = 0
1824 if include is not None:
1825 rinclude = re.compile(include)
1827 def check_include(path):
1828 m = rinclude.search(path)
1829 if not m:
1830 return False
1832 if selector is None:
1833 return True
1835 infos = Anon(**m.groupdict())
1836 return selector(infos)
1837 else:
1838 check_include = None
1840 if exclude is not None:
1841 rexclude = re.compile(exclude)
1843 def check_exclude(path):
1844 return not bool(rexclude.search(path))
1845 else:
1846 check_exclude = None
1848 if check_include and check_exclude:
1850 def check(path):
1851 return check_include(path) and check_exclude(path)
1853 elif check_include:
1854 check = check_include
1856 elif check_exclude:
1857 check = check_exclude
1859 else:
1860 check = None
1862 if isinstance(paths, str):
1863 paths = [paths]
1865 for path in paths:
1866 if pass_through and pass_through(path):
1867 if check is None or check(path):
1868 yield path
1870 elif os.path.isdir(path):
1871 for (dirpath, dirnames, filenames) in os.walk(path):
1872 dirnames.sort()
1873 filenames.sort()
1874 for filename in filenames:
1875 path = op.join(dirpath, filename)
1876 if check is None or check(path):
1877 yield os.path.abspath(path)
1878 ngood += 1
1879 else:
1880 if check is None or check(path):
1881 yield os.path.abspath(path)
1882 ngood += 1
1884 if show_progress:
1885 progress_end('%i file%s selected.' % (ngood, plural_s(ngood)))
1888def select_files(
1889 paths, include=None, exclude=None, selector=None, show_progress=True,
1890 regex=None):
1892 '''
1893 Recursively select files.
1895 :param paths: entry path names
1896 :param include: pattern for conditional inclusion
1897 :param exclude: pattern for conditional exclusion
1898 :param selector: callback for conditional inclusion
1899 :param show_progress: if True, indicate start and stop of processing
1900 :param regex: alias for ``include`` (backwards compatibility)
1901 :returns: list of path names
1903 Recursively finds all files under given entry points ``paths``. If
1904 parameter ``include`` is a regular expression, only files with matching
1905 path names are included. If additionally parameter ``selector`` is given a
1906 callback function, only files for which the callback returns ``True`` are
1907 included. The callback should take a single argument. The callback is
1908 called with a single argument, an object, having as attributes, any named
1909 groups given in ``include``.
1911 Examples
1913 To find all files ending in ``'.mseed'`` or ``'.msd'``::
1915 select_files(paths,
1916 include=r'\\.(mseed|msd)$')
1918 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for
1919 the year::
1921 select_files(paths,
1922 include=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$',
1923 selector=(lambda x: int(x.year) == 2009))
1924 '''
1925 if regex is not None:
1926 assert include is None
1927 include = regex
1929 return list(iter_select_files(
1930 paths, include, exclude, selector, show_progress))
1933def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
1934 '''
1935 Convert positive integer to a base36 string.
1936 '''
1938 if not isinstance(number, (int, long)):
1939 raise TypeError('number must be an integer')
1940 if number < 0:
1941 raise ValueError('number must be positive')
1943 # Special case for small numbers
1944 if number < 36:
1945 return alphabet[number]
1947 base36 = ''
1948 while number != 0:
1949 number, i = divmod(number, 36)
1950 base36 = alphabet[i] + base36
1952 return base36
1955def base36decode(number):
1956 '''
1957 Decode base36 endcoded positive integer.
1958 '''
1960 return int(number, 36)
1963class UnpackError(Exception):
1964 '''
1965 Exception raised when :py:func:`unpack_fixed` encounters an error.
1966 '''
1968 pass
1971ruler = ''.join(['%-10i' % i for i in range(8)]) \
1972 + '\n' + '0123456789' * 8 + '\n'
1975def unpack_fixed(format, line, *callargs):
1976 '''
1977 Unpack fixed format string, as produced by many fortran codes.
1979 :param format: format specification
1980 :param line: string to be processed
1981 :param callargs: callbacks for callback fields in the format
1983 The format is described by a string of comma-separated fields. Each field
1984 is defined by a character for the field type followed by the field width. A
1985 questionmark may be appended to the field description to allow the argument
1986 to be optional (The data string is then allowed to be filled with blanks
1987 and ``None`` is returned in this case).
1989 The following field types are available:
1991 ==== ================================================================
1992 Type Description
1993 ==== ================================================================
1994 A string (full field width is extracted)
1995 a string (whitespace at the beginning and the end is removed)
1996 i integer value
1997 f floating point value
1998 @ special type, a callback must be given for the conversion
1999 x special field type to skip parts of the string
2000 ==== ================================================================
2001 '''
2003 ipos = 0
2004 values = []
2005 icall = 0
2006 for form in format.split(','):
2007 form = form.strip()
2008 optional = form[-1] == '?'
2009 form = form.rstrip('?')
2010 typ = form[0]
2011 ln = int(form[1:])
2012 s = line[ipos:ipos+ln]
2013 cast = {
2014 'x': None,
2015 'A': str,
2016 'a': lambda x: x.strip(),
2017 'i': int,
2018 'f': float,
2019 '@': 'extra'}[typ]
2021 if cast == 'extra':
2022 cast = callargs[icall]
2023 icall += 1
2025 if cast is not None:
2026 if optional and s.strip() == '':
2027 values.append(None)
2028 else:
2029 try:
2030 values.append(cast(s))
2031 except Exception:
2032 mark = [' '] * 80
2033 mark[ipos:ipos+ln] = ['^'] * ln
2034 mark = ''.join(mark)
2035 raise UnpackError(
2036 'Invalid cast to type "%s" at position [%i:%i] of '
2037 'line: \n%s%s\n%s' % (
2038 typ, ipos, ipos+ln, ruler, line.rstrip(), mark))
2040 ipos += ln
2042 return values
2045_pattern_cache = {}
2048def _nslc_pattern(pattern):
2049 if pattern not in _pattern_cache:
2050 rpattern = re.compile(fnmatch.translate(pattern), re.I)
2051 _pattern_cache[pattern] = rpattern
2052 else:
2053 rpattern = _pattern_cache[pattern]
2055 return rpattern
2058def match_nslc(patterns, nslc):
2059 '''
2060 Match network-station-location-channel code against pattern or list of
2061 patterns.
2063 :param patterns: pattern or list of patterns
2064 :param nslc: tuple with (network, station, location, channel) as strings
2066 :returns: ``True`` if the pattern matches or if any of the given patterns
2067 match; or ``False``.
2069 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq].
2071 Example::
2073 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True
2074 '''
2076 if isinstance(patterns, str):
2077 patterns = [patterns]
2079 if not isinstance(nslc, str):
2080 s = '.'.join(nslc)
2081 else:
2082 s = nslc
2084 for pattern in patterns:
2085 if _nslc_pattern(pattern).match(s):
2086 return True
2088 return False
2091def match_nslcs(patterns, nslcs):
2092 '''
2093 Get network-station-location-channel codes that match given pattern or any
2094 of several given patterns.
2096 :param patterns: pattern or list of patterns
2097 :param nslcs: list of (network, station, location, channel) tuples
2099 See also :py:func:`match_nslc`
2100 '''
2102 matching = []
2103 for nslc in nslcs:
2104 if match_nslc(patterns, nslc):
2105 matching.append(nslc)
2107 return matching
2110def glob_filter(patterns, names):
2111 '''
2112 Select names matching any of the given patterns.
2113 '''
2115 if not patterns:
2116 return names
2118 rpattern = re.compile(r'|'.join(
2119 fnmatch.translate(pattern) for pattern in patterns))
2121 return [name for name in names if rpattern.match(name)]
2124class Timeout(Exception):
2125 pass
2128def create_lockfile(fn, timeout=None, timewarn=10.):
2129 t0 = time.time()
2131 while True:
2132 try:
2133 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL)
2134 os.close(f)
2135 return
2137 except OSError as e:
2138 if e.errno in (errno.EEXIST, 13): # 13 occurs on windows
2139 pass # retry
2140 else:
2141 raise
2143 tnow = time.time()
2145 if timeout is not None and tnow - t0 > timeout:
2146 raise Timeout(
2147 'Timeout (%gs) occured while waiting to get exclusive '
2148 'access to: %s' % (timeout, fn))
2150 if timewarn is not None and tnow - t0 > timewarn:
2151 logger.warning(
2152 'Waiting since %gs to get exclusive access to: %s' % (
2153 timewarn, fn))
2155 timewarn *= 2
2157 time.sleep(0.01)
2160def delete_lockfile(fn):
2161 os.unlink(fn)
2164class Lockfile(Exception):
2166 def __init__(self, path, timeout=5, timewarn=10.):
2167 self._path = path
2168 self._timeout = timeout
2169 self._timewarn = timewarn
2171 def __enter__(self):
2172 create_lockfile(
2173 self._path, timeout=self._timeout, timewarn=self._timewarn)
2174 return None
2176 def __exit__(self, type, value, traceback):
2177 delete_lockfile(self._path)
2180class SoleError(Exception):
2181 '''
2182 Exception raised by objects of type :py:class:`Sole`, when an concurrent
2183 instance is running.
2184 '''
2186 pass
2189class Sole(object):
2190 '''
2191 Use POSIX advisory file locking to ensure that only a single instance of a
2192 program is running.
2194 :param pid_path: path to lockfile to be used
2196 Usage::
2198 from pyrocko.util import Sole, SoleError, setup_logging
2199 import os
2201 setup_logging('my_program')
2203 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock')
2204 try:
2205 sole = Sole(pid_path)
2207 except SoleError, e:
2208 logger.fatal( str(e) )
2209 sys.exit(1)
2210 '''
2212 def __init__(self, pid_path):
2213 self._pid_path = pid_path
2214 self._other_running = False
2215 ensuredirs(self._pid_path)
2216 self._lockfile = None
2217 self._os = os
2219 if not fcntl:
2220 raise SoleError(
2221 'Python standard library module "fcntl" not available on '
2222 'this platform.')
2224 self._fcntl = fcntl
2226 try:
2227 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY)
2228 except Exception:
2229 raise SoleError(
2230 'Cannot open lockfile (path = %s)' % self._pid_path)
2232 try:
2233 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB)
2235 except IOError:
2236 self._other_running = True
2237 try:
2238 f = open(self._pid_path, 'r')
2239 pid = f.read().strip()
2240 f.close()
2241 except Exception:
2242 pid = '?'
2244 raise SoleError('Other instance is running (pid = %s)' % pid)
2246 try:
2247 os.ftruncate(self._lockfile, 0)
2248 os.write(self._lockfile, '%i\n' % os.getpid())
2249 os.fsync(self._lockfile)
2251 except Exception:
2252 # the pid is only stored for user information, so this is allowed
2253 # to fail
2254 pass
2256 def __del__(self):
2257 if not self._other_running:
2258 if self._lockfile is not None:
2259 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN)
2260 self._os.close(self._lockfile)
2261 try:
2262 self._os.unlink(self._pid_path)
2263 except Exception:
2264 pass
2267re_escapequotes = re.compile(r"(['\\])")
2270def escapequotes(s):
2271 return re_escapequotes.sub(r'\\\1', s)
2274class TableWriter(object):
2275 '''
2276 Write table of space separated values to a file.
2278 :param f: file like object
2280 Strings containing spaces are quoted on output.
2281 '''
2283 def __init__(self, f):
2284 self._f = f
2286 def writerow(self, row, minfieldwidths=None):
2288 '''
2289 Write one row of values to underlying file.
2291 :param row: iterable of values
2292 :param minfieldwidths: minimum field widths for the values
2294 Each value in in ``row`` is converted to a string and optionally padded
2295 with blanks. The resulting strings are output separated with blanks. If
2296 any values given are strings and if they contain whitespace, they are
2297 quoted with single quotes, and any internal single quotes are
2298 backslash-escaped.
2299 '''
2301 out = []
2303 for i, x in enumerate(row):
2304 w = 0
2305 if minfieldwidths and i < len(minfieldwidths):
2306 w = minfieldwidths[i]
2308 if isinstance(x, str):
2309 if re.search(r"\s|'", x):
2310 x = "'%s'" % escapequotes(x)
2312 x = x.ljust(w)
2313 else:
2314 x = str(x).rjust(w)
2316 out.append(x)
2318 self._f.write(' '.join(out).rstrip() + '\n')
2321class TableReader(object):
2323 '''
2324 Read table of space separated values from a file.
2326 :param f: file-like object
2328 This uses Pythons shlex module to tokenize lines. Should deal correctly
2329 with quoted strings.
2330 '''
2332 def __init__(self, f):
2333 self._f = f
2334 self.eof = False
2336 def readrow(self):
2337 '''
2338 Read one row from the underlying file, tokenize it with shlex.
2340 :returns: tokenized line as a list of strings.
2341 '''
2343 line = self._f.readline()
2344 if not line:
2345 self.eof = True
2346 return []
2347 line.strip()
2348 if line.startswith('#'):
2349 return []
2351 return qsplit(line)
2354def gform(number, significant_digits=3):
2355 '''
2356 Pretty print floating point numbers.
2358 Align floating point numbers at the decimal dot.
2360 ::
2362 | -d.dde+xxx|
2363 | -d.dde+xx |
2364 |-ddd. |
2365 | -dd.d |
2366 | -d.dd |
2367 | -0.ddd |
2368 | -0.0ddd |
2369 | -0.00ddd |
2370 | -d.dde-xx |
2371 | -d.dde-xxx|
2372 | nan|
2375 The formatted string has length ``significant_digits * 2 + 6``.
2376 '''
2378 no_exp_range = (pow(10., -1),
2379 pow(10., significant_digits))
2380 width = significant_digits+significant_digits-1+1+1+5
2382 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.:
2383 s = ('%#.*g' % (significant_digits, number)).rstrip('0')
2384 else:
2385 s = '%.*E' % (significant_digits-1, number)
2386 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width)
2387 if s.strip().lower() == 'nan':
2388 s = 'nan'.rjust(width)
2389 return s
2392def human_bytesize(value):
2394 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split()
2396 if value == 1:
2397 return '1 Byte'
2399 for i, ext in enumerate(exts):
2400 x = float(value) / 1000**i
2401 if round(x) < 10. and not value < 1000:
2402 return '%.1f %s' % (x, ext)
2403 if round(x) < 1000.:
2404 return '%.0f %s' % (x, ext)
2406 return '%i Bytes' % value
2409re_compatibility = re.compile(
2410 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' +
2411 r'(dummy|poel|qseis|qssp))\.'
2412)
2415def pf_is_old(fn):
2416 oldstyle = False
2417 with open(fn, 'r') as f:
2418 for line in f:
2419 if re_compatibility.search(line):
2420 oldstyle = True
2422 return oldstyle
2425def pf_upgrade(fn):
2426 need = pf_is_old(fn)
2427 if need:
2428 fn_temp = fn + '.temp'
2430 with open(fn, 'r') as fin:
2431 with open(fn_temp, 'w') as fout:
2432 for line in fin:
2433 line = re_compatibility.sub('!pf.', line)
2434 fout.write(line)
2436 os.rename(fn_temp, fn)
2438 return need
2441def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'):
2442 '''
2443 Extract leap second information from tzdata.
2445 Based on example at http://stackoverflow.com/questions/19332902/\
2446 extract-historic-leap-seconds-from-tzdata
2448 See also 'man 5 tzfile'.
2449 '''
2451 from struct import unpack, calcsize
2452 out = []
2453 with open(tzfile, 'rb') as f:
2454 # read header
2455 fmt = '>4s c 15x 6l'
2456 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt,
2457 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt)))
2458 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file'
2460 # skip over some uninteresting data
2461 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict(
2462 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt)
2463 f.read(calcsize(fmt))
2465 # read leap-seconds
2466 fmt = '>2l'
2467 for i in range(leapcnt):
2468 tleap, nleap = unpack(fmt, f.read(calcsize(fmt)))
2469 out.append((tleap-nleap+1, nleap))
2471 return out
2474class LeapSecondsError(Exception):
2475 pass
2478class LeapSecondsOutdated(LeapSecondsError):
2479 pass
2482class InvalidLeapSecondsFile(LeapSecondsOutdated):
2483 pass
2486def parse_leap_seconds_list(fn):
2487 data = []
2488 texpires = None
2489 try:
2490 t0 = int(round(str_to_time('1900-01-01 00:00:00')))
2491 except TimeStrError:
2492 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800
2494 tnow = int(round(time.time()))
2496 if not op.exists(fn):
2497 raise LeapSecondsOutdated('no leap seconds file found')
2499 try:
2500 with open(fn, 'rb') as f:
2501 for line in f:
2502 if line.strip().startswith(b'<!DOCTYPE'):
2503 raise InvalidLeapSecondsFile('invalid leap seconds file')
2505 if line.startswith(b'#@'):
2506 texpires = int(line.split()[1]) + t0
2507 elif line.startswith(b'#') or len(line) < 5:
2508 pass
2509 else:
2510 toks = line.split()
2511 t = int(toks[0]) + t0
2512 nleap = int(toks[1]) - 10
2513 data.append((t, nleap))
2515 except IOError:
2516 raise LeapSecondsError('cannot read leap seconds file %s' % fn)
2518 if texpires is None or tnow > texpires:
2519 raise LeapSecondsOutdated('leap seconds list is outdated')
2521 return data
2524def read_leap_seconds2():
2525 from pyrocko import config
2526 conf = config.config()
2527 fn = conf.leapseconds_path
2528 url = conf.leapseconds_url
2529 # check for outdated default URL
2530 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list':
2531 url = 'https://www.ietf.org/timezones/data/leap-seconds.list'
2532 logger.info(
2533 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2535 if url == 'https://www.ietf.org/timezones/data/leap-seconds.list':
2536 url = 'https://hpiers.obspm.fr/iers/bul/bulc/ntp/leap-seconds.list'
2537 logger.info(
2538 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2540 for i in range(3):
2541 try:
2542 return parse_leap_seconds_list(fn)
2544 except LeapSecondsOutdated:
2545 try:
2546 logger.info('updating leap seconds list...')
2547 download_file(url, fn)
2549 except Exception as e:
2550 raise LeapSecondsError(
2551 'cannot download leap seconds list from %s to %s (%s)'
2552 % (url, fn, e))
2554 raise LeapSecondsError('Could not retrieve/read leap seconds file.')
2557def gps_utc_offset(t_utc):
2558 '''
2559 Time offset t_gps - t_utc for a given t_utc.
2560 '''
2561 ls = read_leap_seconds2()
2562 i = 0
2563 if t_utc < ls[0][0]:
2564 return ls[0][1] - 1 - 9
2566 while i < len(ls) - 1:
2567 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]:
2568 return ls[i][1] - 9
2569 i += 1
2571 return ls[-1][1] - 9
2574def utc_gps_offset(t_gps):
2575 '''
2576 Time offset t_utc - t_gps for a given t_gps.
2577 '''
2578 ls = read_leap_seconds2()
2580 if t_gps < ls[0][0] + ls[0][1] - 9:
2581 return - (ls[0][1] - 1 - 9)
2583 i = 0
2584 while i < len(ls) - 1:
2585 if ls[i][0] + ls[i][1] - 9 <= t_gps \
2586 and t_gps < ls[i+1][0] + ls[i+1][1] - 9:
2587 return - (ls[i][1] - 9)
2588 i += 1
2590 return - (ls[-1][1] - 9)
2593def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'):
2594 import itertools
2595 import glob
2596 from pyrocko.io.io_common import FileLoadError
2598 def iload_filename(filename, **kwargs):
2599 try:
2600 with open(filename, 'rb') as f:
2601 for cr in iload_fh(f, **kwargs):
2602 yield cr
2604 except FileLoadError as e:
2605 e.set_context('filename', filename)
2606 raise
2608 def iload_dirname(dirname, **kwargs):
2609 for entry in os.listdir(dirname):
2610 fpath = op.join(dirname, entry)
2611 if op.isfile(fpath):
2612 for cr in iload_filename(fpath, **kwargs):
2613 yield cr
2615 def iload_glob(pattern, **kwargs):
2617 for fn in glob.iglob(pattern):
2618 for cr in iload_filename(fn, **kwargs):
2619 yield cr
2621 def iload(source, **kwargs):
2622 if isinstance(source, str):
2623 if op.isdir(source):
2624 return iload_dirname(source, **kwargs)
2625 elif op.isfile(source):
2626 return iload_filename(source, **kwargs)
2627 else:
2628 return iload_glob(source, **kwargs)
2630 elif hasattr(source, 'read'):
2631 return iload_fh(source, **kwargs)
2632 else:
2633 return itertools.chain.from_iterable(
2634 iload(subsource, **kwargs) for subsource in source)
2636 iload_filename.__doc__ = '''
2637 Read %s information from named file.
2638 ''' % doc_fmt
2640 iload_dirname.__doc__ = '''
2641 Read %s information from directory of %s files.
2642 ''' % (doc_fmt, doc_fmt)
2644 iload_glob.__doc__ = '''
2645 Read %s information from files matching a glob pattern.
2646 ''' % doc_fmt
2648 iload.__doc__ = '''
2649 Load %s information from given source(s)
2651 The ``source`` can be specified as the name of a %s file, the name of a
2652 directory containing %s files, a glob pattern of %s files, an open
2653 filehandle or an iterator yielding any of the forementioned sources.
2655 This function behaves as a generator yielding %s objects.
2656 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects)
2658 for f in iload_filename, iload_dirname, iload_glob, iload:
2659 f.__module__ = iload_fh.__module__
2661 return iload_filename, iload_dirname, iload_glob, iload
2664class Inconsistency(Exception):
2665 pass
2668def consistency_check(list_of_tuples, message='values differ:'):
2669 '''
2670 Check for inconsistencies.
2672 Given a list of tuples, check that all tuple elements except for first one
2673 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be
2674 valid because the coordinates at the two channels are the same.
2675 '''
2677 if len(list_of_tuples) >= 2:
2678 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]):
2679 raise Inconsistency('%s\n' % message + '\n'.join(
2680 ' %s: %s' % (t[0], ', '.join(str(x) for x in t[1:]))
2681 for t in list_of_tuples))
2684class defaultzerodict(dict):
2685 def __missing__(self, k):
2686 return 0
2689def mostfrequent(x):
2690 c = defaultzerodict()
2691 for e in x:
2692 c[e] += 1
2694 return sorted(list(c.keys()), key=lambda k: c[k])[-1]
2697def consistency_merge(list_of_tuples,
2698 message='values differ:',
2699 error='raise',
2700 merge=mostfrequent):
2702 assert error in ('raise', 'warn', 'ignore')
2704 if len(list_of_tuples) == 0:
2705 raise Exception('cannot merge empty sequence')
2707 try:
2708 consistency_check(list_of_tuples, message)
2709 return list_of_tuples[0][1:]
2710 except Inconsistency as e:
2711 if error == 'raise':
2712 raise
2714 elif error == 'warn':
2715 logger.warning(str(e))
2717 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]])
2720def short_to_list(nmax, it):
2721 import itertools
2723 if isinstance(it, list):
2724 return it
2726 li = []
2727 for i in range(nmax+1):
2728 try:
2729 li.append(next(it))
2730 except StopIteration:
2731 return li
2733 return itertools.chain(li, it)
2736def parse_md(f):
2737 try:
2738 with open(op.join(
2739 op.dirname(op.abspath(f)),
2740 'README.md'), 'r') as readme:
2741 mdstr = readme.read()
2742 except IOError as e:
2743 return 'Failed to get README.md: %s' % e
2745 # Remve the title
2746 mdstr = re.sub(r'^# .*\n?', '', mdstr)
2747 # Append sphinx reference to `pyrocko.` modules
2748 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr)
2749 # Convert Subsections to toc-less rubrics
2750 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr)
2751 return mdstr
2754def mpl_show(plt):
2755 import matplotlib
2756 if matplotlib.get_backend().lower() == 'agg':
2757 logger.warning('Cannot show() when using matplotlib "agg" backend')
2758 else:
2759 plt.show()
2762g_re_qsplit = re.compile(
2763 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)')
2764g_re_qsplit_sep = {}
2767def get_re_qsplit(sep):
2768 if sep is None:
2769 return g_re_qsplit
2770 else:
2771 if sep not in g_re_qsplit_sep:
2772 assert len(sep) == 1
2773 assert sep not in '\'"'
2774 esep = re.escape(sep)
2775 g_re_qsplit_sep[sep] = re.compile(
2776 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|'
2777 + r'([^' + esep + r']+|(?<=' + esep + r')(?=' + esep + r')|^(?=' + esep + r')|(?<=' + esep + r')$)') # noqa
2778 return g_re_qsplit_sep[sep]
2781g_re_trivial = re.compile(r'\A[^\'"\s]+\Z')
2782g_re_trivial_sep = {}
2785def get_re_trivial(sep):
2786 if sep is None:
2787 return g_re_trivial
2788 else:
2789 if sep not in g_re_qsplit_sep:
2790 assert len(sep) == 1
2791 assert sep not in '\'"'
2792 esep = re.escape(sep)
2793 g_re_trivial_sep[sep] = re.compile(r'\A[^\'"' + esep + r']+\Z')
2795 return g_re_trivial_sep[sep]
2798g_re_escape_s = re.compile(r'([\\\'])')
2799g_re_unescape_s = re.compile(r'\\([\\\'])')
2800g_re_escape_d = re.compile(r'([\\"])')
2801g_re_unescape_d = re.compile(r'\\([\\"])')
2804def escape_s(s):
2805 '''
2806 Backslash-escape single-quotes and backslashes.
2808 Example: ``Jack's`` => ``Jack\\'s``
2810 '''
2811 return g_re_escape_s.sub(r'\\\1', s)
2814def unescape_s(s):
2815 '''
2816 Unescape backslash-escaped single-quotes and backslashes.
2818 Example: ``Jack\\'s`` => ``Jack's``
2819 '''
2820 return g_re_unescape_s.sub(r'\1', s)
2823def escape_d(s):
2824 '''
2825 Backslash-escape double-quotes and backslashes.
2827 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"``
2828 '''
2829 return g_re_escape_d.sub(r'\\\1', s)
2832def unescape_d(s):
2833 '''
2834 Unescape backslash-escaped double-quotes and backslashes.
2836 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"``
2837 '''
2838 return g_re_unescape_d.sub(r'\1', s)
2841def qjoin_s(it, sep=None):
2842 '''
2843 Join sequence of strings into a line, single-quoting non-trivial strings.
2845 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"``
2846 '''
2847 re_trivial = get_re_trivial(sep)
2849 if sep is None:
2850 sep = ' '
2852 return sep.join(
2853 w if re_trivial.search(w) else "'%s'" % escape_s(w) for w in it)
2856def qjoin_d(it, sep=None):
2857 '''
2858 Join sequence of strings into a line, double-quoting non-trivial strings.
2860 Example: ``['55', 'Pete "The Robot" Smith']`` =>
2861 ``'55' "Pete \\\\"The Robot\\\\" Smith"'``
2862 '''
2863 re_trivial = get_re_trivial(sep)
2864 if sep is None:
2865 sep = ' '
2867 return sep.join(
2868 w if re_trivial.search(w) else '"%s"' % escape_d(w) for w in it)
2871def qsplit(s, sep=None):
2872 '''
2873 Split line into list of strings, allowing for quoted strings.
2875 Example: ``"55 'Sparrow\\\\'s Island'"`` =>
2876 ``["55", "Sparrow's Island"]``,
2877 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` =>
2878 ``['55', 'Pete "The Robot" Smith']``
2879 '''
2880 re_qsplit = get_re_qsplit(sep)
2881 return [
2882 (unescape_d(x[0]) or unescape_s(x[1]) or x[2])
2883 for x in re_qsplit.findall(s)]
2886g_have_warned_threadpoolctl = False
2889class threadpool_limits_dummy(object):
2891 def __init__(self, *args, **kwargs):
2892 pass
2894 def __enter__(self):
2895 global g_have_warned_threadpoolctl
2897 if not g_have_warned_threadpoolctl:
2898 logger.warning(
2899 'Cannot control number of BLAS threads because '
2900 '`threadpoolctl` module is not available. You may want to '
2901 'install `threadpoolctl`.')
2903 g_have_warned_threadpoolctl = True
2905 return self
2907 def __exit__(self, type, value, traceback):
2908 pass
2911def get_threadpool_limits():
2912 '''
2913 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail.
2914 '''
2916 try:
2917 from threadpoolctl import threadpool_limits
2918 return threadpool_limits
2920 except ImportError:
2921 return threadpool_limits_dummy
2924def fmt_summary(entries, widths):
2925 return ' | '.join(
2926 entry.ljust(width) for (entry, width) in zip(entries, widths))
2929def smart_weakref(obj, callback=None):
2930 if inspect.ismethod(obj):
2931 return weakref.WeakMethod(obj, callback)
2932 else:
2933 return weakref.ref(obj, callback)