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 float (:py:class:`numpy.float128`
15or :py:class:`numpy.float96`, whichever is available), aliased here as
16:py:class:`pyrocko.util.hpfloat`. High precision time stamps are required when
17handling data with sub-millisecond precision, i.e. kHz/MHz data streams and
18event catalogs derived from such data.
20Not all functions in Pyrocko and in programs depending on Pyrocko may work
21correctly with high precision times. Therefore, Pyrocko's high precision time
22handling mode has to be actively activated by user config, command line option
23or enforced within a certain script/program.
25The default high precision time handling mode can be configured globally with
26the user config variable
27:py:gattr:`~pyrocko.config.Config.use_high_precision_time`. Calling the
28function :py:func:`use_high_precision_time` overrides the default from the
29config file. This function may be called at startup of a program/script
30requiring a specific time handling mode.
32To create a valid time stamp for use in Pyrocko (e.g. in
33:py:class:`~pyrocko.model.Event` or :py:class:`~pyrocko.trace.Trace` objects),
34use:
36.. code-block :: python
38 import time
39 from pyrocko import util
41 # By default using mode selected in user config, override with:
42 # util.use_high_precision_time(True) # force high precision mode
43 # util.use_high_precision_time(False) # force low precision mode
45 t1 = util.str_to_time('2020-08-27 10:22:00')
46 t2 = util.str_to_time('2020-08-27 10:22:00.111222')
47 t3 = util.to_time_float(time.time())
49 # To get the appropriate float class, use:
51 time_float = util.get_time_float()
52 # -> float, numpy.float128 or numpy.float96
53 [isinstance(t, time_float) for t in [t1, t2, t3]]
54 # -> [True, True, True]
56 # Shortcut:
57 util.check_time_class(t1)
59Module content
60..............
61'''
63import time
64import logging
65import os
66import sys
67import re
68import calendar
69import math
70import fnmatch
71import inspect
72import weakref
73try:
74 import fcntl
75except ImportError:
76 fcntl = None
77import optparse
78import os.path as op
79import errno
81import numpy as num
82from scipy import signal
83import pyrocko
84from pyrocko import dummy_progressbar
87from urllib.parse import urlencode, quote, unquote # noqa
88from urllib.request import Request, build_opener, HTTPDigestAuthHandler # noqa
89from urllib.request import urlopen as _urlopen # noqa
90from urllib.error import HTTPError, URLError # noqa
93try:
94 import certifi
95 import ssl
96 g_ssl_context = ssl.create_default_context(cafile=certifi.where())
97except ImportError:
98 g_ssl_context = None
101class URLErrorSSL(URLError):
103 def __init__(self, urlerror):
104 self.__dict__ = urlerror.__dict__.copy()
106 def __str__(self):
107 return (
108 'Requesting web resource failed and the problem could be '
109 'related to SSL. Python standard libraries on some older '
110 'systems (like Ubuntu 14.04) are known to have trouble '
111 "with some SSL setups of today's servers: %s"
112 % URLError.__str__(self))
115def urlopen_ssl_check(*args, **kwargs):
116 try:
117 return urlopen(*args, **kwargs)
118 except URLError as e:
119 if str(e).find('SSL') != -1:
120 raise URLErrorSSL(e)
121 else:
122 raise
125def urlopen(*args, **kwargs):
127 if 'context' not in kwargs and g_ssl_context is not None:
128 kwargs['context'] = g_ssl_context
130 return _urlopen(*args, **kwargs)
133try:
134 long
135except NameError:
136 long = int
139force_dummy_progressbar = False
142try:
143 from pyrocko import util_ext
144except ImportError:
145 util_ext = None
148logger = logging.getLogger('pyrocko.util')
151try:
152 num_full = num.full
153except AttributeError:
154 def num_full(shape, fill_value, dtype=None, order='C'):
155 a = num.empty(shape, dtype=dtype, order=order)
156 a.fill(fill_value)
157 return a
159try:
160 num_full_like = num.full_like
161except AttributeError:
162 def num_full_like(arr, fill_value, dtype=None, order='K', subok=True):
163 a = num.empty_like(arr, dtype=dtype, order=order, subok=subok)
164 a.fill(fill_value)
165 return a
168g_setup_logging_args = 'pyrocko', 'warning'
171def setup_logging(programname='pyrocko', levelname='warning'):
172 '''
173 Initialize logging.
175 :param programname: program name to be written in log
176 :param levelname: string indicating the logging level ('debug', 'info',
177 'warning', 'error', 'critical')
179 This function is called at startup by most pyrocko programs to set up a
180 consistent logging format. This is simply a shortcut to a call to
181 :py:func:`logging.basicConfig()`.
182 '''
184 global g_setup_logging_args
185 g_setup_logging_args = (programname, levelname)
187 levels = {'debug': logging.DEBUG,
188 'info': logging.INFO,
189 'warning': logging.WARNING,
190 'error': logging.ERROR,
191 'critical': logging.CRITICAL}
193 logging.basicConfig(
194 level=levels[levelname],
195 format=programname+':%(name)-25s - %(levelname)-8s - %(message)s')
198def subprocess_setup_logging_args():
199 '''
200 Get arguments from previous call to setup_logging.
202 These can be sent down to a worker process so it can setup its logging
203 in the same way as the main process.
204 '''
205 return g_setup_logging_args
208def data_file(fn):
209 return os.path.join(os.path.split(__file__)[0], 'data', fn)
212class DownloadError(Exception):
213 pass
216class PathExists(DownloadError):
217 pass
220def _download(url, fpath, username=None, password=None,
221 force=False, method='download', stats=None,
222 status_callback=None, entries_wanted=None,
223 recursive=False, header=None):
225 import requests
226 from requests.auth import HTTPBasicAuth
227 from requests.exceptions import HTTPError as req_HTTPError
229 requests.adapters.DEFAULT_RETRIES = 5
230 urljoin = requests.compat.urljoin
232 session = requests.Session()
233 session.header = header
234 session.auth = None if username is None\
235 else HTTPBasicAuth(username, password)
237 status = {
238 'ntotal_files': 0,
239 'nread_files': 0,
240 'ntotal_bytes_all_files': 0,
241 'nread_bytes_all_files': 0,
242 'ntotal_bytes_current_file': 0,
243 'nread_bytes_current_file': 0,
244 'finished': False
245 }
247 try:
248 url_to_size = {}
250 if callable(status_callback):
251 status_callback(status)
253 if not recursive and url.endswith('/'):
254 raise DownloadError(
255 'URL: %s appears to be a directory'
256 ' but recurvise download is False' % url)
258 if recursive and not url.endswith('/'):
259 url += '/'
261 def parse_directory_tree(url, subdir=''):
262 r = session.get(urljoin(url, subdir))
263 r.raise_for_status()
265 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text)
267 files = sorted(set(subdir + fn for fn in entries
268 if not fn.endswith('/')))
270 if entries_wanted is not None:
271 files = [fn for fn in files
272 if (fn in wanted for wanted in entries_wanted)]
274 status['ntotal_files'] += len(files)
276 dirs = sorted(set(subdir + dn for dn in entries
277 if dn.endswith('/')
278 and dn not in ('./', '../')))
280 for dn in dirs:
281 files.extend(parse_directory_tree(
282 url, subdir=dn))
284 return files
286 def get_content_length(url):
287 if url not in url_to_size:
288 r = session.head(url, headers={'Accept-Encoding': ''})
290 content_length = r.headers.get('content-length', None)
291 if content_length is None:
292 logger.debug('Could not get HTTP header '
293 'Content-Length for %s' % url)
295 content_length = None
297 else:
298 content_length = int(content_length)
299 status['ntotal_bytes_all_files'] += content_length
300 if callable(status_callback):
301 status_callback(status)
303 url_to_size[url] = content_length
305 return url_to_size[url]
307 def download_file(url, fn):
308 logger.info('starting download of %s...' % url)
309 ensuredirs(fn)
311 fsize = get_content_length(url)
312 status['ntotal_bytes_current_file'] = fsize
313 status['nread_bytes_current_file'] = 0
314 if callable(status_callback):
315 status_callback(status)
317 r = session.get(url, stream=True, timeout=5)
318 r.raise_for_status()
320 frx = 0
321 fn_tmp = fn + '.%i.temp' % os.getpid()
322 with open(fn_tmp, 'wb') as f:
323 for d in r.iter_content(chunk_size=1024):
324 f.write(d)
325 frx += len(d)
327 status['nread_bytes_all_files'] += len(d)
328 status['nread_bytes_current_file'] += len(d)
329 if callable(status_callback):
330 status_callback(status)
332 os.rename(fn_tmp, fn)
334 if fsize is not None and frx != fsize:
335 logger.warning(
336 'HTTP header Content-Length: %i bytes does not match '
337 'download size %i bytes' % (fsize, frx))
339 logger.info('finished download of %s' % url)
341 status['nread_files'] += 1
342 if callable(status_callback):
343 status_callback(status)
345 if recursive:
346 if op.exists(fpath) and not force:
347 raise PathExists('path %s already exists' % fpath)
349 files = parse_directory_tree(url)
351 dsize = 0
352 for fn in files:
353 file_url = urljoin(url, fn)
354 dsize += get_content_length(file_url) or 0
356 if method == 'calcsize':
357 return dsize
359 else:
360 for fn in files:
361 file_url = urljoin(url, fn)
362 download_file(file_url, op.join(fpath, fn))
364 else:
365 status['ntotal_files'] += 1
366 if callable(status_callback):
367 status_callback(status)
369 fsize = get_content_length(url)
370 if method == 'calcsize':
371 return fsize
372 else:
373 download_file(url, fpath)
375 except req_HTTPError as e:
376 logging.warning('http error: %s' % e)
377 raise DownloadError('could not download file(s) from: %s' % url)
379 finally:
380 status['finished'] = True
381 if callable(status_callback):
382 status_callback(status)
383 session.close()
386def download_file(
387 url, fpath, username=None, password=None, status_callback=None,
388 **kwargs):
389 return _download(
390 url, fpath, username, password,
391 recursive=False,
392 status_callback=status_callback,
393 **kwargs)
396def download_dir(
397 url, fpath, username=None, password=None, status_callback=None,
398 **kwargs):
400 return _download(
401 url, fpath, username, password,
402 recursive=True,
403 status_callback=status_callback,
404 **kwargs)
407class HPFloatUnavailable(Exception):
408 pass
411class dummy_hpfloat(object):
412 def __init__(self, *args, **kwargs):
413 raise HPFloatUnavailable(
414 'NumPy lacks support for float128 or float96 data type on this '
415 'platform.')
418if hasattr(num, 'float128'):
419 hpfloat = num.float128
421elif hasattr(num, 'float96'):
422 hpfloat = num.float96
424else:
425 hpfloat = dummy_hpfloat
428g_time_float = None
429g_time_dtype = None
432class TimeFloatSettingError(Exception):
433 pass
436def use_high_precision_time(enabled):
437 '''
438 Globally force a specific time handling mode.
440 See :ref:`High precision time handling mode <time-handling-mode>`.
442 :param enabled: enable/disable use of high precision time type
443 :type enabled: bool
445 This function should be called before handling/reading any time data.
446 It can only be called once.
448 Special attention is required when using multiprocessing on a platform
449 which does not use fork under the hood. In such cases, the desired setting
450 must be set also in the subprocess.
451 '''
452 _setup_high_precision_time_mode(enabled_app=enabled)
455def _setup_high_precision_time_mode(enabled_app=False):
456 global g_time_float
457 global g_time_dtype
459 if not (g_time_float is None and g_time_dtype is None):
460 raise TimeFloatSettingError(
461 'Cannot set time handling mode: too late, it has already been '
462 'fixed by an earlier call.')
464 from pyrocko import config
466 conf = config.config()
467 enabled_config = conf.use_high_precision_time
469 enabled_env = os.environ.get('PYROCKO_USE_HIGH_PRECISION_TIME', None)
470 if enabled_env is not None:
471 try:
472 enabled_env = int(enabled_env) == 1
473 except ValueError:
474 raise TimeFloatSettingError(
475 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME '
476 'should be set to 0 or 1.')
478 enabled = enabled_config
479 mode_from = 'config variable `use_high_precision_time`'
480 notify = enabled
482 if enabled_env is not None:
483 if enabled_env != enabled:
484 notify = True
485 enabled = enabled_env
486 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`'
488 if enabled_app is not None:
489 if enabled_app != enabled:
490 notify = True
491 enabled = enabled_app
492 mode_from = 'application override'
494 logger.debug('''
495Pyrocko high precision time mode selection (latter override earlier):
496 config: %s
497 env: %s
498 app: %s
499 -> enabled: %s'''.lstrip() % (
500 enabled_config, enabled_env, enabled_app, enabled))
502 if notify:
503 logger.info('Pyrocko high precision time mode %s by %s.' % (
504 'activated' if enabled else 'deactivated',
505 mode_from))
507 if enabled:
508 g_time_float = hpfloat
509 g_time_dtype = hpfloat
510 else:
511 g_time_float = float
512 g_time_dtype = num.float64
515def get_time_float():
516 '''
517 Get the effective float class for timestamps.
519 See :ref:`High precision time handling mode <time-handling-mode>`.
521 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the
522 current time handling mode
523 '''
524 global g_time_float
526 if g_time_float is None:
527 _setup_high_precision_time_mode()
529 return g_time_float
532def get_time_dtype():
533 '''
534 Get effective NumPy float class to handle timestamps.
536 See :ref:`High precision time handling mode <time-handling-mode>`.
537 '''
539 global g_time_dtype
541 if g_time_dtype is None:
542 _setup_high_precision_time_mode()
544 return g_time_dtype
547def to_time_float(t):
548 '''
549 Convert float to valid time stamp in the current time handling mode.
551 See :ref:`High precision time handling mode <time-handling-mode>`.
552 '''
553 return get_time_float()(t)
556class TimestampTypeError(ValueError):
557 pass
560def check_time_class(t, error='raise'):
561 '''
562 Type-check variable against current time handling mode.
564 See :ref:`High precision time handling mode <time-handling-mode>`.
565 '''
567 if t == 0.0:
568 return
570 if not isinstance(t, get_time_float()):
571 message = (
572 'Timestamp %g is of type %s but should be of type %s with '
573 "Pyrocko's currently selected time handling mode.\n\n"
574 'See https://pyrocko.org/docs/current/library/reference/util.html'
575 '#high-precision-time-handling-mode' % (
576 t, type(t), get_time_float()))
578 if error == 'raise':
579 raise TimestampTypeError(message)
580 elif error == 'warn':
581 logger.warning(message)
582 else:
583 assert False
586class Stopwatch(object):
587 '''
588 Simple stopwatch to measure elapsed wall clock time.
590 Usage::
592 s = Stopwatch()
593 time.sleep(1)
594 print s()
595 time.sleep(1)
596 print s()
597 '''
599 def __init__(self):
600 self.start = time.time()
602 def __call__(self):
603 return time.time() - self.start
606def wrap(text, line_length=80):
607 '''
608 Paragraph and list-aware wrapping of text.
609 '''
611 text = text.strip('\n')
612 at_lineend = re.compile(r' *\n')
613 at_para = re.compile(r'((^|(\n\s*)?\n)(\s+[*] )|\n\s*\n)')
615 paragraphs = at_para.split(text)[::5]
616 listindents = at_para.split(text)[4::5]
617 newlist = at_para.split(text)[3::5]
619 listindents[0:0] = [None]
620 listindents.append(True)
621 newlist.append(None)
623 det_indent = re.compile(r'^ *')
625 outlines = []
626 for ip, p in enumerate(paragraphs):
627 if not p:
628 continue
630 if listindents[ip] is None:
631 _indent = det_indent.findall(p)[0]
632 findent = _indent
633 else:
634 findent = listindents[ip]
635 _indent = ' ' * len(findent)
637 ll = line_length - len(_indent)
638 llf = ll
640 oldlines = [s.strip() for s in at_lineend.split(p.rstrip())]
641 p1 = ' '.join(oldlines)
642 possible = re.compile(r'(^.{1,%i}|.{1,%i})( |$)' % (llf, ll))
643 for imatch, match in enumerate(possible.finditer(p1)):
644 parout = match.group(1)
645 if imatch == 0:
646 outlines.append(findent + parout)
647 else:
648 outlines.append(_indent + parout)
650 if ip != len(paragraphs)-1 and (
651 listindents[ip] is None
652 or newlist[ip] is not None
653 or listindents[ip+1] is None):
655 outlines.append('')
657 return outlines
660def ewrap(lines, width=80, indent=''):
661 lines = list(lines)
662 if not lines:
663 return ''
664 fwidth = max(len(s) for s in lines)
665 nx = max(1, (80-len(indent)) // (fwidth+1))
666 i = 0
667 rows = []
668 while i < len(lines):
669 rows.append(indent + ' '.join(x.ljust(fwidth) for x in lines[i:i+nx]))
670 i += nx
672 return '\n'.join(rows)
675class BetterHelpFormatter(optparse.IndentedHelpFormatter):
677 def __init__(self, *args, **kwargs):
678 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs)
680 def format_option(self, option):
681 '''
682 From IndentedHelpFormatter but using a different wrap method.
683 '''
685 help_text_position = 4 + self.current_indent
686 help_text_width = self.width - help_text_position
688 opts = self.option_strings[option]
689 opts = '%*s%s' % (self.current_indent, '', opts)
690 if option.help:
691 help_text = self.expand_default(option)
693 if self.help_position + len(help_text) + 1 <= self.width:
694 lines = [
695 '',
696 '%-*s %s' % (self.help_position, opts, help_text),
697 '']
698 else:
699 lines = ['']
700 lines.append(opts)
701 lines.append('')
702 if option.help:
703 help_lines = wrap(help_text, help_text_width)
704 lines.extend(['%*s%s' % (help_text_position, '', line)
705 for line in help_lines])
706 lines.append('')
708 return '\n'.join(lines)
710 def format_description(self, description):
711 if not description:
712 return ''
714 if self.current_indent == 0:
715 lines = []
716 else:
717 lines = ['']
719 lines.extend(wrap(description, self.width - self.current_indent))
720 if self.current_indent == 0:
721 lines.append('\n')
723 return '\n'.join(
724 ['%*s%s' % (self.current_indent, '', line) for line in lines])
727class ProgressBar:
728 def __init__(self, label, n):
729 from pyrocko.progress import progress
730 self._context = progress.view()
731 self._context.__enter__()
732 self._task = progress.task(label, n)
734 def update(self, i):
735 self._task.update(i)
737 def finish(self):
738 self._task.done()
739 if self._context:
740 self._context.__exit__()
741 self._context = None
744def progressbar(label, maxval):
745 if force_dummy_progressbar:
746 return dummy_progressbar.ProgressBar(maxval=maxval).start()
748 return ProgressBar(label, maxval)
751def progress_beg(label):
752 '''
753 Notify user that an operation has started.
755 :param label: name of the operation
757 To be used in conjuction with :py:func:`progress_end`.
758 '''
760 sys.stderr.write(label)
761 sys.stderr.flush()
764def progress_end(label=''):
765 '''
766 Notify user that an operation has ended.
768 :param label: name of the operation
770 To be used in conjuction with :py:func:`progress_beg`.
771 '''
773 sys.stderr.write(' done. %s\n' % label)
774 sys.stderr.flush()
777class ArangeError(Exception):
778 pass
781def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'):
782 '''
783 Return evenly spaced numbers over a specified interval.
785 Like :py:func:`numpy.arange` but returning floating point numbers by
786 default and with defined behaviour when stepsize is inconsistent with
787 interval bounds. It is considered inconsistent if the difference between
788 the closest multiple of ``step`` and ``stop`` is larger than ``epsilon *
789 step``. Inconsistencies are handled according to the ``error`` parameter.
790 If it is set to ``'raise'`` an exception of type :py:exc:`ArangeError` is
791 raised. If it is set to ``'round'``, ``'floor'``, or ``'ceil'``, ``stop``
792 is silently changed to the closest, the next smaller, or next larger
793 multiple of ``step``, respectively.
794 '''
796 assert error in ('raise', 'round', 'floor', 'ceil')
798 start = dtype(start)
799 stop = dtype(stop)
800 step = dtype(step)
802 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round)
804 n = int(rnd((stop - start) / step)) + 1
805 stop_check = start + (n-1) * step
807 if error == 'raise' and abs(stop_check - stop) > step * epsilon:
808 raise ArangeError(
809 'inconsistent range specification: start=%g, stop=%g, step=%g'
810 % (start, stop, step))
812 x = num.arange(n, dtype=dtype)
813 x *= step
814 x += start
815 return x
818def polylinefit(x, y, n_or_xnodes):
819 '''
820 Fit piece-wise linear function to data.
822 :param x,y: arrays with coordinates of data
823 :param n_or_xnodes: int, number of segments or x coordinates of polyline
825 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of
826 polyline, root-mean-square error
827 '''
829 x = num.asarray(x)
830 y = num.asarray(y)
832 if isinstance(n_or_xnodes, int):
833 n = n_or_xnodes
834 xmin = x.min()
835 xmax = x.max()
836 xnodes = num.linspace(xmin, xmax, n+1)
837 else:
838 xnodes = num.asarray(n_or_xnodes)
839 n = xnodes.size - 1
841 assert len(x) == len(y)
842 assert n > 0
844 ndata = len(x)
845 a = num.zeros((ndata+(n-1), n*2))
846 for i in range(n):
847 xmin_block = xnodes[i]
848 xmax_block = xnodes[i+1]
849 if i == n-1: # don't loose last point
850 indices = num.where(
851 num.logical_and(xmin_block <= x, x <= xmax_block))[0]
852 else:
853 indices = num.where(
854 num.logical_and(xmin_block <= x, x < xmax_block))[0]
856 a[indices, i*2] = x[indices]
857 a[indices, i*2+1] = 1.0
859 w = float(ndata)*100.
860 if i < n-1:
861 a[ndata+i, i*2] = xmax_block*w
862 a[ndata+i, i*2+1] = 1.0*w
863 a[ndata+i, i*2+2] = -xmax_block*w
864 a[ndata+i, i*2+3] = -1.0*w
866 d = num.concatenate((y, num.zeros(n-1)))
867 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2))
869 ynodes = num.zeros(n+1)
870 ynodes[:n] = model[:, 0]*xnodes[:n] + model[:, 1]
871 ynodes[1:] += model[:, 0]*xnodes[1:] + model[:, 1]
872 ynodes[1:n] *= 0.5
874 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2))
876 return xnodes, ynodes, rms_error
879def plf_integrate_piecewise(x_edges, x, y):
880 '''
881 Calculate definite integral of piece-wise linear function on intervals.
883 Use trapezoidal rule to calculate definite integral of a piece-wise linear
884 function for a series of consecutive intervals. ``x_edges`` and ``x`` must
885 be sorted.
887 :param x_edges: array with edges of the intervals
888 :param x,y: arrays with coordinates of piece-wise linear function's
889 control points
890 '''
892 x_all = num.concatenate((x, x_edges))
893 ii = num.argsort(x_all)
894 y_edges = num.interp(x_edges, x, y)
895 y_all = num.concatenate((y, y_edges))
896 xs = x_all[ii]
897 ys = y_all[ii]
898 y_all[ii[1:]] = num.cumsum((xs[1:] - xs[:-1]) * 0.5 * (ys[1:] + ys[:-1]))
899 return num.diff(y_all[-len(y_edges):])
902def diff_fd_1d_4o(dt, data):
903 '''
904 Approximate first derivative of an array (forth order, central FD).
906 :param dt: sampling interval
907 :param data: NumPy array with data samples
909 :returns: NumPy array with same shape as input
911 Interior points are approximated to fourth order, edge points to first
912 order right- or left-sided respectively, points next to edge to second
913 order central.
914 '''
915 import scipy.signal
917 ddata = num.empty_like(data, dtype=float)
919 if data.size >= 5:
920 ddata[2:-2] = scipy.signal.lfilter(
921 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt)
923 if data.size >= 3:
924 ddata[1] = (data[2] - data[0]) / (2. * dt)
925 ddata[-2] = (data[-1] - data[-3]) / (2. * dt)
927 if data.size >= 2:
928 ddata[0] = (data[1] - data[0]) / dt
929 ddata[-1] = (data[-1] - data[-2]) / dt
931 if data.size == 1:
932 ddata[0] = 0.0
934 return ddata
937def diff_fd_1d_2o(dt, data):
938 '''
939 Approximate first derivative of an array (second order, central FD).
941 :param dt: sampling interval
942 :param data: NumPy array with data samples
944 :returns: NumPy array with same shape as input
946 Interior points are approximated to second order, edge points to first
947 order right- or left-sided respectively.
949 Uses :py:func:`numpy.gradient`.
950 '''
952 return num.gradient(data, dt)
955def diff_fd_2d_4o(dt, data):
956 '''
957 Approximate second derivative of an array (forth order, central FD).
959 :param dt: sampling interval
960 :param data: NumPy array with data samples
962 :returns: NumPy array with same shape as input
964 Interior points are approximated to fourth order, next-to-edge points to
965 second order, edge points repeated.
966 '''
967 import scipy.signal
969 ddata = num.empty_like(data, dtype=float)
971 if data.size >= 5:
972 ddata[2:-2] = scipy.signal.lfilter(
973 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2)
975 if data.size >= 3:
976 ddata[:2] = (data[2] - 2.0 * data[1] + data[0]) / dt**2
977 ddata[-2:] = (data[-1] - 2.0 * data[-2] + data[-3]) / dt**2
979 if data.size < 3:
980 ddata[:] = 0.0
982 return ddata
985def diff_fd_2d_2o(dt, data):
986 '''
987 Approximate second derivative of an array (second order, central FD).
989 :param dt: sampling interval
990 :param data: NumPy array with data samples
992 :returns: NumPy array with same shape as input
994 Interior points are approximated to second order, edge points repeated.
995 '''
996 import scipy.signal
998 ddata = num.empty_like(data, dtype=float)
1000 if data.size >= 3:
1001 ddata[1:-1] = scipy.signal.lfilter(
1002 [1., -2., 1.], [1.], data)[2:] / (dt**2)
1004 ddata[0] = ddata[1]
1005 ddata[-1] = ddata[-2]
1007 if data.size < 3:
1008 ddata[:] = 0.0
1010 return ddata
1013def diff_fd(n, order, dt, data):
1014 '''
1015 Approximate 1st or 2nd derivative of an array.
1017 :param n: 1 for first derivative, 2 for second
1018 :param order: order of the approximation 2 and 4 are supported
1019 :param dt: sampling interval
1020 :param data: NumPy array with data samples
1022 :returns: NumPy array with same shape as input
1024 This is a frontend to the functions :py:func:`diff_fd_1d_2o`,
1025 :py:func:`diff_fd_1d_4o`, :py:func:`diff_fd_2d_2o`, and
1026 :py:func:`diff_fd_2d_4o`.
1028 Raises :py:exc:`ValueError` for unsupported `n` or `order`.
1029 '''
1031 funcs = {
1032 1: {2: diff_fd_1d_2o, 4: diff_fd_1d_4o},
1033 2: {2: diff_fd_2d_2o, 4: diff_fd_2d_4o}}
1035 try:
1036 funcs_n = funcs[n]
1037 except KeyError:
1038 raise ValueError(
1039 'pyrocko.util.diff_fd: '
1040 'Only 1st and 2sd derivatives are supported.')
1042 try:
1043 func = funcs_n[order]
1044 except KeyError:
1045 raise ValueError(
1046 'pyrocko.util.diff_fd: '
1047 'Order %i is not supported for %s derivative. Supported: %s' % (
1048 order, ['', '1st', '2nd'][n],
1049 ', '.join('%i' % order for order in sorted(funcs_n.keys()))))
1051 return func(dt, data)
1054class GlobalVars(object):
1055 reuse_store = dict()
1056 decitab_nmax = 0
1057 decitab = {}
1058 decimate_fir_coeffs = {}
1059 decimate_fir_remez_coeffs = {}
1060 decimate_iir_coeffs = {}
1061 re_frac = None
1064def decimate_coeffs(q, n=None, ftype='iir'):
1066 q = int(q)
1068 if n is None:
1069 if ftype == 'fir':
1070 n = 30
1071 elif ftype == 'fir-remez':
1072 n = 45*q
1073 else:
1074 n = 8
1076 if ftype == 'fir':
1077 coeffs = GlobalVars.decimate_fir_coeffs
1078 if (n, 1./q) not in coeffs:
1079 coeffs[n, 1./q] = signal.firwin(n+1, .75/q, window='hamming')
1081 b = coeffs[n, 1./q]
1082 return b, [1.], n
1084 elif ftype == 'fir-remez':
1085 coeffs = GlobalVars.decimate_fir_remez_coeffs
1086 if (n, 1./q) not in coeffs:
1087 coeffs[n, 1./q] = signal.remez(
1088 n+1, (0., .75/q, 1./q, 1.),
1089 (1., 0.), fs=2, weight=(1, 50))
1090 b = coeffs[n, 1./q]
1091 return b, [1.], n
1093 else:
1094 coeffs = GlobalVars.decimate_iir_coeffs
1095 if (n, 0.05, 0.8/q) not in coeffs:
1096 coeffs[n, 0.05, 0.8/q] = signal.cheby1(n, 0.05, 0.8/q)
1098 b, a = coeffs[n, 0.05, 0.8/q]
1099 return b, a, n
1102def decimate(x, q, n=None, ftype='iir', zi=None, ioff=0):
1103 '''
1104 Downsample the signal x by an integer factor q, using an order n filter
1106 By default, an order 8 Chebyshev type I filter is used or a 30 point FIR
1107 filter with hamming window if ftype is 'fir'.
1109 :param x: the signal to be downsampled (1D :class:`numpy.ndarray`)
1110 :param q: the downsampling factor
1111 :param n: order of the filter (1 less than the length of the filter for a
1112 `fir` filter)
1113 :param ftype: type of the filter; can be `iir`, `fir` or `fir-remez`
1115 :returns: the downsampled signal (1D :class:`numpy.ndarray`)
1117 '''
1119 b, a, n = decimate_coeffs(q, n, ftype)
1121 if zi is None or zi is True:
1122 zi_ = num.zeros(max(len(a), len(b))-1, dtype=float)
1123 else:
1124 zi_ = zi
1126 y, zf = signal.lfilter(b, a, x, zi=zi_)
1128 if zi is not None:
1129 return y[n//2+ioff::q].copy(), zf
1130 else:
1131 return y[n//2+ioff::q].copy()
1134class UnavailableDecimation(Exception):
1135 '''
1136 Exception raised by :py:func:`decitab` for unavailable decimation factors.
1137 '''
1139 pass
1142def gcd(a, b, epsilon=1e-7):
1143 '''
1144 Greatest common divisor.
1145 '''
1147 while b > epsilon*a:
1148 a, b = b, a % b
1150 return a
1153def lcm(a, b):
1154 '''
1155 Least common multiple.
1156 '''
1158 return a*b // gcd(a, b)
1161def mk_decitab(nmax=100):
1162 '''
1163 Make table with decimation sequences.
1165 Decimation from one sampling rate to a lower one is achieved by a
1166 successive application of :py:func:`decimation` with small integer
1167 downsampling factors (because using large downampling factors can make the
1168 decimation unstable or slow). This function sets up a table with downsample
1169 sequences for factors up to ``nmax``.
1170 '''
1172 tab = GlobalVars.decitab
1173 for i in range(1, 10):
1174 for j in range(1, i+1):
1175 for k in range(1, j+1):
1176 for l_ in range(1, k+1):
1177 for m in range(1, l_+1):
1178 p = i*j*k*l_*m
1179 if p > nmax:
1180 break
1181 if p not in tab:
1182 tab[p] = (i, j, k, l_, m)
1183 if i*j*k*l_ > nmax:
1184 break
1185 if i*j*k > nmax:
1186 break
1187 if i*j > nmax:
1188 break
1189 if i > nmax:
1190 break
1192 GlobalVars.decitab_nmax = nmax
1195def zfmt(n):
1196 return '%%0%ii' % (int(math.log10(n - 1)) + 1)
1199def _year_to_time(year):
1200 tt = (year, 1, 1, 0, 0, 0)
1201 return to_time_float(calendar.timegm(tt))
1204def _working_year(year):
1205 try:
1206 tt = (year, 1, 1, 0, 0, 0)
1207 t = calendar.timegm(tt)
1208 tt2_ = time.gmtime(t)
1209 tt2 = tuple(tt2_)[:6]
1210 if tt != tt2:
1211 return False
1213 s = '%i-01-01 00:00:00' % year
1214 s2 = time.strftime('%Y-%m-%d %H:%M:%S', tt2_)
1215 if s != s2:
1216 return False
1218 t3 = str_to_time(s2, format='%Y-%m-%d %H:%M:%S')
1219 s3 = time_to_str(t3, format='%Y-%m-%d %H:%M:%S')
1220 if s3 != s2:
1221 return False
1223 except Exception:
1224 return False
1226 return True
1229def working_system_time_range(year_min_lim=None, year_max_lim=None):
1230 '''
1231 Check time range supported by the systems's time conversion functions.
1233 Returns system time stamps of start of year of first/last fully supported
1234 year span. If this is before 1900 or after 2100, return first/last century
1235 which is fully supported.
1237 :returns: ``(tmin, tmax, year_min, year_max)``
1238 '''
1240 year0 = 2000
1241 year_min = year0
1242 year_max = year0
1244 itests = list(range(101))
1245 for i in range(19):
1246 itests.append(200 + i*100)
1248 for i in itests:
1249 year = year0 - i
1250 if year_min_lim is not None and year < year_min_lim:
1251 year_min = year_min_lim
1252 break
1253 elif not _working_year(year):
1254 break
1255 else:
1256 year_min = year
1258 for i in itests:
1259 year = year0 + i
1260 if year_max_lim is not None and year > year_max_lim:
1261 year_max = year_max_lim
1262 break
1263 elif not _working_year(year + 1):
1264 break
1265 else:
1266 year_max = year
1268 return (
1269 _year_to_time(year_min),
1270 _year_to_time(year_max),
1271 year_min, year_max)
1274g_working_system_time_range = None
1277def get_working_system_time_range():
1278 '''
1279 Caching variant of :py:func:`working_system_time_range`.
1280 '''
1282 global g_working_system_time_range
1283 if g_working_system_time_range is None:
1284 g_working_system_time_range = working_system_time_range()
1286 return g_working_system_time_range
1289def is_working_time(t):
1290 tmin, tmax, _, _ = get_working_system_time_range()
1291 return tmin <= t <= tmax
1294def julian_day_of_year(timestamp):
1295 '''
1296 Get the day number after the 1st of January of year in ``timestamp``.
1298 :returns: day number as int
1299 '''
1301 return time.gmtime(int(timestamp)).tm_yday
1304def hour_start(timestamp):
1305 '''
1306 Get beginning of hour for any point in time.
1308 :param timestamp: time instant as system timestamp (in seconds)
1310 :returns: instant of hour start as system timestamp
1311 '''
1313 tt = time.gmtime(int(timestamp))
1314 tts = tt[0:4] + (0, 0)
1315 return to_time_float(calendar.timegm(tts))
1318def day_start(timestamp):
1319 '''
1320 Get beginning of day for any point in time.
1322 :param timestamp: time instant as system timestamp (in seconds)
1324 :returns: instant of day start as system timestamp
1325 '''
1327 tt = time.gmtime(int(timestamp))
1328 tts = tt[0:3] + (0, 0, 0)
1329 return to_time_float(calendar.timegm(tts))
1332def month_start(timestamp):
1333 '''
1334 Get beginning of month for any point in time.
1336 :param timestamp: time instant as system timestamp (in seconds)
1338 :returns: instant of month start as system timestamp
1339 '''
1341 tt = time.gmtime(int(timestamp))
1342 tts = tt[0:2] + (1, 0, 0, 0)
1343 return to_time_float(calendar.timegm(tts))
1346def year_start(timestamp):
1347 '''
1348 Get beginning of year for any point in time.
1350 :param timestamp: time instant as system timestamp (in seconds)
1352 :returns: instant of year start as system timestamp
1353 '''
1355 tt = time.gmtime(int(timestamp))
1356 tts = tt[0:1] + (1, 1, 0, 0, 0)
1357 return to_time_float(calendar.timegm(tts))
1360def iter_days(tmin, tmax):
1361 '''
1362 Yields begin and end of days until given time span is covered.
1364 :param tmin,tmax: input time span
1366 :yields: tuples with (begin, end) of days as system timestamps
1367 '''
1369 t = day_start(tmin)
1370 while t < tmax:
1371 tend = day_start(t + 26*60*60)
1372 yield t, tend
1373 t = tend
1376def iter_months(tmin, tmax):
1377 '''
1378 Yields begin and end of months until given time span is covered.
1380 :param tmin,tmax: input time span
1382 :yields: tuples with (begin, end) of months as system timestamps
1383 '''
1385 t = month_start(tmin)
1386 while t < tmax:
1387 tend = month_start(t + 24*60*60*33)
1388 yield t, tend
1389 t = tend
1392def iter_years(tmin, tmax):
1393 '''
1394 Yields begin and end of years until given time span is covered.
1396 :param tmin,tmax: input time span
1398 :yields: tuples with (begin, end) of years as system timestamps
1399 '''
1401 t = year_start(tmin)
1402 while t < tmax:
1403 tend = year_start(t + 24*60*60*369)
1404 yield t, tend
1405 t = tend
1408def today():
1409 return day_start(time.time())
1412def tomorrow():
1413 return day_start(time.time() + 24*60*60)
1416def decitab(n):
1417 '''
1418 Get integer decimation sequence for given downampling factor.
1420 :param n: target decimation factor
1422 :returns: tuple with downsampling sequence
1423 '''
1425 if n > GlobalVars.decitab_nmax:
1426 mk_decitab(n*2)
1427 if n not in GlobalVars.decitab:
1428 raise UnavailableDecimation('ratio = %g' % n)
1429 return GlobalVars.decitab[n]
1432def ctimegm(s, format='%Y-%m-%d %H:%M:%S'):
1433 '''
1434 Convert string representing UTC time to system time.
1436 :param s: string to be interpreted
1437 :param format: format string passed to :py:func:`strptime`
1439 :returns: system time stamp
1441 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime.
1443 .. note::
1444 This function is to be replaced by :py:func:`str_to_time`.
1445 '''
1447 return calendar.timegm(time.strptime(s, format))
1450def gmctime(t, format='%Y-%m-%d %H:%M:%S'):
1451 '''
1452 Get string representation from system time, UTC.
1454 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime.
1456 .. note::
1457 This function is to be repaced by :py:func:`time_to_str`.
1458 '''
1460 return time.strftime(format, time.gmtime(t))
1463def gmctime_v(t, format='%a, %d %b %Y %H:%M:%S'):
1464 '''
1465 Get string representation from system time, UTC. Same as
1466 :py:func:`gmctime` but with a more verbose default format.
1468 .. note::
1469 This function is to be replaced by :py:func:`time_to_str`.
1470 '''
1472 return time.strftime(format, time.gmtime(t))
1475def gmctime_fn(t, format='%Y-%m-%d_%H-%M-%S'):
1476 '''
1477 Get string representation from system time, UTC. Same as
1478 :py:func:`gmctime` but with a default usable in filenames.
1480 .. note::
1481 This function is to be replaced by :py:func:`time_to_str`.
1482 '''
1484 return time.strftime(format, time.gmtime(t))
1487class TimeStrError(Exception):
1488 pass
1491class FractionalSecondsMissing(TimeStrError):
1492 '''
1493 Exception raised by :py:func:`str_to_time` when the given string lacks
1494 fractional seconds.
1495 '''
1497 pass
1500class FractionalSecondsWrongNumberOfDigits(TimeStrError):
1501 '''
1502 Exception raised by :py:func:`str_to_time` when the given string has an
1503 incorrect number of digits in the fractional seconds part.
1504 '''
1506 pass
1509def _endswith_n(s, endings):
1510 for ix, x in enumerate(endings):
1511 if s.endswith(x):
1512 return ix
1513 return -1
1516def str_to_time(s, format='%Y-%m-%d %H:%M:%S.OPTFRAC'):
1517 '''
1518 Convert string representing UTC time to floating point system time.
1520 :param s: string representing UTC time
1521 :param format: time string format
1522 :returns: system time stamp as floating point value
1524 Uses the semantics of :py:func:`time.strptime` but allows for fractional
1525 seconds. If the format ends with ``'.FRAC'``, anything after a dot is
1526 interpreted as fractional seconds. If the format ends with ``'.OPTFRAC'``,
1527 the fractional part, including the dot is made optional. The latter has the
1528 consequence, that the time strings and the format may not contain any other
1529 dots. If the format ends with ``'.xFRAC'`` where x is 1, 2, or 3, it is
1530 ensured, that exactly that number of digits are present in the fractional
1531 seconds.
1532 '''
1534 time_float = get_time_float()
1536 if util_ext is not None:
1537 try:
1538 t, tfrac = util_ext.stt(s, format)
1539 except util_ext.UtilExtError as e:
1540 raise TimeStrError(
1541 '%s, string=%s, format=%s' % (str(e), s, format))
1543 return time_float(t)+tfrac
1545 fracsec = 0.
1546 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC'
1548 iend = _endswith_n(format, fixed_endings)
1549 if iend != -1:
1550 dotpos = s.rfind('.')
1551 if dotpos == -1:
1552 raise FractionalSecondsMissing(
1553 'string=%s, format=%s' % (s, format))
1555 if iend > 0 and iend != (len(s)-dotpos-1):
1556 raise FractionalSecondsWrongNumberOfDigits(
1557 'string=%s, format=%s' % (s, format))
1559 format = format[:-len(fixed_endings[iend])]
1560 fracsec = float(s[dotpos:])
1561 s = s[:dotpos]
1563 elif format.endswith('.OPTFRAC'):
1564 dotpos = s.rfind('.')
1565 format = format[:-8]
1566 if dotpos != -1 and len(s[dotpos:]) > 1:
1567 fracsec = float(s[dotpos:])
1569 if dotpos != -1:
1570 s = s[:dotpos]
1572 try:
1573 return time_float(calendar.timegm(time.strptime(s, format))) \
1574 + fracsec
1575 except ValueError as e:
1576 raise TimeStrError('%s, string=%s, format=%s' % (str(e), s, format))
1579stt = str_to_time
1582def str_to_time_fillup(s):
1583 '''
1584 Default :py:func:`str_to_time` with filling in of missing values.
1586 Allows e.g. `'2010-01-01 00:00:00'` as `'2010-01-01 00:00'`,
1587 `'2010-01-01 00'`, ..., or `'2010'`.
1588 '''
1590 if s == 'now':
1591 return time.time()
1593 if len(s) in (4, 7, 10, 13, 16):
1594 s += '0000-01-01 00:00:00'[len(s):]
1596 return str_to_time(s)
1599def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'):
1600 '''
1601 Get string representation for floating point system time.
1603 :param t: floating point system time
1604 :param format: time string format
1605 :returns: string representing UTC time
1607 Uses the semantics of :py:func:`time.strftime` but additionally allows for
1608 fractional seconds. If ``format`` contains ``'.xFRAC'``, where ``x`` is a
1609 digit between 1 and 9, this is replaced with the fractional part of ``t``
1610 with ``x`` digits precision.
1611 '''
1613 if pyrocko.grumpy > 0:
1614 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise')
1616 if isinstance(format, int):
1617 if format > 0:
1618 format = '%Y-%m-%d %H:%M:%S.' + '%iFRAC' % format
1619 else:
1620 format = '%Y-%m-%d %H:%M:%S'
1622 if util_ext is not None:
1623 t0 = num.floor(t)
1624 try:
1625 return util_ext.tts(int(t0), float(t - t0), format)
1626 except util_ext.UtilExtError as e:
1627 raise TimeStrError(
1628 '%s, timestamp=%f, format=%s' % (str(e), t, format))
1630 if not GlobalVars.re_frac:
1631 GlobalVars.re_frac = re.compile(r'\.[1-9]FRAC')
1632 GlobalVars.frac_formats = dict(
1633 [('.%sFRAC' % x, '%.'+x+'f') for x in '123456789'])
1635 ts = float(num.floor(t))
1636 tfrac = t-ts
1638 m = GlobalVars.re_frac.search(format)
1639 if m:
1640 sfrac = (GlobalVars.frac_formats[m.group(0)] % tfrac)
1641 if sfrac[0] == '1':
1642 ts += 1.
1644 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1)
1646 return time.strftime(format, time.gmtime(ts))
1649tts = time_to_str
1650_abbr_weekday = 'Mon Tue Wed Thu Fri Sat Sun'.split()
1651_abbr_month = 'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split()
1654def mystrftime(fmt=None, tt=None, milliseconds=0):
1655 # Needed by Snuffler for the time axis. In other cases `time_to_str`
1656 # should be used.
1658 if fmt is None:
1659 fmt = '%Y-%m-%d %H:%M:%S .%r'
1661 # Get these two locale independent, needed by Snuffler.
1662 # Setting the locale seems to have no effect.
1663 fmt = fmt.replace('%a', _abbr_weekday[tt.tm_wday])
1664 fmt = fmt.replace('%b', _abbr_month[tt.tm_mon-1])
1666 fmt = fmt.replace('%r', '%03i' % int(round(milliseconds)))
1667 fmt = fmt.replace('%u', '%06i' % int(round(milliseconds*1000)))
1668 fmt = fmt.replace('%n', '%09i' % int(round(milliseconds*1000000)))
1669 return time.strftime(fmt, tt)
1672def gmtime_x(timestamp):
1673 etimestamp = float(num.floor(timestamp))
1674 tt = time.gmtime(etimestamp)
1675 ms = (timestamp-etimestamp)*1000
1676 return tt, ms
1679def plural_s(n):
1680 if not isinstance(n, int):
1681 n = len(n)
1683 return 's' if n != 1 else ''
1686def ensuredirs(dst):
1687 '''
1688 Create all intermediate path components for a target path.
1690 :param dst: target path
1692 The leaf part of the target path is not created (use :py:func:`ensuredir`
1693 if a the target path is a directory to be created).
1694 '''
1696 d, x = os.path.split(dst.rstrip(os.sep))
1697 dirs = []
1698 while d and not os.path.exists(d):
1699 dirs.append(d)
1700 d, x = os.path.split(d)
1702 dirs.reverse()
1704 for d in dirs:
1705 try:
1706 os.mkdir(d)
1707 except OSError as e:
1708 if not e.errno == errno.EEXIST:
1709 raise
1712def ensuredir(dst):
1713 '''
1714 Create directory and all intermediate path components to it as needed.
1716 :param dst: directory name
1718 Nothing is done if the given target already exists.
1719 '''
1721 if os.path.exists(dst):
1722 return
1724 dst.rstrip(os.sep)
1726 ensuredirs(dst)
1727 try:
1728 os.mkdir(dst)
1729 except OSError as e:
1730 if not e.errno == errno.EEXIST:
1731 raise
1734def reuse(x):
1735 '''
1736 Get unique instance of an object.
1738 :param x: hashable object
1739 :returns: reference to x or an equivalent object
1741 Cache object ``x`` in a global dict for reuse, or if x already
1742 is in that dict, return a reference to it.
1743 '''
1745 grs = GlobalVars.reuse_store
1746 if x not in grs:
1747 grs[x] = x
1748 return grs[x]
1751def deuse(x):
1752 grs = GlobalVars.reuse_store
1753 if x in grs:
1754 del grs[x]
1757class Anon(object):
1758 '''
1759 Dict-to-object utility.
1761 Any given arguments are stored as attributes.
1763 Example::
1765 a = Anon(x=1, y=2)
1766 print a.x, a.y
1767 '''
1769 def __init__(self, **dict):
1770 for k in dict:
1771 self.__dict__[k] = dict[k]
1774def iter_select_files(
1775 paths,
1776 include=None,
1777 exclude=None,
1778 selector=None,
1779 show_progress=True,
1780 pass_through=None):
1782 '''
1783 Recursively select files (generator variant).
1785 See :py:func:`select_files`.
1786 '''
1788 if show_progress:
1789 progress_beg('selecting files...')
1791 ngood = 0
1792 check_include = None
1793 if include is not None:
1794 rinclude = re.compile(include)
1796 def check_include(path):
1797 m = rinclude.search(path)
1798 if not m:
1799 return False
1801 if selector is None:
1802 return True
1804 infos = Anon(**m.groupdict())
1805 return selector(infos)
1807 check_exclude = None
1808 if exclude is not None:
1809 rexclude = re.compile(exclude)
1811 def check_exclude(path):
1812 return not bool(rexclude.search(path))
1814 if check_include and check_exclude:
1816 def check(path):
1817 return check_include(path) and check_exclude(path)
1819 elif check_include:
1820 check = check_include
1822 elif check_exclude:
1823 check = check_exclude
1825 else:
1826 check = None
1828 if isinstance(paths, str):
1829 paths = [paths]
1831 for path in paths:
1832 if pass_through and pass_through(path):
1833 if check is None or check(path):
1834 yield path
1836 elif os.path.isdir(path):
1837 for (dirpath, dirnames, filenames) in os.walk(path):
1838 dirnames.sort()
1839 filenames.sort()
1840 for filename in filenames:
1841 path = op.join(dirpath, filename)
1842 if check is None or check(path):
1843 yield os.path.abspath(path)
1844 ngood += 1
1845 else:
1846 if check is None or check(path):
1847 yield os.path.abspath(path)
1848 ngood += 1
1850 if show_progress:
1851 progress_end('%i file%s selected.' % (ngood, plural_s(ngood)))
1854def select_files(
1855 paths, include=None, exclude=None, selector=None, show_progress=True,
1856 regex=None):
1858 '''
1859 Recursively select files.
1861 :param paths: entry path names
1862 :param include: pattern for conditional inclusion
1863 :param exclude: pattern for conditional exclusion
1864 :param selector: callback for conditional inclusion
1865 :param show_progress: if True, indicate start and stop of processing
1866 :param regex: alias for ``include`` (backwards compatibility)
1867 :returns: list of path names
1869 Recursively finds all files under given entry points ``paths``. If
1870 parameter ``include`` is a regular expression, only files with matching
1871 path names are included. If additionally parameter ``selector`` is given a
1872 callback function, only files for which the callback returns ``True`` are
1873 included. The callback should take a single argument. The callback is
1874 called with a single argument, an object, having as attributes, any named
1875 groups given in ``include``.
1877 Examples
1879 To find all files ending in ``'.mseed'`` or ``'.msd'``::
1881 select_files(paths,
1882 include=r'\\.(mseed|msd)$')
1884 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for
1885 the year::
1887 select_files(paths,
1888 include=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$',
1889 selector=(lambda x: int(x.year) == 2009))
1890 '''
1891 if regex is not None:
1892 assert include is None
1893 include = regex
1895 return list(iter_select_files(
1896 paths, include, exclude, selector, show_progress))
1899def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
1900 '''
1901 Convert positive integer to a base36 string.
1902 '''
1904 if not isinstance(number, (int, long)):
1905 raise TypeError('number must be an integer')
1906 if number < 0:
1907 raise ValueError('number must be positive')
1909 # Special case for small numbers
1910 if number < 36:
1911 return alphabet[number]
1913 base36 = ''
1914 while number != 0:
1915 number, i = divmod(number, 36)
1916 base36 = alphabet[i] + base36
1918 return base36
1921def base36decode(number):
1922 '''
1923 Decode base36 endcoded positive integer.
1924 '''
1926 return int(number, 36)
1929class UnpackError(Exception):
1930 '''
1931 Exception raised when :py:func:`unpack_fixed` encounters an error.
1932 '''
1934 pass
1937ruler = ''.join(['%-10i' % i for i in range(8)]) \
1938 + '\n' + '0123456789' * 8 + '\n'
1941def unpack_fixed(format, line, *callargs):
1942 '''
1943 Unpack fixed format string, as produced by many fortran codes.
1945 :param format: format specification
1946 :param line: string to be processed
1947 :param callargs: callbacks for callback fields in the format
1949 The format is described by a string of comma-separated fields. Each field
1950 is defined by a character for the field type followed by the field width. A
1951 questionmark may be appended to the field description to allow the argument
1952 to be optional (The data string is then allowed to be filled with blanks
1953 and ``None`` is returned in this case).
1955 The following field types are available:
1957 ==== ================================================================
1958 Type Description
1959 ==== ================================================================
1960 A string (full field width is extracted)
1961 a string (whitespace at the beginning and the end is removed)
1962 i integer value
1963 f floating point value
1964 @ special type, a callback must be given for the conversion
1965 x special field type to skip parts of the string
1966 ==== ================================================================
1967 '''
1969 ipos = 0
1970 values = []
1971 icall = 0
1972 for form in format.split(','):
1973 form = form.strip()
1974 optional = form[-1] == '?'
1975 form = form.rstrip('?')
1976 typ = form[0]
1977 ln = int(form[1:])
1978 s = line[ipos:ipos+ln]
1979 cast = {
1980 'x': None,
1981 'A': str,
1982 'a': lambda x: x.strip(),
1983 'i': int,
1984 'f': float,
1985 '@': 'extra'}[typ]
1987 if cast == 'extra':
1988 cast = callargs[icall]
1989 icall += 1
1991 if cast is not None:
1992 if optional and s.strip() == '':
1993 values.append(None)
1994 else:
1995 try:
1996 values.append(cast(s))
1997 except Exception:
1998 mark = [' '] * 80
1999 mark[ipos:ipos+ln] = ['^'] * ln
2000 mark = ''.join(mark)
2001 raise UnpackError(
2002 'Invalid cast to type "%s" at position [%i:%i] of '
2003 'line: \n%s%s\n%s' % (
2004 typ, ipos, ipos+ln, ruler, line.rstrip(), mark))
2006 ipos += ln
2008 return values
2011_pattern_cache = {}
2014def _nslc_pattern(pattern):
2015 if pattern not in _pattern_cache:
2016 rpattern = re.compile(fnmatch.translate(pattern), re.I)
2017 _pattern_cache[pattern] = rpattern
2018 else:
2019 rpattern = _pattern_cache[pattern]
2021 return rpattern
2024def match_nslc(patterns, nslc):
2025 '''
2026 Match network-station-location-channel code against pattern or list of
2027 patterns.
2029 :param patterns: pattern or list of patterns
2030 :param nslc: tuple with (network, station, location, channel) as strings
2032 :returns: ``True`` if the pattern matches or if any of the given patterns
2033 match; or ``False``.
2035 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq].
2037 Example::
2039 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True
2040 '''
2042 if isinstance(patterns, str):
2043 patterns = [patterns]
2045 if not isinstance(nslc, str):
2046 s = '.'.join(nslc)
2047 else:
2048 s = nslc
2050 for pattern in patterns:
2051 if _nslc_pattern(pattern).match(s):
2052 return True
2054 return False
2057def match_nslcs(patterns, nslcs):
2058 '''
2059 Get network-station-location-channel codes that match given pattern or any
2060 of several given patterns.
2062 :param patterns: pattern or list of patterns
2063 :param nslcs: list of (network, station, location, channel) tuples
2065 See also :py:func:`match_nslc`
2066 '''
2068 matching = []
2069 for nslc in nslcs:
2070 if match_nslc(patterns, nslc):
2071 matching.append(nslc)
2073 return matching
2076class Timeout(Exception):
2077 pass
2080def create_lockfile(fn, timeout=None, timewarn=10.):
2081 t0 = time.time()
2083 while True:
2084 try:
2085 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL)
2086 os.close(f)
2087 return
2089 except OSError as e:
2090 if e.errno in (errno.EEXIST, 13): # 13 occurs on windows
2091 pass # retry
2092 else:
2093 raise
2095 tnow = time.time()
2097 if timeout is not None and tnow - t0 > timeout:
2098 raise Timeout(
2099 'Timeout (%gs) occured while waiting to get exclusive '
2100 'access to: %s' % (timeout, fn))
2102 if timewarn is not None and tnow - t0 > timewarn:
2103 logger.warning(
2104 'Waiting since %gs to get exclusive access to: %s' % (
2105 timewarn, fn))
2107 timewarn *= 2
2109 time.sleep(0.01)
2112def delete_lockfile(fn):
2113 os.unlink(fn)
2116class Lockfile(Exception):
2118 def __init__(self, path, timeout=5, timewarn=10.):
2119 self._path = path
2120 self._timeout = timeout
2121 self._timewarn = timewarn
2123 def __enter__(self):
2124 create_lockfile(
2125 self._path, timeout=self._timeout, timewarn=self._timewarn)
2126 return None
2128 def __exit__(self, type, value, traceback):
2129 delete_lockfile(self._path)
2132class SoleError(Exception):
2133 '''
2134 Exception raised by objects of type :py:class:`Sole`, when an concurrent
2135 instance is running.
2136 '''
2138 pass
2141class Sole(object):
2142 '''
2143 Use POSIX advisory file locking to ensure that only a single instance of a
2144 program is running.
2146 :param pid_path: path to lockfile to be used
2148 Usage::
2150 from pyrocko.util import Sole, SoleError, setup_logging
2151 import os
2153 setup_logging('my_program')
2155 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock')
2156 try:
2157 sole = Sole(pid_path)
2159 except SoleError, e:
2160 logger.fatal( str(e) )
2161 sys.exit(1)
2162 '''
2164 def __init__(self, pid_path):
2165 self._pid_path = pid_path
2166 self._other_running = False
2167 ensuredirs(self._pid_path)
2168 self._lockfile = None
2169 self._os = os
2171 if not fcntl:
2172 raise SoleError(
2173 'Python standard library module "fcntl" not available on '
2174 'this platform.')
2176 self._fcntl = fcntl
2178 try:
2179 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY)
2180 except Exception:
2181 raise SoleError(
2182 'Cannot open lockfile (path = %s)' % self._pid_path)
2184 try:
2185 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB)
2187 except IOError:
2188 self._other_running = True
2189 try:
2190 f = open(self._pid_path, 'r')
2191 pid = f.read().strip()
2192 f.close()
2193 except Exception:
2194 pid = '?'
2196 raise SoleError('Other instance is running (pid = %s)' % pid)
2198 try:
2199 os.ftruncate(self._lockfile, 0)
2200 os.write(self._lockfile, '%i\n' % os.getpid())
2201 os.fsync(self._lockfile)
2203 except Exception:
2204 # the pid is only stored for user information, so this is allowed
2205 # to fail
2206 pass
2208 def __del__(self):
2209 if not self._other_running:
2210 if self._lockfile is not None:
2211 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN)
2212 self._os.close(self._lockfile)
2213 try:
2214 self._os.unlink(self._pid_path)
2215 except Exception:
2216 pass
2219re_escapequotes = re.compile(r"(['\\])")
2222def escapequotes(s):
2223 return re_escapequotes.sub(r'\\\1', s)
2226class TableWriter(object):
2227 '''
2228 Write table of space separated values to a file.
2230 :param f: file like object
2232 Strings containing spaces are quoted on output.
2233 '''
2235 def __init__(self, f):
2236 self._f = f
2238 def writerow(self, row, minfieldwidths=None):
2240 '''
2241 Write one row of values to underlying file.
2243 :param row: iterable of values
2244 :param minfieldwidths: minimum field widths for the values
2246 Each value in in ``row`` is converted to a string and optionally padded
2247 with blanks. The resulting strings are output separated with blanks. If
2248 any values given are strings and if they contain whitespace, they are
2249 quoted with single quotes, and any internal single quotes are
2250 backslash-escaped.
2251 '''
2253 out = []
2255 for i, x in enumerate(row):
2256 w = 0
2257 if minfieldwidths and i < len(minfieldwidths):
2258 w = minfieldwidths[i]
2260 if isinstance(x, str):
2261 if re.search(r"\s|'", x):
2262 x = "'%s'" % escapequotes(x)
2264 x = x.ljust(w)
2265 else:
2266 x = str(x).rjust(w)
2268 out.append(x)
2270 self._f.write(' '.join(out).rstrip() + '\n')
2273class TableReader(object):
2275 '''
2276 Read table of space separated values from a file.
2278 :param f: file-like object
2280 This uses Pythons shlex module to tokenize lines. Should deal correctly
2281 with quoted strings.
2282 '''
2284 def __init__(self, f):
2285 self._f = f
2286 self.eof = False
2288 def readrow(self):
2289 '''
2290 Read one row from the underlying file, tokenize it with shlex.
2292 :returns: tokenized line as a list of strings.
2293 '''
2295 line = self._f.readline()
2296 if not line:
2297 self.eof = True
2298 return []
2299 line.strip()
2300 if line.startswith('#'):
2301 return []
2303 return qsplit(line)
2306def gform(number, significant_digits=3):
2307 '''
2308 Pretty print floating point numbers.
2310 Align floating point numbers at the decimal dot.
2312 ::
2314 | -d.dde+xxx|
2315 | -d.dde+xx |
2316 |-ddd. |
2317 | -dd.d |
2318 | -d.dd |
2319 | -0.ddd |
2320 | -0.0ddd |
2321 | -0.00ddd |
2322 | -d.dde-xx |
2323 | -d.dde-xxx|
2324 | nan|
2327 The formatted string has length ``significant_digits * 2 + 6``.
2328 '''
2330 no_exp_range = (pow(10., -1),
2331 pow(10., significant_digits))
2332 width = significant_digits+significant_digits-1+1+1+5
2334 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.:
2335 s = ('%#.*g' % (significant_digits, number)).rstrip('0')
2336 else:
2337 s = '%.*E' % (significant_digits-1, number)
2338 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width)
2339 if s.strip().lower() == 'nan':
2340 s = 'nan'.rjust(width)
2341 return s
2344def human_bytesize(value):
2346 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split()
2348 if value == 1:
2349 return '1 Byte'
2351 for i, ext in enumerate(exts):
2352 x = float(value) / 1000**i
2353 if round(x) < 10. and not value < 1000:
2354 return '%.1f %s' % (x, ext)
2355 if round(x) < 1000.:
2356 return '%.0f %s' % (x, ext)
2358 return '%i Bytes' % value
2361re_compatibility = re.compile(
2362 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' +
2363 r'(dummy|poel|qseis|qssp))\.'
2364)
2367def pf_is_old(fn):
2368 oldstyle = False
2369 with open(fn, 'r') as f:
2370 for line in f:
2371 if re_compatibility.search(line):
2372 oldstyle = True
2374 return oldstyle
2377def pf_upgrade(fn):
2378 need = pf_is_old(fn)
2379 if need:
2380 fn_temp = fn + '.temp'
2382 with open(fn, 'r') as fin:
2383 with open(fn_temp, 'w') as fout:
2384 for line in fin:
2385 line = re_compatibility.sub('!pf.', line)
2386 fout.write(line)
2388 os.rename(fn_temp, fn)
2390 return need
2393def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'):
2394 '''
2395 Extract leap second information from tzdata.
2397 Based on example at http://stackoverflow.com/questions/19332902/\
2398 extract-historic-leap-seconds-from-tzdata
2400 See also 'man 5 tzfile'.
2401 '''
2403 from struct import unpack, calcsize
2404 out = []
2405 with open(tzfile, 'rb') as f:
2406 # read header
2407 fmt = '>4s c 15x 6l'
2408 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt,
2409 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt)))
2410 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file'
2412 # skip over some uninteresting data
2413 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict(
2414 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt)
2415 f.read(calcsize(fmt))
2417 # read leap-seconds
2418 fmt = '>2l'
2419 for i in range(leapcnt):
2420 tleap, nleap = unpack(fmt, f.read(calcsize(fmt)))
2421 out.append((tleap-nleap+1, nleap))
2423 return out
2426class LeapSecondsError(Exception):
2427 pass
2430class LeapSecondsOutdated(LeapSecondsError):
2431 pass
2434class InvalidLeapSecondsFile(LeapSecondsOutdated):
2435 pass
2438def parse_leap_seconds_list(fn):
2439 data = []
2440 texpires = None
2441 try:
2442 t0 = int(round(str_to_time('1900-01-01 00:00:00')))
2443 except TimeStrError:
2444 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800
2446 tnow = int(round(time.time()))
2448 if not op.exists(fn):
2449 raise LeapSecondsOutdated('no leap seconds file found')
2451 try:
2452 with open(fn, 'rb') as f:
2453 for line in f:
2454 if line.strip().startswith(b'<!DOCTYPE'):
2455 raise InvalidLeapSecondsFile('invalid leap seconds file')
2457 if line.startswith(b'#@'):
2458 texpires = int(line.split()[1]) + t0
2459 elif line.startswith(b'#') or len(line) < 5:
2460 pass
2461 else:
2462 toks = line.split()
2463 t = int(toks[0]) + t0
2464 nleap = int(toks[1]) - 10
2465 data.append((t, nleap))
2467 except IOError:
2468 raise LeapSecondsError('cannot read leap seconds file %s' % fn)
2470 if texpires is None or tnow > texpires:
2471 raise LeapSecondsOutdated('leap seconds list is outdated')
2473 return data
2476def read_leap_seconds2():
2477 from pyrocko import config
2478 conf = config.config()
2479 fn = conf.leapseconds_path
2480 url = conf.leapseconds_url
2481 # check for outdated default URL
2482 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list':
2483 url = 'https://www.ietf.org/timezones/data/leap-seconds.list'
2484 logger.info(
2485 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2487 for i in range(3):
2488 try:
2489 return parse_leap_seconds_list(fn)
2491 except LeapSecondsOutdated:
2492 try:
2493 logger.info('updating leap seconds list...')
2494 download_file(url, fn)
2496 except Exception as e:
2497 raise LeapSecondsError(
2498 'cannot download leap seconds list from %s to %s (%s)'
2499 % (url, fn, e))
2501 raise LeapSecondsError('Could not retrieve/read leap seconds file.')
2504def gps_utc_offset(t_utc):
2505 '''
2506 Time offset t_gps - t_utc for a given t_utc.
2507 '''
2508 ls = read_leap_seconds2()
2509 i = 0
2510 if t_utc < ls[0][0]:
2511 return ls[0][1] - 1 - 9
2513 while i < len(ls) - 1:
2514 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]:
2515 return ls[i][1] - 9
2516 i += 1
2518 return ls[-1][1] - 9
2521def utc_gps_offset(t_gps):
2522 '''
2523 Time offset t_utc - t_gps for a given t_gps.
2524 '''
2525 ls = read_leap_seconds2()
2527 if t_gps < ls[0][0] + ls[0][1] - 9:
2528 return - (ls[0][1] - 1 - 9)
2530 i = 0
2531 while i < len(ls) - 1:
2532 if ls[i][0] + ls[i][1] - 9 <= t_gps \
2533 and t_gps < ls[i+1][0] + ls[i+1][1] - 9:
2534 return - (ls[i][1] - 9)
2535 i += 1
2537 return - (ls[-1][1] - 9)
2540def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'):
2541 import itertools
2542 import glob
2543 from pyrocko.io.io_common import FileLoadError
2545 def iload_filename(filename, **kwargs):
2546 try:
2547 with open(filename, 'rb') as f:
2548 for cr in iload_fh(f, **kwargs):
2549 yield cr
2551 except FileLoadError as e:
2552 e.set_context('filename', filename)
2553 raise
2555 def iload_dirname(dirname, **kwargs):
2556 for entry in os.listdir(dirname):
2557 fpath = op.join(dirname, entry)
2558 if op.isfile(fpath):
2559 for cr in iload_filename(fpath, **kwargs):
2560 yield cr
2562 def iload_glob(pattern, **kwargs):
2564 for fn in glob.iglob(pattern):
2565 for cr in iload_filename(fn, **kwargs):
2566 yield cr
2568 def iload(source, **kwargs):
2569 if isinstance(source, str):
2570 if op.isdir(source):
2571 return iload_dirname(source, **kwargs)
2572 elif op.isfile(source):
2573 return iload_filename(source, **kwargs)
2574 else:
2575 return iload_glob(source, **kwargs)
2577 elif hasattr(source, 'read'):
2578 return iload_fh(source, **kwargs)
2579 else:
2580 return itertools.chain.from_iterable(
2581 iload(subsource, **kwargs) for subsource in source)
2583 iload_filename.__doc__ = '''
2584 Read %s information from named file.
2585 ''' % doc_fmt
2587 iload_dirname.__doc__ = '''
2588 Read %s information from directory of %s files.
2589 ''' % (doc_fmt, doc_fmt)
2591 iload_glob.__doc__ = '''
2592 Read %s information from files matching a glob pattern.
2593 ''' % doc_fmt
2595 iload.__doc__ = '''
2596 Load %s information from given source(s)
2598 The ``source`` can be specified as the name of a %s file, the name of a
2599 directory containing %s files, a glob pattern of %s files, an open
2600 filehandle or an iterator yielding any of the forementioned sources.
2602 This function behaves as a generator yielding %s objects.
2603 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects)
2605 for f in iload_filename, iload_dirname, iload_glob, iload:
2606 f.__module__ = iload_fh.__module__
2608 return iload_filename, iload_dirname, iload_glob, iload
2611class Inconsistency(Exception):
2612 pass
2615def consistency_check(list_of_tuples, message='values differ:'):
2616 '''
2617 Check for inconsistencies.
2619 Given a list of tuples, check that all tuple elements except for first one
2620 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be
2621 valid because the coordinates at the two channels are the same.
2622 '''
2624 if len(list_of_tuples) >= 2:
2625 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]):
2626 raise Inconsistency('%s\n' % message + '\n'.join(
2627 ' %s: %s' % (t[0], ', '.join(str(x) for x in t[1:]))
2628 for t in list_of_tuples))
2631class defaultzerodict(dict):
2632 def __missing__(self, k):
2633 return 0
2636def mostfrequent(x):
2637 c = defaultzerodict()
2638 for e in x:
2639 c[e] += 1
2641 return sorted(list(c.keys()), key=lambda k: c[k])[-1]
2644def consistency_merge(list_of_tuples,
2645 message='values differ:',
2646 error='raise',
2647 merge=mostfrequent):
2649 assert error in ('raise', 'warn', 'ignore')
2651 if len(list_of_tuples) == 0:
2652 raise Exception('cannot merge empty sequence')
2654 try:
2655 consistency_check(list_of_tuples, message)
2656 return list_of_tuples[0][1:]
2657 except Inconsistency as e:
2658 if error == 'raise':
2659 raise
2661 elif error == 'warn':
2662 logger.warning(str(e))
2664 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]])
2667def short_to_list(nmax, it):
2668 import itertools
2670 if isinstance(it, list):
2671 return it
2673 li = []
2674 for i in range(nmax+1):
2675 try:
2676 li.append(next(it))
2677 except StopIteration:
2678 return li
2680 return itertools.chain(li, it)
2683def parse_md(f):
2684 try:
2685 with open(op.join(
2686 op.dirname(op.abspath(f)),
2687 'README.md'), 'r') as readme:
2688 mdstr = readme.read()
2689 except IOError as e:
2690 return 'Failed to get README.md: %s' % e
2692 # Remve the title
2693 mdstr = re.sub(r'^# .*\n?', '', mdstr)
2694 # Append sphinx reference to `pyrocko.` modules
2695 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr)
2696 # Convert Subsections to toc-less rubrics
2697 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr)
2698 return mdstr
2701def mpl_show(plt):
2702 import matplotlib
2703 if matplotlib.get_backend().lower() == 'agg':
2704 logger.warning('Cannot show() when using matplotlib "agg" backend')
2705 else:
2706 plt.show()
2709g_re_qsplit = re.compile(
2710 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)')
2711g_re_qsplit_sep = {}
2714def get_re_qsplit(sep):
2715 if sep is None:
2716 return g_re_qsplit
2717 else:
2718 if sep not in g_re_qsplit_sep:
2719 assert len(sep) == 1
2720 assert sep not in '\'"'
2721 esep = re.escape(sep)
2722 g_re_qsplit_sep[sep] = re.compile(
2723 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|'
2724 + r'([^' + esep + r']+|(?<=' + esep + r')(?=' + esep + r')|^(?=' + esep + r')|(?<=' + esep + r')$)') # noqa
2725 return g_re_qsplit_sep[sep]
2728g_re_trivial = re.compile(r'\A[^\'"\s]+\Z')
2729g_re_trivial_sep = {}
2732def get_re_trivial(sep):
2733 if sep is None:
2734 return g_re_trivial
2735 else:
2736 if sep not in g_re_qsplit_sep:
2737 assert len(sep) == 1
2738 assert sep not in '\'"'
2739 esep = re.escape(sep)
2740 g_re_trivial_sep[sep] = re.compile(r'\A[^\'"' + esep + r']+\Z')
2742 return g_re_trivial_sep[sep]
2745g_re_escape_s = re.compile(r'([\\\'])')
2746g_re_unescape_s = re.compile(r'\\([\\\'])')
2747g_re_escape_d = re.compile(r'([\\"])')
2748g_re_unescape_d = re.compile(r'\\([\\"])')
2751def escape_s(s):
2752 '''
2753 Backslash-escape single-quotes and backslashes.
2755 Example: ``Jack's`` => ``Jack\\'s``
2757 '''
2758 return g_re_escape_s.sub(r'\\\1', s)
2761def unescape_s(s):
2762 '''
2763 Unescape backslash-escaped single-quotes and backslashes.
2765 Example: ``Jack\\'s`` => ``Jack's``
2766 '''
2767 return g_re_unescape_s.sub(r'\1', s)
2770def escape_d(s):
2771 '''
2772 Backslash-escape double-quotes and backslashes.
2774 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"``
2775 '''
2776 return g_re_escape_d.sub(r'\\\1', s)
2779def unescape_d(s):
2780 '''
2781 Unescape backslash-escaped double-quotes and backslashes.
2783 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"``
2784 '''
2785 return g_re_unescape_d.sub(r'\1', s)
2788def qjoin_s(it, sep=None):
2789 '''
2790 Join sequence of strings into a line, single-quoting non-trivial strings.
2792 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"``
2793 '''
2794 re_trivial = get_re_trivial(sep)
2796 if sep is None:
2797 sep = ' '
2799 return sep.join(
2800 w if re_trivial.search(w) else "'%s'" % escape_s(w) for w in it)
2803def qjoin_d(it, sep=None):
2804 '''
2805 Join sequence of strings into a line, double-quoting non-trivial strings.
2807 Example: ``['55', 'Pete "The Robot" Smith']`` =>
2808 ``'55' "Pete \\\\"The Robot\\\\" Smith"'``
2809 '''
2810 re_trivial = get_re_trivial(sep)
2811 if sep is None:
2812 sep = ' '
2814 return sep.join(
2815 w if re_trivial.search(w) else '"%s"' % escape_d(w) for w in it)
2818def qsplit(s, sep=None):
2819 '''
2820 Split line into list of strings, allowing for quoted strings.
2822 Example: ``"55 'Sparrow\\\\'s Island'"`` =>
2823 ``["55", "Sparrow's Island"]``,
2824 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` =>
2825 ``['55', 'Pete "The Robot" Smith']``
2826 '''
2827 re_qsplit = get_re_qsplit(sep)
2828 return [
2829 (unescape_d(x[0]) or unescape_s(x[1]) or x[2])
2830 for x in re_qsplit.findall(s)]
2833g_have_warned_threadpoolctl = False
2836class threadpool_limits_dummy(object):
2838 def __init__(self, *args, **kwargs):
2839 pass
2841 def __enter__(self):
2842 global g_have_warned_threadpoolctl
2844 if not g_have_warned_threadpoolctl:
2845 logger.warning(
2846 'Cannot control number of BLAS threads because '
2847 '`threadpoolctl` module is not available. You may want to '
2848 'install `threadpoolctl`.')
2850 g_have_warned_threadpoolctl = True
2852 return self
2854 def __exit__(self, type, value, traceback):
2855 pass
2858def get_threadpool_limits():
2859 '''
2860 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail.
2861 '''
2863 try:
2864 from threadpoolctl import threadpool_limits
2865 return threadpool_limits
2867 except ImportError:
2868 return threadpool_limits_dummy
2871def fmt_summary(entries, widths):
2872 return ' | '.join(
2873 entry.ljust(width) for (entry, width) in zip(entries, widths))
2876def smart_weakref(obj, callback=None):
2877 if inspect.ismethod(obj):
2878 return weakref.WeakMethod(obj, callback)
2879 else:
2880 return weakref.ref(obj, callback)