1# http://pyrocko.org - GPLv3
3#
4# The Pyrocko Developers, 21st Century
5# ---|P------/S----------~Lg----------
6'''
7Utility functions for Pyrocko.
9.. _time-handling-mode:
11High precision time handling mode
12.................................
14Pyrocko can treat timestamps either as standard double precision (64 bit)
15floating point values, or as high precision float (:py:class:`numpy.float128`
16or :py:class:`numpy.float96`, whichever is available), aliased here as
17:py:class:`pyrocko.util.hpfloat`. High precision time stamps are required when
18handling data with sub-millisecond precision, i.e. kHz/MHz data streams and
19event catalogs derived from such data.
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 config variable
28:py:gattr:`~pyrocko.config.Config.use_high_precision_time`. Calling the
29function :py:func:`use_high_precision_time` overrides the default from the
30config file. This function may be called at startup of a program/script
31requiring a specific time handling mode.
33To create a valid time stamp for use in Pyrocko (e.g. in
34:py:class:`~pyrocko.model.Event` or :py:class:`~pyrocko.trace.Trace` objects),
35use:
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..............
62'''
64from __future__ import division, print_function
66import time
67import logging
68import os
69import sys
70import re
71import calendar
72import math
73import fnmatch
74try:
75 import fcntl
76except ImportError:
77 fcntl = None
78import optparse
79import os.path as op
80import errno
82import numpy as num
83from scipy import signal
84import pyrocko
85from pyrocko import dummy_progressbar
88try:
89 from urllib.parse import urlencode, quote, unquote # noqa
90 from urllib.request import (
91 Request, build_opener, HTTPDigestAuthHandler, urlopen as _urlopen) # noqa
92 from urllib.error import HTTPError, URLError # noqa
94except ImportError:
95 from urllib import urlencode, quote, unquote # noqa
96 from urllib2 import (Request, build_opener, HTTPDigestAuthHandler, # noqa
97 HTTPError, URLError, urlopen as _urlopen) # noqa
100class URLErrorSSL(URLError):
102 def __init__(self, urlerror):
103 self.__dict__ = urlerror.__dict__.copy()
105 def __str__(self):
106 return (
107 'Requesting web resource failed and the problem could be '
108 'related to SSL. Python standard libraries on some older '
109 'systems (like Ubuntu 14.04) are known to have trouble '
110 'with some SSL setups of today\'s servers: %s'
111 % URLError.__str__(self))
114def urlopen_ssl_check(*args, **kwargs):
115 try:
116 return urlopen(*args, **kwargs)
117 except URLError as e:
118 if str(e).find('SSL') != -1:
119 raise URLErrorSSL(e)
120 else:
121 raise
124def urlopen(*args, **kwargs):
125 if 'cafile' not in kwargs:
126 try:
127 import certifi
128 kwargs['cafile'] = certifi.where()
129 except ImportError:
130 pass
132 return _urlopen(*args, **kwargs)
135try:
136 long
137except NameError:
138 long = int
141force_dummy_progressbar = False
144try:
145 from pyrocko import util_ext
146except ImportError:
147 util_ext = None
150logger = logging.getLogger('pyrocko.util')
152try:
153 import progressbar as progressbar_mod
154except ImportError:
155 from pyrocko import dummy_progressbar as progressbar_mod
158try:
159 num_full = num.full
160except AttributeError:
161 def num_full(shape, fill_value, dtype=None, order='C'):
162 a = num.empty(shape, dtype=dtype, order=order)
163 a.fill(fill_value)
164 return a
166try:
167 num_full_like = num.full_like
168except AttributeError:
169 def num_full_like(arr, fill_value, dtype=None, order='K', subok=True):
170 a = num.empty_like(arr, dtype=dtype, order=order, subok=subok)
171 a.fill(fill_value)
172 return a
175def progressbar_module():
176 return progressbar_mod
179g_setup_logging_args = 'pyrocko', 'warning'
182def setup_logging(programname='pyrocko', levelname='warning'):
183 '''
184 Initialize logging.
186 :param programname: program name to be written in log
187 :param levelname: string indicating the logging level ('debug', 'info',
188 'warning', 'error', 'critical')
190 This function is called at startup by most pyrocko programs to set up a
191 consistent logging format. This is simply a shortcut to a call to
192 :py:func:`logging.basicConfig()`.
193 '''
195 global g_setup_logging_args
196 g_setup_logging_args = (programname, levelname)
198 levels = {'debug': logging.DEBUG,
199 'info': logging.INFO,
200 'warning': logging.WARNING,
201 'error': logging.ERROR,
202 'critical': logging.CRITICAL}
204 logging.basicConfig(
205 level=levels[levelname],
206 format=programname+':%(name)-25s - %(levelname)-8s - %(message)s')
209def subprocess_setup_logging_args():
210 '''
211 Get arguments from previous call to setup_logging.
213 These can be sent down to a worker process so it can setup its logging
214 in the same way as the main process.
215 '''
216 return g_setup_logging_args
219def data_file(fn):
220 return os.path.join(os.path.split(__file__)[0], 'data', fn)
223class DownloadError(Exception):
224 pass
227class PathExists(DownloadError):
228 pass
231def _download(url, fpath, username=None, password=None,
232 force=False, method='download', stats=None,
233 status_callback=None, entries_wanted=None,
234 recursive=False, header=None):
236 import requests
237 from requests.auth import HTTPBasicAuth
238 from requests.exceptions import HTTPError as req_HTTPError
240 requests.adapters.DEFAULT_RETRIES = 5
241 urljoin = requests.compat.urljoin
243 session = requests.Session()
244 session.header = header
245 session.auth = None if username is None\
246 else HTTPBasicAuth(username, password)
248 status = {
249 'ntotal_files': 0,
250 'nread_files': 0,
251 'ntotal_bytes_all_files': 0,
252 'nread_bytes_all_files': 0,
253 'ntotal_bytes_current_file': 0,
254 'nread_bytes_current_file': 0,
255 'finished': False
256 }
258 try:
259 url_to_size = {}
261 if callable(status_callback):
262 status_callback(status)
264 if not recursive and url.endswith('/'):
265 raise DownloadError(
266 'URL: %s appears to be a directory'
267 ' but recurvise download is False' % url)
269 if recursive and not url.endswith('/'):
270 url += '/'
272 def parse_directory_tree(url, subdir=''):
273 r = session.get(urljoin(url, subdir))
274 r.raise_for_status()
276 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text)
278 files = sorted(set(subdir + fn for fn in entries
279 if not fn.endswith('/')))
281 if entries_wanted is not None:
282 files = [fn for fn in files
283 if (fn in wanted for wanted in entries_wanted)]
285 status['ntotal_files'] += len(files)
287 dirs = sorted(set(subdir + dn for dn in entries
288 if dn.endswith('/')
289 and dn not in ('./', '../')))
291 for dn in dirs:
292 files.extend(parse_directory_tree(
293 url, subdir=dn))
295 return files
297 def get_content_length(url):
298 if url not in url_to_size:
299 r = session.head(url, headers={'Accept-Encoding': ''})
301 content_length = r.headers.get('content-length', None)
302 if content_length is None:
303 logger.debug('Could not get HTTP header '
304 'Content-Length for %s' % url)
306 content_length = None
308 else:
309 content_length = int(content_length)
310 status['ntotal_bytes_all_files'] += content_length
311 if callable(status_callback):
312 status_callback(status)
314 url_to_size[url] = content_length
316 return url_to_size[url]
318 def download_file(url, fn):
319 logger.info('starting download of %s...' % url)
320 ensuredirs(fn)
322 fsize = get_content_length(url)
323 status['ntotal_bytes_current_file'] = fsize
324 status['nread_bytes_current_file'] = 0
325 if callable(status_callback):
326 status_callback(status)
328 r = session.get(url, stream=True, timeout=5)
329 r.raise_for_status()
331 frx = 0
332 fn_tmp = fn + '.%i.temp' % os.getpid()
333 with open(fn_tmp, 'wb') as f:
334 for d in r.iter_content(chunk_size=1024):
335 f.write(d)
336 frx += len(d)
338 status['nread_bytes_all_files'] += len(d)
339 status['nread_bytes_current_file'] += len(d)
340 if callable(status_callback):
341 status_callback(status)
343 os.rename(fn_tmp, fn)
345 if fsize is not None and frx != fsize:
346 logger.warning(
347 'HTTP header Content-Length: %i bytes does not match '
348 'download size %i bytes' % (fsize, frx))
350 logger.info('finished download of %s' % url)
352 status['nread_files'] += 1
353 if callable(status_callback):
354 status_callback(status)
356 if recursive:
357 if op.exists(fpath) and not force:
358 raise PathExists('path %s already exists' % fpath)
360 files = parse_directory_tree(url)
362 dsize = 0
363 for fn in files:
364 file_url = urljoin(url, fn)
365 dsize += get_content_length(file_url) or 0
367 if method == 'calcsize':
368 return dsize
370 else:
371 for fn in files:
372 file_url = urljoin(url, fn)
373 download_file(file_url, op.join(fpath, fn))
375 else:
376 status['ntotal_files'] += 1
377 if callable(status_callback):
378 status_callback(status)
380 fsize = get_content_length(url)
381 if method == 'calcsize':
382 return fsize
383 else:
384 download_file(url, fpath)
386 except req_HTTPError as e:
387 logging.warning("http error: %s" % e)
388 raise DownloadError('could not download file(s) from: %s' % url)
390 finally:
391 status['finished'] = True
392 if callable(status_callback):
393 status_callback(status)
394 session.close()
397def download_file(
398 url, fpath, username=None, password=None, status_callback=None,
399 **kwargs):
400 return _download(
401 url, fpath, username, password,
402 recursive=False,
403 status_callback=status_callback,
404 **kwargs)
407def download_dir(
408 url, fpath, username=None, password=None, status_callback=None,
409 **kwargs):
411 return _download(
412 url, fpath, username, password,
413 recursive=True,
414 status_callback=status_callback,
415 **kwargs)
418class HPFloatUnavailable(Exception):
419 pass
422class dummy_hpfloat(object):
423 def __init__(self, *args, **kwargs):
424 raise HPFloatUnavailable(
425 'NumPy lacks support for float128 or float96 data type on this '
426 'platform.')
429if hasattr(num, 'float128'):
430 hpfloat = num.float128
432elif hasattr(num, 'float96'):
433 hpfloat = num.float96
435else:
436 hpfloat = dummy_hpfloat
439g_time_float = None
440g_time_dtype = None
443class TimeFloatSettingError(Exception):
444 pass
447def use_high_precision_time(enabled):
448 '''
449 Globally force a specific time handling mode.
451 See :ref:`High precision time handling mode <time-handling-mode>`.
453 :param enabled: enable/disable use of high precision time type
454 :type enabled: bool
456 This function should be called before handling/reading any time data.
457 It can only be called once.
459 Special attention is required when using multiprocessing on a platform
460 which does not use fork under the hood. In such cases, the desired setting
461 must be set also in the subprocess.
462 '''
463 _setup_high_precision_time_mode(enabled_app=enabled)
466def _setup_high_precision_time_mode(enabled_app=False):
467 global g_time_float
468 global g_time_dtype
470 if not (g_time_float is None and g_time_dtype is None):
471 raise TimeFloatSettingError(
472 'Cannot set time handling mode: too late, it has already been '
473 'fixed by an earlier call.')
475 from pyrocko import config
477 conf = config.config()
478 enabled_config = conf.use_high_precision_time
480 enabled_env = os.environ.get('PYROCKO_USE_HIGH_PRECISION_TIME', None)
481 if enabled_env is not None:
482 try:
483 enabled_env = int(enabled_env) == 1
484 except ValueError:
485 raise TimeFloatSettingError(
486 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME '
487 'should be set to 0 or 1.')
489 enabled = enabled_config
490 mode_from = 'config variable `use_high_precision_time`'
491 notify = enabled
493 if enabled_env is not None:
494 if enabled_env != enabled:
495 notify = True
496 enabled = enabled_env
497 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`'
499 if enabled_app is not None:
500 if enabled_app != enabled:
501 notify = True
502 enabled = enabled_app
503 mode_from = 'application override'
505 logger.debug('''
506Pyrocko high precision time mode selection (latter override earlier):
507 config: %s
508 env: %s
509 app: %s
510 -> enabled: %s'''.lstrip() % (
511 enabled_config, enabled_env, enabled_app, enabled))
513 if notify:
514 logger.info('Pyrocko high precision time mode %s by %s.' % (
515 'activated' if enabled else 'deactivated',
516 mode_from))
518 if enabled:
519 g_time_float = hpfloat
520 g_time_dtype = hpfloat
521 else:
522 g_time_float = float
523 g_time_dtype = num.float64
526def get_time_float():
527 '''
528 Get the effective float class for timestamps.
530 See :ref:`High precision time handling mode <time-handling-mode>`.
532 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the
533 current time handling mode
534 '''
535 global g_time_float
537 if g_time_float is None:
538 _setup_high_precision_time_mode()
540 return g_time_float
543def get_time_dtype():
544 '''
545 Get effective NumPy float class to handle timestamps.
547 See :ref:`High precision time handling mode <time-handling-mode>`.
548 '''
550 global g_time_dtype
552 if g_time_dtype is None:
553 _setup_high_precision_time_mode()
555 return g_time_dtype
558def to_time_float(t):
559 '''
560 Convert float to valid time stamp in the current time handling mode.
562 See :ref:`High precision time handling mode <time-handling-mode>`.
563 '''
564 return get_time_float()(t)
567class TimestampTypeError(ValueError):
568 pass
571def check_time_class(t, error='raise'):
572 '''
573 Type-check variable against current time handling mode.
575 See :ref:`High precision time handling mode <time-handling-mode>`.
576 '''
578 if t == 0.0:
579 return
581 if not isinstance(t, get_time_float()):
582 message = (
583 'Timestamp %g is of type %s but should be of type %s with '
584 'Pyrocko\'s currently selected time handling mode.\n\n'
585 'See https://pyrocko.org/docs/current/library/reference/util.html'
586 '#high-precision-time-handling-mode' % (
587 t, type(t), get_time_float()))
589 if error == 'raise':
590 raise TimestampTypeError(message)
591 elif error == 'warn':
592 logger.warning(message)
593 else:
594 assert False
597class Stopwatch(object):
598 '''
599 Simple stopwatch to measure elapsed wall clock time.
601 Usage::
603 s = Stopwatch()
604 time.sleep(1)
605 print s()
606 time.sleep(1)
607 print s()
608 '''
610 def __init__(self):
611 self.start = time.time()
613 def __call__(self):
614 return time.time() - self.start
617def wrap(text, line_length=80):
618 '''
619 Paragraph and list-aware wrapping of text.
620 '''
622 text = text.strip('\n')
623 at_lineend = re.compile(r' *\n')
624 at_para = re.compile(r'((^|(\n\s*)?\n)(\s+[*] )|\n\s*\n)')
626 paragraphs = at_para.split(text)[::5]
627 listindents = at_para.split(text)[4::5]
628 newlist = at_para.split(text)[3::5]
630 listindents[0:0] = [None]
631 listindents.append(True)
632 newlist.append(None)
634 det_indent = re.compile(r'^ *')
636 outlines = []
637 for ip, p in enumerate(paragraphs):
638 if not p:
639 continue
641 if listindents[ip] is None:
642 _indent = det_indent.findall(p)[0]
643 findent = _indent
644 else:
645 findent = listindents[ip]
646 _indent = ' ' * len(findent)
648 ll = line_length - len(_indent)
649 llf = ll
651 oldlines = [s.strip() for s in at_lineend.split(p.rstrip())]
652 p1 = ' '.join(oldlines)
653 possible = re.compile(r'(^.{1,%i}|.{1,%i})( |$)' % (llf, ll))
654 for imatch, match in enumerate(possible.finditer(p1)):
655 parout = match.group(1)
656 if imatch == 0:
657 outlines.append(findent + parout)
658 else:
659 outlines.append(_indent + parout)
661 if ip != len(paragraphs)-1 and (
662 listindents[ip] is None
663 or newlist[ip] is not None
664 or listindents[ip+1] is None):
666 outlines.append('')
668 return outlines
671def ewrap(lines, width=80, indent=''):
672 lines = list(lines)
673 if not lines:
674 return ''
675 fwidth = max(len(s) for s in lines)
676 nx = max(1, (80-len(indent)) // (fwidth+1))
677 i = 0
678 rows = []
679 while i < len(lines):
680 rows.append(indent + ' '.join(x.ljust(fwidth) for x in lines[i:i+nx]))
681 i += nx
683 return '\n'.join(rows)
686class BetterHelpFormatter(optparse.IndentedHelpFormatter):
688 def __init__(self, *args, **kwargs):
689 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs)
691 def format_option(self, option):
692 '''
693 From IndentedHelpFormatter but using a different wrap method.
694 '''
696 help_text_position = 4 + self.current_indent
697 help_text_width = self.width - help_text_position
699 opts = self.option_strings[option]
700 opts = "%*s%s" % (self.current_indent, "", opts)
701 if option.help:
702 help_text = self.expand_default(option)
704 if self.help_position + len(help_text) + 1 <= self.width:
705 lines = [
706 '',
707 '%-*s %s' % (self.help_position, opts, help_text),
708 '']
709 else:
710 lines = ['']
711 lines.append(opts)
712 lines.append('')
713 if option.help:
714 help_lines = wrap(help_text, help_text_width)
715 lines.extend(["%*s%s" % (help_text_position, "", line)
716 for line in help_lines])
717 lines.append('')
719 return "\n".join(lines)
721 def format_description(self, description):
722 if not description:
723 return ''
725 if self.current_indent == 0:
726 lines = []
727 else:
728 lines = ['']
730 lines.extend(wrap(description, self.width - self.current_indent))
731 if self.current_indent == 0:
732 lines.append('\n')
734 return '\n'.join(
735 ['%*s%s' % (self.current_indent, '', line) for line in lines])
738def progressbar(label, maxval):
739 progressbar_mod = progressbar_module()
740 if force_dummy_progressbar:
741 progressbar_mod = dummy_progressbar
743 widgets = [
744 label, ' ',
745 progressbar_mod.Bar(marker='-', left='[', right=']'), ' ',
746 progressbar_mod.Percentage(), ' ',
747 progressbar_mod.ETA()]
749 pbar = progressbar_mod.ProgressBar(widgets=widgets, maxval=maxval).start()
750 return pbar
753def progress_beg(label):
754 '''
755 Notify user that an operation has started.
757 :param label: name of the operation
759 To be used in conjuction with :py:func:`progress_end`.
760 '''
762 sys.stderr.write(label)
763 sys.stderr.flush()
766def progress_end(label=''):
767 '''
768 Notify user that an operation has ended.
770 :param label: name of the operation
772 To be used in conjuction with :py:func:`progress_beg`.
773 '''
775 sys.stderr.write(' done. %s\n' % label)
776 sys.stderr.flush()
779class ArangeError(Exception):
780 pass
783def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'):
784 '''
785 Return evenly spaced numbers over a specified interval.
787 Like :py:func:`numpy.arange` but returning floating point numbers by
788 default and with defined behaviour when stepsize is inconsistent with
789 interval bounds. It is considered inconsistent if the difference between
790 the closest multiple of ``step`` and ``stop`` is larger than ``epsilon *
791 step``. Inconsistencies are handled according to the ``error`` parameter.
792 If it is set to ``'raise'`` an exception of type :py:exc:`ArangeError` is
793 raised. If it is set to ``'round'``, ``'floor'``, or ``'ceil'``, ``stop``
794 is silently changed to the closest, the next smaller, or next larger
795 multiple of ``step``, respectively.
796 '''
798 assert error in ('raise', 'round', 'floor', 'ceil')
800 start = dtype(start)
801 stop = dtype(stop)
802 step = dtype(step)
804 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round)
806 n = int(rnd((stop - start) / step)) + 1
807 stop_check = start + (n-1) * step
809 if error == 'raise' and abs(stop_check - stop) > step * epsilon:
810 raise ArangeError(
811 'inconsistent range specification: start=%g, stop=%g, step=%g'
812 % (start, stop, step))
814 x = num.arange(n, dtype=dtype)
815 x *= step
816 x += start
817 return x
820def polylinefit(x, y, n_or_xnodes):
821 '''
822 Fit piece-wise linear function to data.
824 :param x,y: arrays with coordinates of data
825 :param n_or_xnodes: int, number of segments or x coordinates of polyline
827 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of
828 polyline, root-mean-square error
829 '''
831 x = num.asarray(x)
832 y = num.asarray(y)
834 if isinstance(n_or_xnodes, int):
835 n = n_or_xnodes
836 xmin = x.min()
837 xmax = x.max()
838 xnodes = num.linspace(xmin, xmax, n+1)
839 else:
840 xnodes = num.asarray(n_or_xnodes)
841 n = xnodes.size - 1
843 assert len(x) == len(y)
844 assert n > 0
846 ndata = len(x)
847 a = num.zeros((ndata+(n-1), n*2))
848 for i in range(n):
849 xmin_block = xnodes[i]
850 xmax_block = xnodes[i+1]
851 if i == n-1: # don't loose last point
852 indices = num.where(
853 num.logical_and(xmin_block <= x, x <= xmax_block))[0]
854 else:
855 indices = num.where(
856 num.logical_and(xmin_block <= x, x < xmax_block))[0]
858 a[indices, i*2] = x[indices]
859 a[indices, i*2+1] = 1.0
861 w = float(ndata)*100.
862 if i < n-1:
863 a[ndata+i, i*2] = xmax_block*w
864 a[ndata+i, i*2+1] = 1.0*w
865 a[ndata+i, i*2+2] = -xmax_block*w
866 a[ndata+i, i*2+3] = -1.0*w
868 d = num.concatenate((y, num.zeros(n-1)))
869 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2))
871 ynodes = num.zeros(n+1)
872 ynodes[:n] = model[:, 0]*xnodes[:n] + model[:, 1]
873 ynodes[1:] += model[:, 0]*xnodes[1:] + model[:, 1]
874 ynodes[1:n] *= 0.5
876 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2))
878 return xnodes, ynodes, rms_error
881def plf_integrate_piecewise(x_edges, x, y):
882 '''
883 Calculate definite integral of piece-wise linear function on intervals.
885 Use trapezoidal rule to calculate definite integral of a piece-wise linear
886 function for a series of consecutive intervals. ``x_edges`` and ``x`` must
887 be sorted.
889 :param x_edges: array with edges of the intervals
890 :param x,y: arrays with coordinates of piece-wise linear function's
891 control points
892 '''
894 x_all = num.concatenate((x, x_edges))
895 ii = num.argsort(x_all)
896 y_edges = num.interp(x_edges, x, y)
897 y_all = num.concatenate((y, y_edges))
898 xs = x_all[ii]
899 ys = y_all[ii]
900 y_all[ii[1:]] = num.cumsum((xs[1:] - xs[:-1]) * 0.5 * (ys[1:] + ys[:-1]))
901 return num.diff(y_all[-len(y_edges):])
904def diff_fd_1d_4o(dt, data):
905 '''
906 Approximate first derivative of an array (forth order, central FD).
908 :param dt: sampling interval
909 :param data: NumPy array with data samples
911 :returns: NumPy array with same shape as input
913 Interior points are approximated to fourth order, edge points to first
914 order right- or left-sided respectively, points next to edge to second
915 order central.
916 '''
917 import scipy.signal
919 ddata = num.empty_like(data, dtype=float)
921 if data.size >= 5:
922 ddata[2:-2] = scipy.signal.lfilter(
923 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt)
925 if data.size >= 3:
926 ddata[1] = (data[2] - data[0]) / (2. * dt)
927 ddata[-2] = (data[-1] - data[-3]) / (2. * dt)
929 if data.size >= 2:
930 ddata[0] = (data[1] - data[0]) / dt
931 ddata[-1] = (data[-1] - data[-2]) / dt
933 if data.size == 1:
934 ddata[0] = 0.0
936 return ddata
939def diff_fd_1d_2o(dt, data):
940 '''
941 Approximate first derivative of an array (second order, central FD).
943 :param dt: sampling interval
944 :param data: NumPy array with data samples
946 :returns: NumPy array with same shape as input
948 Interior points are approximated to second order, edge points to first
949 order right- or left-sided respectively.
951 Uses :py:func:`numpy.gradient`.
952 '''
954 return num.gradient(data, dt)
957def diff_fd_2d_4o(dt, data):
958 '''
959 Approximate second derivative of an array (forth order, central FD).
961 :param dt: sampling interval
962 :param data: NumPy array with data samples
964 :returns: NumPy array with same shape as input
966 Interior points are approximated to fourth order, next-to-edge points to
967 second order, edge points repeated.
968 '''
969 import scipy.signal
971 ddata = num.empty_like(data, dtype=float)
973 if data.size >= 5:
974 ddata[2:-2] = scipy.signal.lfilter(
975 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2)
977 if data.size >= 3:
978 ddata[:2] = (data[2] - 2.0 * data[1] + data[0]) / dt**2
979 ddata[-2:] = (data[-1] - 2.0 * data[-2] + data[-3]) / dt**2
981 if data.size < 3:
982 ddata[:] = 0.0
984 return ddata
987def diff_fd_2d_2o(dt, data):
988 '''
989 Approximate second derivative of an array (second order, central FD).
991 :param dt: sampling interval
992 :param data: NumPy array with data samples
994 :returns: NumPy array with same shape as input
996 Interior points are approximated to second order, edge points repeated.
997 '''
998 import scipy.signal
1000 ddata = num.empty_like(data, dtype=float)
1002 if data.size >= 3:
1003 ddata[1:-1] = scipy.signal.lfilter(
1004 [1., -2., 1.], [1.], data)[2:] / (dt**2)
1006 ddata[0] = ddata[1]
1007 ddata[-1] = ddata[-2]
1009 if data.size < 3:
1010 ddata[:] = 0.0
1012 return ddata
1015def diff_fd(n, order, dt, data):
1016 '''
1017 Approximate 1st or 2nd derivative of an array.
1019 :param n: 1 for first derivative, 2 for second
1020 :param order: order of the approximation 2 and 4 are supported
1021 :param dt: sampling interval
1022 :param data: NumPy array with data samples
1024 :returns: NumPy array with same shape as input
1026 This is a frontend to the functions :py:func:`diff_fd_1d_2o`,
1027 :py:func:`diff_fd_1d_4o`, :py:func:`diff_fd_2d_2o`, and
1028 :py:func:`diff_fd_2d_4o`.
1030 Raises :py:exc:`ValueError` for unsupported `n` or `order`.
1031 '''
1033 funcs = {
1034 1: {2: diff_fd_1d_2o, 4: diff_fd_1d_4o},
1035 2: {2: diff_fd_2d_2o, 4: diff_fd_2d_4o}}
1037 try:
1038 funcs_n = funcs[n]
1039 except KeyError:
1040 raise ValueError(
1041 'pyrocko.util.diff_fd: '
1042 'Only 1st and 2sd derivatives are supported.')
1044 try:
1045 func = funcs_n[order]
1046 except KeyError:
1047 raise ValueError(
1048 'pyrocko.util.diff_fd: '
1049 'Order %i is not supported for %s derivative. Supported: %s' % (
1050 order, ['', '1st', '2nd'][n],
1051 ', '.join('%i' % order for order in sorted(funcs_n.keys()))))
1053 return func(dt, data)
1056class GlobalVars(object):
1057 reuse_store = dict()
1058 decitab_nmax = 0
1059 decitab = {}
1060 decimate_fir_coeffs = {}
1061 decimate_fir_remez_coeffs = {}
1062 decimate_iir_coeffs = {}
1063 re_frac = None
1066def decimate_coeffs(q, n=None, ftype='iir'):
1068 q = int(q)
1070 if n is None:
1071 if ftype == 'fir':
1072 n = 30
1073 elif ftype == 'fir-remez':
1074 n = 40*q
1075 else:
1076 n = 8
1078 if ftype == 'fir':
1079 coeffs = GlobalVars.decimate_fir_coeffs
1080 if (n, 1./q) not in coeffs:
1081 coeffs[n, 1./q] = signal.firwin(n+1, .75/q, window='hamming', fs=2)
1083 b = coeffs[n, 1./q]
1084 return b, [1.], n
1086 elif ftype == 'fir-remez':
1087 coeffs = GlobalVars.decimate_fir_remez_coeffs
1088 if (n, 1./q) not in coeffs:
1089 coeffs[n, 1./q] = signal.remez(
1090 n+1, (0., .75/q, 1./q, 1.),
1091 (1., 0.), fs=2, weight=(1, 50))
1092 b = coeffs[n, 1./q]
1093 return b, [1.], n
1095 else:
1096 coeffs = GlobalVars.decimate_iir_coeffs
1097 if (n, 0.05, 0.8/q) not in coeffs:
1098 coeffs[n, 0.05, 0.8/q] = signal.cheby1(n, 0.05, 0.8/q)
1100 b, a = coeffs[n, 0.05, 0.8/q]
1101 return b, a, n
1104def decimate(x, q, n=None, ftype='iir', zi=None, ioff=0):
1105 '''
1106 Downsample the signal x by an integer factor q, using an order n filter
1108 By default, an order 8 Chebyshev type I filter is used or a 30 point FIR
1109 filter with hamming window if ftype is 'fir'.
1111 :param x: the signal to be downsampled (1D :class:`numpy.ndarray`)
1112 :param q: the downsampling factor
1113 :param n: order of the filter (1 less than the length of the filter for a
1114 `fir` filter)
1115 :param ftype: type of the filter; can be `iir`, `fir` or `fir-remez`
1117 :returns: the downsampled signal (1D :class:`numpy.ndarray`)
1119 '''
1121 b, a, n = decimate_coeffs(q, n, ftype)
1123 if zi is None or zi is True:
1124 zi_ = num.zeros(max(len(a), len(b))-1, dtype=float)
1125 else:
1126 zi_ = zi
1128 y, zf = signal.lfilter(b, a, x, zi=zi_)
1130 if zi is not None:
1131 return y[n//2+ioff::q].copy(), zf
1132 else:
1133 return y[n//2+ioff::q].copy()
1136class UnavailableDecimation(Exception):
1137 '''
1138 Exception raised by :py:func:`decitab` for unavailable decimation factors.
1139 '''
1141 pass
1144def gcd(a, b, epsilon=1e-7):
1145 '''
1146 Greatest common divisor.
1147 '''
1149 while b > epsilon*a:
1150 a, b = b, a % b
1152 return a
1155def lcm(a, b):
1156 '''
1157 Least common multiple.
1158 '''
1160 return a*b // gcd(a, b)
1163def mk_decitab(nmax=100):
1164 '''
1165 Make table with decimation sequences.
1167 Decimation from one sampling rate to a lower one is achieved by a
1168 successive application of :py:func:`decimation` with small integer
1169 downsampling factors (because using large downampling factors can make the
1170 decimation unstable or slow). This function sets up a table with downsample
1171 sequences for factors up to ``nmax``.
1172 '''
1174 tab = GlobalVars.decitab
1175 for i in range(1, 10):
1176 for j in range(1, i+1):
1177 for k in range(1, j+1):
1178 for l_ in range(1, k+1):
1179 for m in range(1, l_+1):
1180 p = i*j*k*l_*m
1181 if p > nmax:
1182 break
1183 if p not in tab:
1184 tab[p] = (i, j, k, l_, m)
1185 if i*j*k*l_ > nmax:
1186 break
1187 if i*j*k > nmax:
1188 break
1189 if i*j > nmax:
1190 break
1191 if i > nmax:
1192 break
1194 GlobalVars.decitab_nmax = nmax
1197def zfmt(n):
1198 return '%%0%ii' % (int(math.log10(n - 1)) + 1)
1201def _year_to_time(year):
1202 tt = (year, 1, 1, 0, 0, 0)
1203 return to_time_float(calendar.timegm(tt))
1206def _working_year(year):
1207 try:
1208 tt = (year, 1, 1, 0, 0, 0)
1209 t = calendar.timegm(tt)
1210 tt2_ = time.gmtime(t)
1211 tt2 = tuple(tt2_)[:6]
1212 if tt != tt2:
1213 return False
1215 s = '%i-01-01 00:00:00' % year
1216 s2 = time.strftime('%Y-%m-%d %H:%M:%S', tt2_)
1217 if s != s2:
1218 return False
1220 t3 = str_to_time(s2, format='%Y-%m-%d %H:%M:%S')
1221 s3 = time_to_str(t3, format='%Y-%m-%d %H:%M:%S')
1222 if s3 != s2:
1223 return False
1225 except Exception:
1226 return False
1228 return True
1231def working_system_time_range(year_min_lim=None, year_max_lim=None):
1232 '''
1233 Check time range supported by the systems's time conversion functions.
1235 Returns system time stamps of start of year of first/last fully supported
1236 year span. If this is before 1900 or after 2100, return first/last century
1237 which is fully supported.
1239 :returns: ``(tmin, tmax, year_min, year_max)``
1240 '''
1242 year0 = 2000
1243 year_min = year0
1244 year_max = year0
1246 itests = list(range(101))
1247 for i in range(19):
1248 itests.append(200 + i*100)
1250 for i in itests:
1251 year = year0 - i
1252 if year_min_lim is not None and year < year_min_lim:
1253 year_min = year_min_lim
1254 break
1255 elif not _working_year(year):
1256 break
1257 else:
1258 year_min = year
1260 for i in itests:
1261 year = year0 + i
1262 if year_max_lim is not None and year > year_max_lim:
1263 year_max = year_max_lim
1264 break
1265 elif not _working_year(year + 1):
1266 break
1267 else:
1268 year_max = year
1270 return (
1271 _year_to_time(year_min),
1272 _year_to_time(year_max),
1273 year_min, year_max)
1276g_working_system_time_range = None
1279def get_working_system_time_range():
1280 '''
1281 Caching variant of :py:func:`working_system_time_range`.
1282 '''
1284 global g_working_system_time_range
1285 if g_working_system_time_range is None:
1286 g_working_system_time_range = working_system_time_range()
1288 return g_working_system_time_range
1291def is_working_time(t):
1292 tmin, tmax, _, _ = get_working_system_time_range()
1293 return tmin <= t <= tmax
1296def julian_day_of_year(timestamp):
1297 '''
1298 Get the day number after the 1st of January of year in ``timestamp``.
1300 :returns: day number as int
1301 '''
1303 return time.gmtime(int(timestamp)).tm_yday
1306def hour_start(timestamp):
1307 '''
1308 Get beginning of hour for any point in time.
1310 :param timestamp: time instant as system timestamp (in seconds)
1312 :returns: instant of hour start as system timestamp
1313 '''
1315 tt = time.gmtime(int(timestamp))
1316 tts = tt[0:4] + (0, 0)
1317 return to_time_float(calendar.timegm(tts))
1320def day_start(timestamp):
1321 '''
1322 Get beginning of day for any point in time.
1324 :param timestamp: time instant as system timestamp (in seconds)
1326 :returns: instant of day start as system timestamp
1327 '''
1329 tt = time.gmtime(int(timestamp))
1330 tts = tt[0:3] + (0, 0, 0)
1331 return to_time_float(calendar.timegm(tts))
1334def month_start(timestamp):
1335 '''
1336 Get beginning of month for any point in time.
1338 :param timestamp: time instant as system timestamp (in seconds)
1340 :returns: instant of month start as system timestamp
1341 '''
1343 tt = time.gmtime(int(timestamp))
1344 tts = tt[0:2] + (1, 0, 0, 0)
1345 return to_time_float(calendar.timegm(tts))
1348def year_start(timestamp):
1349 '''
1350 Get beginning of year for any point in time.
1352 :param timestamp: time instant as system timestamp (in seconds)
1354 :returns: instant of year start as system timestamp
1355 '''
1357 tt = time.gmtime(int(timestamp))
1358 tts = tt[0:1] + (1, 1, 0, 0, 0)
1359 return to_time_float(calendar.timegm(tts))
1362def iter_days(tmin, tmax):
1363 '''
1364 Yields begin and end of days until given time span is covered.
1366 :param tmin,tmax: input time span
1368 :yields: tuples with (begin, end) of days as system timestamps
1369 '''
1371 t = day_start(tmin)
1372 while t < tmax:
1373 tend = day_start(t + 26*60*60)
1374 yield t, tend
1375 t = tend
1378def iter_months(tmin, tmax):
1379 '''
1380 Yields begin and end of months until given time span is covered.
1382 :param tmin,tmax: input time span
1384 :yields: tuples with (begin, end) of months as system timestamps
1385 '''
1387 t = month_start(tmin)
1388 while t < tmax:
1389 tend = month_start(t + 24*60*60*33)
1390 yield t, tend
1391 t = tend
1394def iter_years(tmin, tmax):
1395 '''
1396 Yields begin and end of years until given time span is covered.
1398 :param tmin,tmax: input time span
1400 :yields: tuples with (begin, end) of years as system timestamps
1401 '''
1403 t = year_start(tmin)
1404 while t < tmax:
1405 tend = year_start(t + 24*60*60*369)
1406 yield t, tend
1407 t = tend
1410def today():
1411 return day_start(time.time())
1414def tomorrow():
1415 return day_start(time.time() + 24*60*60)
1418def decitab(n):
1419 '''
1420 Get integer decimation sequence for given downampling factor.
1422 :param n: target decimation factor
1424 :returns: tuple with downsampling sequence
1425 '''
1427 if n > GlobalVars.decitab_nmax:
1428 mk_decitab(n*2)
1429 if n not in GlobalVars.decitab:
1430 raise UnavailableDecimation('ratio = %g' % n)
1431 return GlobalVars.decitab[n]
1434def ctimegm(s, format="%Y-%m-%d %H:%M:%S"):
1435 '''
1436 Convert string representing UTC time to system time.
1438 :param s: string to be interpreted
1439 :param format: format string passed to :py:func:`strptime`
1441 :returns: system time stamp
1443 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime.
1445 .. note::
1446 This function is to be replaced by :py:func:`str_to_time`.
1447 '''
1449 return calendar.timegm(time.strptime(s, format))
1452def gmctime(t, format="%Y-%m-%d %H:%M:%S"):
1453 '''
1454 Get string representation from system time, UTC.
1456 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime.
1458 .. note::
1459 This function is to be repaced by :py:func:`time_to_str`.
1460 '''
1462 return time.strftime(format, time.gmtime(t))
1465def gmctime_v(t, format="%a, %d %b %Y %H:%M:%S"):
1466 '''
1467 Get string representation from system time, UTC. Same as
1468 :py:func:`gmctime` but with a more verbose default format.
1470 .. note::
1471 This function is to be replaced by :py:func:`time_to_str`.
1472 '''
1474 return time.strftime(format, time.gmtime(t))
1477def gmctime_fn(t, format="%Y-%m-%d_%H-%M-%S"):
1478 '''
1479 Get string representation from system time, UTC. Same as
1480 :py:func:`gmctime` but with a default usable in filenames.
1482 .. note::
1483 This function is to be replaced by :py:func:`time_to_str`.
1484 '''
1486 return time.strftime(format, time.gmtime(t))
1489class TimeStrError(Exception):
1490 pass
1493class FractionalSecondsMissing(TimeStrError):
1494 '''
1495 Exception raised by :py:func:`str_to_time` when the given string lacks
1496 fractional seconds.
1497 '''
1499 pass
1502class FractionalSecondsWrongNumberOfDigits(TimeStrError):
1503 '''
1504 Exception raised by :py:func:`str_to_time` when the given string has an
1505 incorrect number of digits in the fractional seconds part.
1506 '''
1508 pass
1511def _endswith_n(s, endings):
1512 for ix, x in enumerate(endings):
1513 if s.endswith(x):
1514 return ix
1515 return -1
1518def str_to_time(s, format='%Y-%m-%d %H:%M:%S.OPTFRAC'):
1519 '''
1520 Convert string representing UTC time to floating point system time.
1522 :param s: string representing UTC time
1523 :param format: time string format
1524 :returns: system time stamp as floating point value
1526 Uses the semantics of :py:func:`time.strptime` but allows for fractional
1527 seconds. If the format ends with ``'.FRAC'``, anything after a dot is
1528 interpreted as fractional seconds. If the format ends with ``'.OPTFRAC'``,
1529 the fractional part, including the dot is made optional. The latter has the
1530 consequence, that the time strings and the format may not contain any other
1531 dots. If the format ends with ``'.xFRAC'`` where x is 1, 2, or 3, it is
1532 ensured, that exactly that number of digits are present in the fractional
1533 seconds.
1534 '''
1536 time_float = get_time_float()
1538 if util_ext is not None:
1539 try:
1540 t, tfrac = util_ext.stt(s, format)
1541 except util_ext.UtilExtError as e:
1542 raise TimeStrError(
1543 '%s, string=%s, format=%s' % (str(e), s, format))
1545 return time_float(t)+tfrac
1547 fracsec = 0.
1548 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC'
1550 iend = _endswith_n(format, fixed_endings)
1551 if iend != -1:
1552 dotpos = s.rfind('.')
1553 if dotpos == -1:
1554 raise FractionalSecondsMissing(
1555 'string=%s, format=%s' % (s, format))
1557 if iend > 0 and iend != (len(s)-dotpos-1):
1558 raise FractionalSecondsWrongNumberOfDigits(
1559 'string=%s, format=%s' % (s, format))
1561 format = format[:-len(fixed_endings[iend])]
1562 fracsec = float(s[dotpos:])
1563 s = s[:dotpos]
1565 elif format.endswith('.OPTFRAC'):
1566 dotpos = s.rfind('.')
1567 format = format[:-8]
1568 if dotpos != -1 and len(s[dotpos:]) > 1:
1569 fracsec = float(s[dotpos:])
1571 if dotpos != -1:
1572 s = s[:dotpos]
1574 try:
1575 return time_float(calendar.timegm(time.strptime(s, format))) \
1576 + fracsec
1577 except ValueError as e:
1578 raise TimeStrError('%s, string=%s, format=%s' % (str(e), s, format))
1581stt = str_to_time
1584def str_to_time_fillup(s):
1585 '''
1586 Default :py:func:`str_to_time` with filling in of missing values.
1588 Allows e.g. `'2010-01-01 00:00:00'` as `'2010-01-01 00:00'`,
1589 `'2010-01-01 00'`, ..., or `'2010'`.
1590 '''
1592 if len(s) in (4, 7, 10, 13, 16):
1593 s += '0000-01-01 00:00:00'[len(s):]
1595 return str_to_time(s)
1598def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'):
1599 '''
1600 Get string representation for floating point system time.
1602 :param t: floating point system time
1603 :param format: time string format
1604 :returns: string representing UTC time
1606 Uses the semantics of :py:func:`time.strftime` but additionally allows for
1607 fractional seconds. If ``format`` contains ``'.xFRAC'``, where ``x`` is a
1608 digit between 1 and 9, this is replaced with the fractional part of ``t``
1609 with ``x`` digits precision.
1610 '''
1612 if pyrocko.grumpy > 0:
1613 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise')
1615 if isinstance(format, int):
1616 if format > 0:
1617 format = '%Y-%m-%d %H:%M:%S.' + '%iFRAC' % format
1618 else:
1619 format = '%Y-%m-%d %H:%M:%S'
1621 if util_ext is not None:
1622 t0 = num.floor(t)
1623 try:
1624 return util_ext.tts(int(t0), float(t - t0), format)
1625 except util_ext.UtilExtError as e:
1626 raise TimeStrError(
1627 '%s, timestamp=%f, format=%s' % (str(e), t, format))
1629 if not GlobalVars.re_frac:
1630 GlobalVars.re_frac = re.compile(r'\.[1-9]FRAC')
1631 GlobalVars.frac_formats = dict(
1632 [('.%sFRAC' % x, '%.'+x+'f') for x in '123456789'])
1634 ts = float(num.floor(t))
1635 tfrac = t-ts
1637 m = GlobalVars.re_frac.search(format)
1638 if m:
1639 sfrac = (GlobalVars.frac_formats[m.group(0)] % tfrac)
1640 if sfrac[0] == '1':
1641 ts += 1.
1643 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1)
1645 return time.strftime(format, time.gmtime(ts))
1648tts = time_to_str
1651def mystrftime(fmt=None, tt=None, milliseconds=0):
1653 if fmt is None:
1654 fmt = '%Y-%m-%d %H:%M:%S .%r'
1656 if tt is None:
1657 tt = time.time()
1659 fmt = fmt.replace('%r', '%03i' % int(round(milliseconds)))
1660 fmt = fmt.replace('%u', '%06i' % int(round(milliseconds*1000)))
1661 fmt = fmt.replace('%n', '%09i' % int(round(milliseconds*1000000)))
1662 return time.strftime(fmt, tt)
1665def gmtime_x(timestamp):
1666 etimestamp = float(num.floor(timestamp))
1667 tt = time.gmtime(etimestamp)
1668 ms = (timestamp-etimestamp)*1000
1669 return tt, ms
1672def plural_s(n):
1673 if n == 1:
1674 return ''
1675 else:
1676 return 's'
1679def ensuredirs(dst):
1680 '''
1681 Create all intermediate path components for a target path.
1683 :param dst: target path
1685 The leaf part of the target path is not created (use :py:func:`ensuredir`
1686 if a the target path is a directory to be created).
1687 '''
1689 d, x = os.path.split(dst.rstrip(os.sep))
1690 dirs = []
1691 while d and not os.path.exists(d):
1692 dirs.append(d)
1693 d, x = os.path.split(d)
1695 dirs.reverse()
1697 for d in dirs:
1698 try:
1699 os.mkdir(d)
1700 except OSError as e:
1701 if not e.errno == errno.EEXIST:
1702 raise
1705def ensuredir(dst):
1706 '''
1707 Create directory and all intermediate path components to it as needed.
1709 :param dst: directory name
1711 Nothing is done if the given target already exists.
1712 '''
1714 if os.path.exists(dst):
1715 return
1717 dst.rstrip(os.sep)
1719 ensuredirs(dst)
1720 try:
1721 os.mkdir(dst)
1722 except OSError as e:
1723 if not e.errno == errno.EEXIST:
1724 raise
1727def reuse(x):
1728 '''
1729 Get unique instance of an object.
1731 :param x: hashable object
1732 :returns: reference to x or an equivalent object
1734 Cache object ``x`` in a global dict for reuse, or if x already
1735 is in that dict, return a reference to it.
1736 '''
1738 grs = GlobalVars.reuse_store
1739 if x not in grs:
1740 grs[x] = x
1741 return grs[x]
1744def deuse(x):
1745 grs = GlobalVars.reuse_store
1746 if x in grs:
1747 del grs[x]
1750class Anon(object):
1751 '''
1752 Dict-to-object utility.
1754 Any given arguments are stored as attributes.
1756 Example::
1758 a = Anon(x=1, y=2)
1759 print a.x, a.y
1760 '''
1762 def __init__(self, **dict):
1763 for k in dict:
1764 self.__dict__[k] = dict[k]
1767def iter_select_files(
1768 paths,
1769 selector=None,
1770 regex=None,
1771 show_progress=True,
1772 pass_through=None):
1774 '''
1775 Recursively select files (generator variant).
1777 See :py:func:`select_files`.
1778 '''
1780 if show_progress:
1781 progress_beg('selecting files...')
1782 if logger.isEnabledFor(logging.DEBUG):
1783 sys.stderr.write('\n')
1785 ngood = 0
1786 if regex:
1787 rselector = re.compile(regex)
1789 if regex:
1790 def check(path):
1791 logger.debug("looking at filename: '%s'" % path)
1792 m = rselector.search(path)
1793 if not m:
1794 logger.debug(" regex '%s' does not match." % regex)
1795 return False
1797 infos = Anon(**m.groupdict())
1798 logger.debug(" regex '%s' matches." % regex)
1799 for k, v in m.groupdict().items():
1800 logger.debug(
1801 " attribute '%s' has value '%s'" % (k, v))
1802 if selector is None or selector(infos):
1803 return True
1804 else:
1805 logger.debug(' not selected.')
1806 return False
1808 else:
1809 def check(path):
1810 return True
1812 if isinstance(paths, str):
1813 paths = [paths]
1815 for path in paths:
1816 if pass_through and pass_through(path):
1817 if check(path):
1818 yield path
1820 elif os.path.isdir(path):
1821 for (dirpath, dirnames, filenames) in os.walk(path):
1822 dirnames.sort()
1823 filenames.sort()
1824 for filename in filenames:
1825 path = op.join(dirpath, filename)
1826 if check(path):
1827 yield os.path.abspath(path)
1828 ngood += 1
1829 else:
1830 if check(path):
1831 yield os.path.abspath(path)
1832 ngood += 1
1834 if show_progress:
1835 progress_end('%i file%s selected.' % (ngood, plural_s(ngood)))
1838def select_files(paths, selector=None, regex=None, show_progress=True):
1839 '''
1840 Recursively select files.
1842 :param paths: entry path names
1843 :param selector: callback for conditional inclusion
1844 :param regex: pattern for conditional inclusion
1845 :param show_progress: if True, indicate start and stop of processing
1846 :returns: list of path names
1848 Recursively finds all files under given entry points ``paths``. If
1849 parameter ``regex`` is a regular expression, only files with matching path
1850 names are included. If additionally parameter ``selector`` is given a
1851 callback function, only files for which the callback returns ``True`` are
1852 included. The callback should take a single argument. The callback is
1853 called with a single argument, an object, having as attributes, any named
1854 groups given in ``regex``.
1856 Examples
1858 To find all files ending in ``'.mseed'`` or ``'.msd'``::
1860 select_files(paths,
1861 regex=r'\\.(mseed|msd)$')
1863 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for
1864 the year::
1866 select_files(paths,
1867 regex=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$',
1868 selector=(lambda x: int(x.year) == 2009))
1869 '''
1870 return list(iter_select_files(paths, selector, regex, show_progress))
1873def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
1874 '''
1875 Convert positive integer to a base36 string.
1876 '''
1878 if not isinstance(number, (int, long)):
1879 raise TypeError('number must be an integer')
1880 if number < 0:
1881 raise ValueError('number must be positive')
1883 # Special case for small numbers
1884 if number < 36:
1885 return alphabet[number]
1887 base36 = ''
1888 while number != 0:
1889 number, i = divmod(number, 36)
1890 base36 = alphabet[i] + base36
1892 return base36
1895def base36decode(number):
1896 '''
1897 Decode base36 endcoded positive integer.
1898 '''
1900 return int(number, 36)
1903class UnpackError(Exception):
1904 '''
1905 Exception raised when :py:func:`unpack_fixed` encounters an error.
1906 '''
1908 pass
1911ruler = ''.join(['%-10i' % i for i in range(8)]) \
1912 + '\n' + '0123456789' * 8 + '\n'
1915def unpack_fixed(format, line, *callargs):
1916 '''
1917 Unpack fixed format string, as produced by many fortran codes.
1919 :param format: format specification
1920 :param line: string to be processed
1921 :param callargs: callbacks for callback fields in the format
1923 The format is described by a string of comma-separated fields. Each field
1924 is defined by a character for the field type followed by the field width. A
1925 questionmark may be appended to the field description to allow the argument
1926 to be optional (The data string is then allowed to be filled with blanks
1927 and ``None`` is returned in this case).
1929 The following field types are available:
1931 ==== ================================================================
1932 Type Description
1933 ==== ================================================================
1934 A string (full field width is extracted)
1935 a string (whitespace at the beginning and the end is removed)
1936 i integer value
1937 f floating point value
1938 @ special type, a callback must be given for the conversion
1939 x special field type to skip parts of the string
1940 ==== ================================================================
1941 '''
1943 ipos = 0
1944 values = []
1945 icall = 0
1946 for form in format.split(','):
1947 form = form.strip()
1948 optional = form[-1] == '?'
1949 form = form.rstrip('?')
1950 typ = form[0]
1951 ln = int(form[1:])
1952 s = line[ipos:ipos+ln]
1953 cast = {
1954 'x': None,
1955 'A': str,
1956 'a': lambda x: x.strip(),
1957 'i': int,
1958 'f': float,
1959 '@': 'extra'}[typ]
1961 if cast == 'extra':
1962 cast = callargs[icall]
1963 icall += 1
1965 if cast is not None:
1966 if optional and s.strip() == '':
1967 values.append(None)
1968 else:
1969 try:
1970 values.append(cast(s))
1971 except Exception:
1972 mark = [' '] * 80
1973 mark[ipos:ipos+ln] = ['^'] * ln
1974 mark = ''.join(mark)
1975 raise UnpackError(
1976 'Invalid cast to type "%s" at position [%i:%i] of '
1977 'line: \n%s%s\n%s' % (
1978 typ, ipos, ipos+ln, ruler, line.rstrip(), mark))
1980 ipos += ln
1982 return values
1985_pattern_cache = {}
1988def _nslc_pattern(pattern):
1989 if pattern not in _pattern_cache:
1990 rpattern = re.compile(fnmatch.translate(pattern), re.I)
1991 _pattern_cache[pattern] = rpattern
1992 else:
1993 rpattern = _pattern_cache[pattern]
1995 return rpattern
1998def match_nslc(patterns, nslc):
1999 '''
2000 Match network-station-location-channel code against pattern or list of
2001 patterns.
2003 :param patterns: pattern or list of patterns
2004 :param nslc: tuple with (network, station, location, channel) as strings
2006 :returns: ``True`` if the pattern matches or if any of the given patterns
2007 match; or ``False``.
2009 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq].
2011 Example::
2013 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True
2014 '''
2016 if isinstance(patterns, str):
2017 patterns = [patterns]
2019 if not isinstance(nslc, str):
2020 s = '.'.join(nslc)
2021 else:
2022 s = nslc
2024 for pattern in patterns:
2025 if _nslc_pattern(pattern).match(s):
2026 return True
2028 return False
2031def match_nslcs(patterns, nslcs):
2032 '''
2033 Get network-station-location-channel codes that match given pattern or any
2034 of several given patterns.
2036 :param patterns: pattern or list of patterns
2037 :param nslcs: list of (network, station, location, channel) tuples
2039 See also :py:func:`match_nslc`
2040 '''
2042 matching = []
2043 for nslc in nslcs:
2044 if match_nslc(patterns, nslc):
2045 matching.append(nslc)
2047 return matching
2050class Timeout(Exception):
2051 pass
2054def create_lockfile(fn, timeout=None, timewarn=10.):
2055 t0 = time.time()
2057 while True:
2058 try:
2059 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL)
2060 os.close(f)
2061 return
2063 except OSError as e:
2064 if e.errno == errno.EEXIST:
2065 tnow = time.time()
2067 if timeout is not None and tnow - t0 > timeout:
2068 raise Timeout(
2069 'Timeout (%gs) occured while waiting to get exclusive '
2070 'access to: %s' % (timeout, fn))
2072 if timewarn is not None and tnow - t0 > timewarn:
2073 logger.warning(
2074 'Waiting since %gs to get exclusive access to: %s' % (
2075 timewarn, fn))
2077 timewarn *= 2
2079 time.sleep(0.01)
2080 else:
2081 raise
2084def delete_lockfile(fn):
2085 os.unlink(fn)
2088class Lockfile(Exception):
2090 def __init__(self, path, timeout=5, timewarn=10.):
2091 self._path = path
2092 self._timeout = timeout
2093 self._timewarn = timewarn
2095 def __enter__(self):
2096 create_lockfile(
2097 self._path, timeout=self._timeout, timewarn=self._timewarn)
2098 return None
2100 def __exit__(self, type, value, traceback):
2101 delete_lockfile(self._path)
2104class SoleError(Exception):
2105 '''
2106 Exception raised by objects of type :py:class:`Sole`, when an concurrent
2107 instance is running.
2108 '''
2110 pass
2113class Sole(object):
2114 '''
2115 Use POSIX advisory file locking to ensure that only a single instance of a
2116 program is running.
2118 :param pid_path: path to lockfile to be used
2120 Usage::
2122 from pyrocko.util import Sole, SoleError, setup_logging
2123 import os
2125 setup_logging('my_program')
2127 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock')
2128 try:
2129 sole = Sole(pid_path)
2131 except SoleError, e:
2132 logger.fatal( str(e) )
2133 sys.exit(1)
2134 '''
2136 def __init__(self, pid_path):
2137 self._pid_path = pid_path
2138 self._other_running = False
2139 ensuredirs(self._pid_path)
2140 self._lockfile = None
2141 self._os = os
2143 if not fcntl:
2144 raise SoleError(
2145 'Python standard library module "fcntl" not available on '
2146 'this platform.')
2148 self._fcntl = fcntl
2150 try:
2151 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY)
2152 except Exception:
2153 raise SoleError(
2154 'Cannot open lockfile (path = %s)' % self._pid_path)
2156 try:
2157 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB)
2159 except IOError:
2160 self._other_running = True
2161 try:
2162 f = open(self._pid_path, 'r')
2163 pid = f.read().strip()
2164 f.close()
2165 except Exception:
2166 pid = '?'
2168 raise SoleError('Other instance is running (pid = %s)' % pid)
2170 try:
2171 os.ftruncate(self._lockfile, 0)
2172 os.write(self._lockfile, '%i\n' % os.getpid())
2173 os.fsync(self._lockfile)
2175 except Exception:
2176 # the pid is only stored for user information, so this is allowed
2177 # to fail
2178 pass
2180 def __del__(self):
2181 if not self._other_running:
2182 if self._lockfile is not None:
2183 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN)
2184 self._os.close(self._lockfile)
2185 try:
2186 self._os.unlink(self._pid_path)
2187 except Exception:
2188 pass
2191re_escapequotes = re.compile(r"(['\\])")
2194def escapequotes(s):
2195 return re_escapequotes.sub(r"\\\1", s)
2198class TableWriter(object):
2199 '''
2200 Write table of space separated values to a file.
2202 :param f: file like object
2204 Strings containing spaces are quoted on output.
2205 '''
2207 def __init__(self, f):
2208 self._f = f
2210 def writerow(self, row, minfieldwidths=None):
2212 '''
2213 Write one row of values to underlying file.
2215 :param row: iterable of values
2216 :param minfieldwidths: minimum field widths for the values
2218 Each value in in ``row`` is converted to a string and optionally padded
2219 with blanks. The resulting strings are output separated with blanks. If
2220 any values given are strings and if they contain whitespace, they are
2221 quoted with single quotes, and any internal single quotes are
2222 backslash-escaped.
2223 '''
2225 out = []
2227 for i, x in enumerate(row):
2228 w = 0
2229 if minfieldwidths and i < len(minfieldwidths):
2230 w = minfieldwidths[i]
2232 if isinstance(x, str):
2233 if re.search(r"\s|'", x):
2234 x = "'%s'" % escapequotes(x)
2236 x = x.ljust(w)
2237 else:
2238 x = str(x).rjust(w)
2240 out.append(x)
2242 self._f.write(' '.join(out).rstrip() + '\n')
2245class TableReader(object):
2247 '''
2248 Read table of space separated values from a file.
2250 :param f: file-like object
2252 This uses Pythons shlex module to tokenize lines. Should deal correctly
2253 with quoted strings.
2254 '''
2256 def __init__(self, f):
2257 self._f = f
2258 self.eof = False
2260 def readrow(self):
2261 '''
2262 Read one row from the underlying file, tokenize it with shlex.
2264 :returns: tokenized line as a list of strings.
2265 '''
2267 line = self._f.readline()
2268 if not line:
2269 self.eof = True
2270 return []
2271 line.strip()
2272 if line.startswith('#'):
2273 return []
2275 return qsplit(line)
2278def gform(number, significant_digits=3):
2279 '''
2280 Pretty print floating point numbers.
2282 Align floating point numbers at the decimal dot.
2284 ::
2286 | -d.dde+xxx|
2287 | -d.dde+xx |
2288 |-ddd. |
2289 | -dd.d |
2290 | -d.dd |
2291 | -0.ddd |
2292 | -0.0ddd |
2293 | -0.00ddd |
2294 | -d.dde-xx |
2295 | -d.dde-xxx|
2296 | nan|
2299 The formatted string has length ``significant_digits * 2 + 6``.
2300 '''
2302 no_exp_range = (pow(10., -1),
2303 pow(10., significant_digits))
2304 width = significant_digits+significant_digits-1+1+1+5
2306 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.:
2307 s = ('%#.*g' % (significant_digits, number)).rstrip('0')
2308 else:
2309 s = '%.*E' % (significant_digits-1, number)
2310 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width)
2311 if s.strip().lower() == 'nan':
2312 s = 'nan'.rjust(width)
2313 return s
2316def human_bytesize(value):
2318 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split()
2320 if value == 1:
2321 return '1 Byte'
2323 for i, ext in enumerate(exts):
2324 x = float(value) / 1000**i
2325 if round(x) < 10. and not value < 1000:
2326 return '%.1f %s' % (x, ext)
2327 if round(x) < 1000.:
2328 return '%.0f %s' % (x, ext)
2330 return '%i Bytes' % value
2333re_compatibility = re.compile(
2334 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' +
2335 r'(dummy|poel|qseis|qssp))\.'
2336)
2339def pf_is_old(fn):
2340 oldstyle = False
2341 with open(fn, 'r') as f:
2342 for line in f:
2343 if re_compatibility.search(line):
2344 oldstyle = True
2346 return oldstyle
2349def pf_upgrade(fn):
2350 need = pf_is_old(fn)
2351 if need:
2352 fn_temp = fn + '.temp'
2354 with open(fn, 'r') as fin:
2355 with open(fn_temp, 'w') as fout:
2356 for line in fin:
2357 line = re_compatibility.sub('!pf.', line)
2358 fout.write(line)
2360 os.rename(fn_temp, fn)
2362 return need
2365def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'):
2366 '''
2367 Extract leap second information from tzdata.
2369 Based on example at http://stackoverflow.com/questions/19332902/\
2370 extract-historic-leap-seconds-from-tzdata
2372 See also 'man 5 tzfile'.
2373 '''
2375 from struct import unpack, calcsize
2376 out = []
2377 with open(tzfile, 'rb') as f:
2378 # read header
2379 fmt = '>4s c 15x 6l'
2380 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt,
2381 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt)))
2382 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file'
2384 # skip over some uninteresting data
2385 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict(
2386 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt)
2387 f.read(calcsize(fmt))
2389 # read leap-seconds
2390 fmt = '>2l'
2391 for i in range(leapcnt):
2392 tleap, nleap = unpack(fmt, f.read(calcsize(fmt)))
2393 out.append((tleap-nleap+1, nleap))
2395 return out
2398class LeapSecondsError(Exception):
2399 pass
2402class LeapSecondsOutdated(LeapSecondsError):
2403 pass
2406class InvalidLeapSecondsFile(LeapSecondsOutdated):
2407 pass
2410def parse_leap_seconds_list(fn):
2411 data = []
2412 texpires = None
2413 try:
2414 t0 = int(round(str_to_time('1900-01-01 00:00:00')))
2415 except TimeStrError:
2416 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800
2418 tnow = int(round(time.time()))
2420 if not op.exists(fn):
2421 raise LeapSecondsOutdated('no leap seconds file found')
2423 try:
2424 with open(fn, 'rb') as f:
2425 for line in f:
2426 if line.strip().startswith(b'<!DOCTYPE'):
2427 raise InvalidLeapSecondsFile('invalid leap seconds file')
2429 if line.startswith(b'#@'):
2430 texpires = int(line.split()[1]) + t0
2431 elif line.startswith(b'#') or len(line) < 5:
2432 pass
2433 else:
2434 toks = line.split()
2435 t = int(toks[0]) + t0
2436 nleap = int(toks[1]) - 10
2437 data.append((t, nleap))
2439 except IOError:
2440 raise LeapSecondsError('cannot read leap seconds file %s' % fn)
2442 if texpires is None or tnow > texpires:
2443 raise LeapSecondsOutdated('leap seconds list is outdated')
2445 return data
2448def read_leap_seconds2():
2449 from pyrocko import config
2450 conf = config.config()
2451 fn = conf.leapseconds_path
2452 url = conf.leapseconds_url
2453 # check for outdated default URL
2454 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list':
2455 url = 'https://www.ietf.org/timezones/data/leap-seconds.list'
2456 logger.info(
2457 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2459 for i in range(3):
2460 try:
2461 return parse_leap_seconds_list(fn)
2463 except LeapSecondsOutdated:
2464 try:
2465 logger.info('updating leap seconds list...')
2466 download_file(url, fn)
2468 except Exception as e:
2469 raise LeapSecondsError(
2470 'cannot download leap seconds list from %s to %s (%s)'
2471 % (url, fn, e))
2473 raise LeapSecondsError('Could not retrieve/read leap seconds file.')
2476def gps_utc_offset(t_utc):
2477 '''
2478 Time offset t_gps - t_utc for a given t_utc.
2479 '''
2480 ls = read_leap_seconds2()
2481 i = 0
2482 if t_utc < ls[0][0]:
2483 return ls[0][1] - 1 - 9
2485 while i < len(ls) - 1:
2486 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]:
2487 return ls[i][1] - 9
2488 i += 1
2490 return ls[-1][1] - 9
2493def utc_gps_offset(t_gps):
2494 '''
2495 Time offset t_utc - t_gps for a given t_gps.
2496 '''
2497 ls = read_leap_seconds2()
2499 if t_gps < ls[0][0] + ls[0][1] - 9:
2500 return - (ls[0][1] - 1 - 9)
2502 i = 0
2503 while i < len(ls) - 1:
2504 if ls[i][0] + ls[i][1] - 9 <= t_gps \
2505 and t_gps < ls[i+1][0] + ls[i+1][1] - 9:
2506 return - (ls[i][1] - 9)
2507 i += 1
2509 return - (ls[-1][1] - 9)
2512def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'):
2513 import itertools
2514 import glob
2515 from pyrocko.io.io_common import FileLoadError
2517 def iload_filename(filename, **kwargs):
2518 try:
2519 with open(filename, 'rb') as f:
2520 for cr in iload_fh(f, **kwargs):
2521 yield cr
2523 except FileLoadError as e:
2524 e.set_context('filename', filename)
2525 raise
2527 def iload_dirname(dirname, **kwargs):
2528 for entry in os.listdir(dirname):
2529 fpath = op.join(dirname, entry)
2530 if op.isfile(fpath):
2531 for cr in iload_filename(fpath, **kwargs):
2532 yield cr
2534 def iload_glob(pattern, **kwargs):
2536 for fn in glob.iglob(pattern):
2537 for cr in iload_filename(fn, **kwargs):
2538 yield cr
2540 def iload(source, **kwargs):
2541 if isinstance(source, str):
2542 if op.isdir(source):
2543 return iload_dirname(source, **kwargs)
2544 elif op.isfile(source):
2545 return iload_filename(source, **kwargs)
2546 else:
2547 return iload_glob(source, **kwargs)
2549 elif hasattr(source, 'read'):
2550 return iload_fh(source, **kwargs)
2551 else:
2552 return itertools.chain.from_iterable(
2553 iload(subsource, **kwargs) for subsource in source)
2555 iload_filename.__doc__ = '''
2556 Read %s information from named file.
2557 ''' % doc_fmt
2559 iload_dirname.__doc__ = '''
2560 Read %s information from directory of %s files.
2561 ''' % (doc_fmt, doc_fmt)
2563 iload_glob.__doc__ = '''
2564 Read %s information from files matching a glob pattern.
2565 ''' % doc_fmt
2567 iload.__doc__ = '''
2568 Load %s information from given source(s)
2570 The ``source`` can be specified as the name of a %s file, the name of a
2571 directory containing %s files, a glob pattern of %s files, an open
2572 filehandle or an iterator yielding any of the forementioned sources.
2574 This function behaves as a generator yielding %s objects.
2575 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects)
2577 for f in iload_filename, iload_dirname, iload_glob, iload:
2578 f.__module__ = iload_fh.__module__
2580 return iload_filename, iload_dirname, iload_glob, iload
2583class Inconsistency(Exception):
2584 pass
2587def consistency_check(list_of_tuples, message='values differ:'):
2588 '''
2589 Check for inconsistencies.
2591 Given a list of tuples, check that all tuple elements except for first one
2592 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be
2593 valid because the coordinates at the two channels are the same.
2594 '''
2596 if len(list_of_tuples) >= 2:
2597 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]):
2598 raise Inconsistency('%s\n' % message + '\n'.join(
2599 ' %s: %s' % (t[0], ', '.join('%g' % x for x in t[1:]))
2600 for t in list_of_tuples))
2603class defaultzerodict(dict):
2604 def __missing__(self, k):
2605 return 0
2608def mostfrequent(x):
2609 c = defaultzerodict()
2610 for e in x:
2611 c[e] += 1
2613 return sorted(list(c.keys()), key=lambda k: c[k])[-1]
2616def consistency_merge(list_of_tuples,
2617 message='values differ:',
2618 error='raise',
2619 merge=mostfrequent):
2621 assert error in ('raise', 'warn', 'ignore')
2623 if len(list_of_tuples) == 0:
2624 raise Exception('cannot merge empty sequence')
2626 try:
2627 consistency_check(list_of_tuples, message)
2628 return list_of_tuples[0][1:]
2629 except Inconsistency as e:
2630 if error == 'raise':
2631 raise
2633 elif error == 'warn':
2634 logger.warning(str(e))
2636 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]])
2639def short_to_list(nmax, it):
2640 import itertools
2642 if isinstance(it, list):
2643 return it
2645 li = []
2646 for i in range(nmax+1):
2647 try:
2648 li.append(next(it))
2649 except StopIteration:
2650 return li
2652 return itertools.chain(li, it)
2655def parse_md(f):
2656 try:
2657 with open(op.join(
2658 op.dirname(op.abspath(f)),
2659 'README.md'), 'r') as readme:
2660 mdstr = readme.read()
2661 except IOError as e:
2662 return 'Failed to get README.md: %s' % e
2664 # Remve the title
2665 mdstr = re.sub(r'^# .*\n?', '', mdstr)
2666 # Append sphinx reference to `pyrocko.` modules
2667 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr)
2668 # Convert Subsections to toc-less rubrics
2669 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr)
2670 return mdstr
2673def mpl_show(plt):
2674 import matplotlib
2675 if matplotlib.get_backend().lower() == 'agg':
2676 logger.warning('Cannot show() when using matplotlib "agg" backend')
2677 else:
2678 plt.show()
2681g_re_qsplit = re.compile(
2682 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)')
2685g_re_trivial = re.compile(r'\A[^\'"\s]+\Z')
2686g_re_escape_s = re.compile(r'([\\\'])')
2687g_re_unescape_s = re.compile(r'\\([\\\'])')
2688g_re_escape_d = re.compile(r'([\\"])')
2689g_re_unescape_d = re.compile(r'\\([\\"])')
2692def escape_s(s):
2693 '''
2694 Backslash-escape single-quotes and backslashes.
2696 Example: ``Jack's`` => ``Jack\\'s``
2698 '''
2699 return g_re_escape_s.sub(r'\\\1', s)
2702def unescape_s(s):
2703 '''
2704 Unescape backslash-escaped single-quotes and backslashes.
2706 Example: ``Jack\\'s`` => ``Jack's``
2707 '''
2708 return g_re_unescape_s.sub(r'\1', s)
2711def escape_d(s):
2712 '''
2713 Backslash-escape double-quotes and backslashes.
2715 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"``
2716 '''
2717 return g_re_escape_d.sub(r'\\\1', s)
2720def unescape_d(s):
2721 '''
2722 Unescape backslash-escaped double-quotes and backslashes.
2724 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"``
2725 '''
2726 return g_re_unescape_d.sub(r'\1', s)
2729def qjoin_s(it):
2730 '''
2731 Join sequence of strings into a line, single-quoting non-trivial strings.
2733 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"``
2734 '''
2735 return ' '.join(
2736 w if g_re_trivial.search(w) else "'%s'" % escape_s(w) for w in it)
2739def qjoin_d(it):
2740 '''
2741 Join sequence of strings into a line, double-quoting non-trivial strings.
2743 Example: ``['55', 'Pete "The Robot" Smith']`` =>
2744 ``'55' "Pete \\\\"The Robot\\\\" Smith"'``
2745 '''
2746 return ' '.join(
2747 w if g_re_trivial.search(w) else '"%s"' % escape_d(w) for w in it)
2750def qsplit(s):
2751 '''
2752 Split line into list of strings, allowing for quoted strings.
2754 Example: ``"55 'Sparrow\\\\'s Island'"`` =>
2755 ``["55", "Sparrow's Island"]``,
2756 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` =>
2757 ``['55', 'Pete "The Robot" Smith']``
2758 '''
2759 return [
2760 (unescape_d(x[0]) or unescape_s(x[1]) or x[2])
2761 for x in g_re_qsplit.findall(s)]
2764g_have_warned_threadpoolctl = False
2767class threadpool_limits_dummy(object):
2769 def __init__(self, *args, **kwargs):
2770 pass
2772 def __enter__(self):
2773 global g_have_warned_threadpoolctl
2775 if not g_have_warned_threadpoolctl:
2776 logger.warning(
2777 'Cannot control number of BLAS threads because '
2778 '`threadpoolctl` module is not available. You may want to '
2779 'install `threadpoolctl`.')
2781 g_have_warned_threadpoolctl = True
2783 return self
2785 def __exit__(self, type, value, traceback):
2786 pass
2789def get_threadpool_limits():
2790 '''
2791 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail.
2792 '''
2794 try:
2795 from threadpoolctl import threadpool_limits
2796 return threadpool_limits
2798 except ImportError:
2799 return threadpool_limits_dummy