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 include=None,
1770 exclude=None,
1771 selector=None,
1772 show_progress=True,
1773 pass_through=None):
1775 '''
1776 Recursively select files (generator variant).
1778 See :py:func:`select_files`.
1779 '''
1781 if show_progress:
1782 progress_beg('selecting files...')
1784 ngood = 0
1785 check_include = None
1786 if include is not None:
1787 rinclude = re.compile(include)
1789 def check_include(path):
1790 m = rinclude.search(path)
1791 if not m:
1792 return False
1794 if selector is None:
1795 return True
1797 infos = Anon(**m.groupdict())
1798 return selector(infos)
1800 check_exclude = None
1801 if exclude is not None:
1802 rexclude = re.compile(exclude)
1804 def check_exclude(path):
1805 return not bool(rexclude.search(path))
1807 if check_include and check_exclude:
1809 def check(path):
1810 return check_include(path) and check_exclude(path)
1812 elif check_include:
1813 check = check_include
1815 elif check_exclude:
1816 check = check_exclude
1818 else:
1819 check = None
1821 if isinstance(paths, str):
1822 paths = [paths]
1824 for path in paths:
1825 if pass_through and pass_through(path):
1826 if check is None or check(path):
1827 yield path
1829 elif os.path.isdir(path):
1830 for (dirpath, dirnames, filenames) in os.walk(path):
1831 dirnames.sort()
1832 filenames.sort()
1833 for filename in filenames:
1834 path = op.join(dirpath, filename)
1835 if check is None or check(path):
1836 yield os.path.abspath(path)
1837 ngood += 1
1838 else:
1839 if check is None or check(path):
1840 yield os.path.abspath(path)
1841 ngood += 1
1843 if show_progress:
1844 progress_end('%i file%s selected.' % (ngood, plural_s(ngood)))
1847def select_files(
1848 paths, include=None, exclude=None, selector=None, show_progress=True):
1850 '''
1851 Recursively select files.
1853 :param paths: entry path names
1854 :param include: pattern for conditional inclusion
1855 :param exclude: pattern for conditional exclusion
1856 :param selector: callback for conditional inclusion
1857 :param show_progress: if True, indicate start and stop of processing
1858 :returns: list of path names
1860 Recursively finds all files under given entry points ``paths``. If
1861 parameter ``include`` is a regular expression, only files with matching
1862 path names are included. If additionally parameter ``selector`` is given a
1863 callback function, only files for which the callback returns ``True`` are
1864 included. The callback should take a single argument. The callback is
1865 called with a single argument, an object, having as attributes, any named
1866 groups given in ``include``.
1868 Examples
1870 To find all files ending in ``'.mseed'`` or ``'.msd'``::
1872 select_files(paths,
1873 include=r'\\.(mseed|msd)$')
1875 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for
1876 the year::
1878 select_files(paths,
1879 include=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$',
1880 selector=(lambda x: int(x.year) == 2009))
1881 '''
1882 return list(iter_select_files(
1883 paths, include, exclude, selector, show_progress))
1886def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
1887 '''
1888 Convert positive integer to a base36 string.
1889 '''
1891 if not isinstance(number, (int, long)):
1892 raise TypeError('number must be an integer')
1893 if number < 0:
1894 raise ValueError('number must be positive')
1896 # Special case for small numbers
1897 if number < 36:
1898 return alphabet[number]
1900 base36 = ''
1901 while number != 0:
1902 number, i = divmod(number, 36)
1903 base36 = alphabet[i] + base36
1905 return base36
1908def base36decode(number):
1909 '''
1910 Decode base36 endcoded positive integer.
1911 '''
1913 return int(number, 36)
1916class UnpackError(Exception):
1917 '''
1918 Exception raised when :py:func:`unpack_fixed` encounters an error.
1919 '''
1921 pass
1924ruler = ''.join(['%-10i' % i for i in range(8)]) \
1925 + '\n' + '0123456789' * 8 + '\n'
1928def unpack_fixed(format, line, *callargs):
1929 '''
1930 Unpack fixed format string, as produced by many fortran codes.
1932 :param format: format specification
1933 :param line: string to be processed
1934 :param callargs: callbacks for callback fields in the format
1936 The format is described by a string of comma-separated fields. Each field
1937 is defined by a character for the field type followed by the field width. A
1938 questionmark may be appended to the field description to allow the argument
1939 to be optional (The data string is then allowed to be filled with blanks
1940 and ``None`` is returned in this case).
1942 The following field types are available:
1944 ==== ================================================================
1945 Type Description
1946 ==== ================================================================
1947 A string (full field width is extracted)
1948 a string (whitespace at the beginning and the end is removed)
1949 i integer value
1950 f floating point value
1951 @ special type, a callback must be given for the conversion
1952 x special field type to skip parts of the string
1953 ==== ================================================================
1954 '''
1956 ipos = 0
1957 values = []
1958 icall = 0
1959 for form in format.split(','):
1960 form = form.strip()
1961 optional = form[-1] == '?'
1962 form = form.rstrip('?')
1963 typ = form[0]
1964 ln = int(form[1:])
1965 s = line[ipos:ipos+ln]
1966 cast = {
1967 'x': None,
1968 'A': str,
1969 'a': lambda x: x.strip(),
1970 'i': int,
1971 'f': float,
1972 '@': 'extra'}[typ]
1974 if cast == 'extra':
1975 cast = callargs[icall]
1976 icall += 1
1978 if cast is not None:
1979 if optional and s.strip() == '':
1980 values.append(None)
1981 else:
1982 try:
1983 values.append(cast(s))
1984 except Exception:
1985 mark = [' '] * 80
1986 mark[ipos:ipos+ln] = ['^'] * ln
1987 mark = ''.join(mark)
1988 raise UnpackError(
1989 'Invalid cast to type "%s" at position [%i:%i] of '
1990 'line: \n%s%s\n%s' % (
1991 typ, ipos, ipos+ln, ruler, line.rstrip(), mark))
1993 ipos += ln
1995 return values
1998_pattern_cache = {}
2001def _nslc_pattern(pattern):
2002 if pattern not in _pattern_cache:
2003 rpattern = re.compile(fnmatch.translate(pattern), re.I)
2004 _pattern_cache[pattern] = rpattern
2005 else:
2006 rpattern = _pattern_cache[pattern]
2008 return rpattern
2011def match_nslc(patterns, nslc):
2012 '''
2013 Match network-station-location-channel code against pattern or list of
2014 patterns.
2016 :param patterns: pattern or list of patterns
2017 :param nslc: tuple with (network, station, location, channel) as strings
2019 :returns: ``True`` if the pattern matches or if any of the given patterns
2020 match; or ``False``.
2022 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq].
2024 Example::
2026 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True
2027 '''
2029 if isinstance(patterns, str):
2030 patterns = [patterns]
2032 if not isinstance(nslc, str):
2033 s = '.'.join(nslc)
2034 else:
2035 s = nslc
2037 for pattern in patterns:
2038 if _nslc_pattern(pattern).match(s):
2039 return True
2041 return False
2044def match_nslcs(patterns, nslcs):
2045 '''
2046 Get network-station-location-channel codes that match given pattern or any
2047 of several given patterns.
2049 :param patterns: pattern or list of patterns
2050 :param nslcs: list of (network, station, location, channel) tuples
2052 See also :py:func:`match_nslc`
2053 '''
2055 matching = []
2056 for nslc in nslcs:
2057 if match_nslc(patterns, nslc):
2058 matching.append(nslc)
2060 return matching
2063class Timeout(Exception):
2064 pass
2067def create_lockfile(fn, timeout=None, timewarn=10.):
2068 t0 = time.time()
2070 while True:
2071 try:
2072 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL)
2073 os.close(f)
2074 return
2076 except OSError as e:
2077 if e.errno == errno.EEXIST:
2078 tnow = time.time()
2080 if timeout is not None and tnow - t0 > timeout:
2081 raise Timeout(
2082 'Timeout (%gs) occured while waiting to get exclusive '
2083 'access to: %s' % (timeout, fn))
2085 if timewarn is not None and tnow - t0 > timewarn:
2086 logger.warning(
2087 'Waiting since %gs to get exclusive access to: %s' % (
2088 timewarn, fn))
2090 timewarn *= 2
2092 time.sleep(0.01)
2093 else:
2094 raise
2097def delete_lockfile(fn):
2098 os.unlink(fn)
2101class Lockfile(Exception):
2103 def __init__(self, path, timeout=5, timewarn=10.):
2104 self._path = path
2105 self._timeout = timeout
2106 self._timewarn = timewarn
2108 def __enter__(self):
2109 create_lockfile(
2110 self._path, timeout=self._timeout, timewarn=self._timewarn)
2111 return None
2113 def __exit__(self, type, value, traceback):
2114 delete_lockfile(self._path)
2117class SoleError(Exception):
2118 '''
2119 Exception raised by objects of type :py:class:`Sole`, when an concurrent
2120 instance is running.
2121 '''
2123 pass
2126class Sole(object):
2127 '''
2128 Use POSIX advisory file locking to ensure that only a single instance of a
2129 program is running.
2131 :param pid_path: path to lockfile to be used
2133 Usage::
2135 from pyrocko.util import Sole, SoleError, setup_logging
2136 import os
2138 setup_logging('my_program')
2140 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock')
2141 try:
2142 sole = Sole(pid_path)
2144 except SoleError, e:
2145 logger.fatal( str(e) )
2146 sys.exit(1)
2147 '''
2149 def __init__(self, pid_path):
2150 self._pid_path = pid_path
2151 self._other_running = False
2152 ensuredirs(self._pid_path)
2153 self._lockfile = None
2154 self._os = os
2156 if not fcntl:
2157 raise SoleError(
2158 'Python standard library module "fcntl" not available on '
2159 'this platform.')
2161 self._fcntl = fcntl
2163 try:
2164 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY)
2165 except Exception:
2166 raise SoleError(
2167 'Cannot open lockfile (path = %s)' % self._pid_path)
2169 try:
2170 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB)
2172 except IOError:
2173 self._other_running = True
2174 try:
2175 f = open(self._pid_path, 'r')
2176 pid = f.read().strip()
2177 f.close()
2178 except Exception:
2179 pid = '?'
2181 raise SoleError('Other instance is running (pid = %s)' % pid)
2183 try:
2184 os.ftruncate(self._lockfile, 0)
2185 os.write(self._lockfile, '%i\n' % os.getpid())
2186 os.fsync(self._lockfile)
2188 except Exception:
2189 # the pid is only stored for user information, so this is allowed
2190 # to fail
2191 pass
2193 def __del__(self):
2194 if not self._other_running:
2195 if self._lockfile is not None:
2196 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN)
2197 self._os.close(self._lockfile)
2198 try:
2199 self._os.unlink(self._pid_path)
2200 except Exception:
2201 pass
2204re_escapequotes = re.compile(r"(['\\])")
2207def escapequotes(s):
2208 return re_escapequotes.sub(r"\\\1", s)
2211class TableWriter(object):
2212 '''
2213 Write table of space separated values to a file.
2215 :param f: file like object
2217 Strings containing spaces are quoted on output.
2218 '''
2220 def __init__(self, f):
2221 self._f = f
2223 def writerow(self, row, minfieldwidths=None):
2225 '''
2226 Write one row of values to underlying file.
2228 :param row: iterable of values
2229 :param minfieldwidths: minimum field widths for the values
2231 Each value in in ``row`` is converted to a string and optionally padded
2232 with blanks. The resulting strings are output separated with blanks. If
2233 any values given are strings and if they contain whitespace, they are
2234 quoted with single quotes, and any internal single quotes are
2235 backslash-escaped.
2236 '''
2238 out = []
2240 for i, x in enumerate(row):
2241 w = 0
2242 if minfieldwidths and i < len(minfieldwidths):
2243 w = minfieldwidths[i]
2245 if isinstance(x, str):
2246 if re.search(r"\s|'", x):
2247 x = "'%s'" % escapequotes(x)
2249 x = x.ljust(w)
2250 else:
2251 x = str(x).rjust(w)
2253 out.append(x)
2255 self._f.write(' '.join(out).rstrip() + '\n')
2258class TableReader(object):
2260 '''
2261 Read table of space separated values from a file.
2263 :param f: file-like object
2265 This uses Pythons shlex module to tokenize lines. Should deal correctly
2266 with quoted strings.
2267 '''
2269 def __init__(self, f):
2270 self._f = f
2271 self.eof = False
2273 def readrow(self):
2274 '''
2275 Read one row from the underlying file, tokenize it with shlex.
2277 :returns: tokenized line as a list of strings.
2278 '''
2280 line = self._f.readline()
2281 if not line:
2282 self.eof = True
2283 return []
2284 line.strip()
2285 if line.startswith('#'):
2286 return []
2288 return qsplit(line)
2291def gform(number, significant_digits=3):
2292 '''
2293 Pretty print floating point numbers.
2295 Align floating point numbers at the decimal dot.
2297 ::
2299 | -d.dde+xxx|
2300 | -d.dde+xx |
2301 |-ddd. |
2302 | -dd.d |
2303 | -d.dd |
2304 | -0.ddd |
2305 | -0.0ddd |
2306 | -0.00ddd |
2307 | -d.dde-xx |
2308 | -d.dde-xxx|
2309 | nan|
2312 The formatted string has length ``significant_digits * 2 + 6``.
2313 '''
2315 no_exp_range = (pow(10., -1),
2316 pow(10., significant_digits))
2317 width = significant_digits+significant_digits-1+1+1+5
2319 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.:
2320 s = ('%#.*g' % (significant_digits, number)).rstrip('0')
2321 else:
2322 s = '%.*E' % (significant_digits-1, number)
2323 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width)
2324 if s.strip().lower() == 'nan':
2325 s = 'nan'.rjust(width)
2326 return s
2329def human_bytesize(value):
2331 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split()
2333 if value == 1:
2334 return '1 Byte'
2336 for i, ext in enumerate(exts):
2337 x = float(value) / 1000**i
2338 if round(x) < 10. and not value < 1000:
2339 return '%.1f %s' % (x, ext)
2340 if round(x) < 1000.:
2341 return '%.0f %s' % (x, ext)
2343 return '%i Bytes' % value
2346re_compatibility = re.compile(
2347 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' +
2348 r'(dummy|poel|qseis|qssp))\.'
2349)
2352def pf_is_old(fn):
2353 oldstyle = False
2354 with open(fn, 'r') as f:
2355 for line in f:
2356 if re_compatibility.search(line):
2357 oldstyle = True
2359 return oldstyle
2362def pf_upgrade(fn):
2363 need = pf_is_old(fn)
2364 if need:
2365 fn_temp = fn + '.temp'
2367 with open(fn, 'r') as fin:
2368 with open(fn_temp, 'w') as fout:
2369 for line in fin:
2370 line = re_compatibility.sub('!pf.', line)
2371 fout.write(line)
2373 os.rename(fn_temp, fn)
2375 return need
2378def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'):
2379 '''
2380 Extract leap second information from tzdata.
2382 Based on example at http://stackoverflow.com/questions/19332902/\
2383 extract-historic-leap-seconds-from-tzdata
2385 See also 'man 5 tzfile'.
2386 '''
2388 from struct import unpack, calcsize
2389 out = []
2390 with open(tzfile, 'rb') as f:
2391 # read header
2392 fmt = '>4s c 15x 6l'
2393 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt,
2394 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt)))
2395 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file'
2397 # skip over some uninteresting data
2398 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict(
2399 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt)
2400 f.read(calcsize(fmt))
2402 # read leap-seconds
2403 fmt = '>2l'
2404 for i in range(leapcnt):
2405 tleap, nleap = unpack(fmt, f.read(calcsize(fmt)))
2406 out.append((tleap-nleap+1, nleap))
2408 return out
2411class LeapSecondsError(Exception):
2412 pass
2415class LeapSecondsOutdated(LeapSecondsError):
2416 pass
2419class InvalidLeapSecondsFile(LeapSecondsOutdated):
2420 pass
2423def parse_leap_seconds_list(fn):
2424 data = []
2425 texpires = None
2426 try:
2427 t0 = int(round(str_to_time('1900-01-01 00:00:00')))
2428 except TimeStrError:
2429 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800
2431 tnow = int(round(time.time()))
2433 if not op.exists(fn):
2434 raise LeapSecondsOutdated('no leap seconds file found')
2436 try:
2437 with open(fn, 'rb') as f:
2438 for line in f:
2439 if line.strip().startswith(b'<!DOCTYPE'):
2440 raise InvalidLeapSecondsFile('invalid leap seconds file')
2442 if line.startswith(b'#@'):
2443 texpires = int(line.split()[1]) + t0
2444 elif line.startswith(b'#') or len(line) < 5:
2445 pass
2446 else:
2447 toks = line.split()
2448 t = int(toks[0]) + t0
2449 nleap = int(toks[1]) - 10
2450 data.append((t, nleap))
2452 except IOError:
2453 raise LeapSecondsError('cannot read leap seconds file %s' % fn)
2455 if texpires is None or tnow > texpires:
2456 raise LeapSecondsOutdated('leap seconds list is outdated')
2458 return data
2461def read_leap_seconds2():
2462 from pyrocko import config
2463 conf = config.config()
2464 fn = conf.leapseconds_path
2465 url = conf.leapseconds_url
2466 # check for outdated default URL
2467 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list':
2468 url = 'https://www.ietf.org/timezones/data/leap-seconds.list'
2469 logger.info(
2470 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2472 for i in range(3):
2473 try:
2474 return parse_leap_seconds_list(fn)
2476 except LeapSecondsOutdated:
2477 try:
2478 logger.info('updating leap seconds list...')
2479 download_file(url, fn)
2481 except Exception as e:
2482 raise LeapSecondsError(
2483 'cannot download leap seconds list from %s to %s (%s)'
2484 % (url, fn, e))
2486 raise LeapSecondsError('Could not retrieve/read leap seconds file.')
2489def gps_utc_offset(t_utc):
2490 '''
2491 Time offset t_gps - t_utc for a given t_utc.
2492 '''
2493 ls = read_leap_seconds2()
2494 i = 0
2495 if t_utc < ls[0][0]:
2496 return ls[0][1] - 1 - 9
2498 while i < len(ls) - 1:
2499 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]:
2500 return ls[i][1] - 9
2501 i += 1
2503 return ls[-1][1] - 9
2506def utc_gps_offset(t_gps):
2507 '''
2508 Time offset t_utc - t_gps for a given t_gps.
2509 '''
2510 ls = read_leap_seconds2()
2512 if t_gps < ls[0][0] + ls[0][1] - 9:
2513 return - (ls[0][1] - 1 - 9)
2515 i = 0
2516 while i < len(ls) - 1:
2517 if ls[i][0] + ls[i][1] - 9 <= t_gps \
2518 and t_gps < ls[i+1][0] + ls[i+1][1] - 9:
2519 return - (ls[i][1] - 9)
2520 i += 1
2522 return - (ls[-1][1] - 9)
2525def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'):
2526 import itertools
2527 import glob
2528 from pyrocko.io.io_common import FileLoadError
2530 def iload_filename(filename, **kwargs):
2531 try:
2532 with open(filename, 'rb') as f:
2533 for cr in iload_fh(f, **kwargs):
2534 yield cr
2536 except FileLoadError as e:
2537 e.set_context('filename', filename)
2538 raise
2540 def iload_dirname(dirname, **kwargs):
2541 for entry in os.listdir(dirname):
2542 fpath = op.join(dirname, entry)
2543 if op.isfile(fpath):
2544 for cr in iload_filename(fpath, **kwargs):
2545 yield cr
2547 def iload_glob(pattern, **kwargs):
2549 for fn in glob.iglob(pattern):
2550 for cr in iload_filename(fn, **kwargs):
2551 yield cr
2553 def iload(source, **kwargs):
2554 if isinstance(source, str):
2555 if op.isdir(source):
2556 return iload_dirname(source, **kwargs)
2557 elif op.isfile(source):
2558 return iload_filename(source, **kwargs)
2559 else:
2560 return iload_glob(source, **kwargs)
2562 elif hasattr(source, 'read'):
2563 return iload_fh(source, **kwargs)
2564 else:
2565 return itertools.chain.from_iterable(
2566 iload(subsource, **kwargs) for subsource in source)
2568 iload_filename.__doc__ = '''
2569 Read %s information from named file.
2570 ''' % doc_fmt
2572 iload_dirname.__doc__ = '''
2573 Read %s information from directory of %s files.
2574 ''' % (doc_fmt, doc_fmt)
2576 iload_glob.__doc__ = '''
2577 Read %s information from files matching a glob pattern.
2578 ''' % doc_fmt
2580 iload.__doc__ = '''
2581 Load %s information from given source(s)
2583 The ``source`` can be specified as the name of a %s file, the name of a
2584 directory containing %s files, a glob pattern of %s files, an open
2585 filehandle or an iterator yielding any of the forementioned sources.
2587 This function behaves as a generator yielding %s objects.
2588 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects)
2590 for f in iload_filename, iload_dirname, iload_glob, iload:
2591 f.__module__ = iload_fh.__module__
2593 return iload_filename, iload_dirname, iload_glob, iload
2596class Inconsistency(Exception):
2597 pass
2600def consistency_check(list_of_tuples, message='values differ:'):
2601 '''
2602 Check for inconsistencies.
2604 Given a list of tuples, check that all tuple elements except for first one
2605 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be
2606 valid because the coordinates at the two channels are the same.
2607 '''
2609 if len(list_of_tuples) >= 2:
2610 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]):
2611 raise Inconsistency('%s\n' % message + '\n'.join(
2612 ' %s: %s' % (t[0], ', '.join('%g' % x for x in t[1:]))
2613 for t in list_of_tuples))
2616class defaultzerodict(dict):
2617 def __missing__(self, k):
2618 return 0
2621def mostfrequent(x):
2622 c = defaultzerodict()
2623 for e in x:
2624 c[e] += 1
2626 return sorted(list(c.keys()), key=lambda k: c[k])[-1]
2629def consistency_merge(list_of_tuples,
2630 message='values differ:',
2631 error='raise',
2632 merge=mostfrequent):
2634 assert error in ('raise', 'warn', 'ignore')
2636 if len(list_of_tuples) == 0:
2637 raise Exception('cannot merge empty sequence')
2639 try:
2640 consistency_check(list_of_tuples, message)
2641 return list_of_tuples[0][1:]
2642 except Inconsistency as e:
2643 if error == 'raise':
2644 raise
2646 elif error == 'warn':
2647 logger.warning(str(e))
2649 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]])
2652def short_to_list(nmax, it):
2653 import itertools
2655 if isinstance(it, list):
2656 return it
2658 li = []
2659 for i in range(nmax+1):
2660 try:
2661 li.append(next(it))
2662 except StopIteration:
2663 return li
2665 return itertools.chain(li, it)
2668def parse_md(f):
2669 try:
2670 with open(op.join(
2671 op.dirname(op.abspath(f)),
2672 'README.md'), 'r') as readme:
2673 mdstr = readme.read()
2674 except IOError as e:
2675 return 'Failed to get README.md: %s' % e
2677 # Remve the title
2678 mdstr = re.sub(r'^# .*\n?', '', mdstr)
2679 # Append sphinx reference to `pyrocko.` modules
2680 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr)
2681 # Convert Subsections to toc-less rubrics
2682 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr)
2683 return mdstr
2686def mpl_show(plt):
2687 import matplotlib
2688 if matplotlib.get_backend().lower() == 'agg':
2689 logger.warning('Cannot show() when using matplotlib "agg" backend')
2690 else:
2691 plt.show()
2694g_re_qsplit = re.compile(
2695 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)')
2696g_re_qsplit_sep = {}
2699def get_re_qsplit(sep):
2700 if sep is None:
2701 return g_re_qsplit
2702 else:
2703 if sep not in g_re_qsplit_sep:
2704 assert len(sep) == 1
2705 assert sep not in '\'"'
2706 esep = re.escape(sep)
2707 g_re_qsplit_sep[sep] = re.compile(
2708 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|'
2709 + r'([^' + esep + r']+|(?<=' + esep + r')(?=' + esep + r')|^(?=' + esep + r')|(?<=' + esep + r')$)') # noqa
2710 return g_re_qsplit_sep[sep]
2713g_re_trivial = re.compile(r'\A[^\'"\s]+\Z')
2714g_re_trivial_sep = {}
2717def get_re_trivial(sep):
2718 if sep is None:
2719 return g_re_trivial
2720 else:
2721 if sep not in g_re_qsplit_sep:
2722 assert len(sep) == 1
2723 assert sep not in '\'"'
2724 esep = re.escape(sep)
2725 g_re_trivial_sep[sep] = re.compile(r'\A[^\'"' + esep + r']+\Z')
2727 return g_re_trivial_sep[sep]
2730g_re_escape_s = re.compile(r'([\\\'])')
2731g_re_unescape_s = re.compile(r'\\([\\\'])')
2732g_re_escape_d = re.compile(r'([\\"])')
2733g_re_unescape_d = re.compile(r'\\([\\"])')
2736def escape_s(s):
2737 '''
2738 Backslash-escape single-quotes and backslashes.
2740 Example: ``Jack's`` => ``Jack\\'s``
2742 '''
2743 return g_re_escape_s.sub(r'\\\1', s)
2746def unescape_s(s):
2747 '''
2748 Unescape backslash-escaped single-quotes and backslashes.
2750 Example: ``Jack\\'s`` => ``Jack's``
2751 '''
2752 return g_re_unescape_s.sub(r'\1', s)
2755def escape_d(s):
2756 '''
2757 Backslash-escape double-quotes and backslashes.
2759 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"``
2760 '''
2761 return g_re_escape_d.sub(r'\\\1', s)
2764def unescape_d(s):
2765 '''
2766 Unescape backslash-escaped double-quotes and backslashes.
2768 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"``
2769 '''
2770 return g_re_unescape_d.sub(r'\1', s)
2773def qjoin_s(it, sep=None):
2774 '''
2775 Join sequence of strings into a line, single-quoting non-trivial strings.
2777 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"``
2778 '''
2779 re_trivial = get_re_trivial(sep)
2781 if sep is None:
2782 sep = ' '
2784 return sep.join(
2785 w if re_trivial.search(w) else "'%s'" % escape_s(w) for w in it)
2788def qjoin_d(it, sep=None):
2789 '''
2790 Join sequence of strings into a line, double-quoting non-trivial strings.
2792 Example: ``['55', 'Pete "The Robot" Smith']`` =>
2793 ``'55' "Pete \\\\"The Robot\\\\" Smith"'``
2794 '''
2795 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_d(w) for w in it)
2803def qsplit(s, sep=None):
2804 '''
2805 Split line into list of strings, allowing for quoted strings.
2807 Example: ``"55 'Sparrow\\\\'s Island'"`` =>
2808 ``["55", "Sparrow's Island"]``,
2809 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` =>
2810 ``['55', 'Pete "The Robot" Smith']``
2811 '''
2812 re_qsplit = get_re_qsplit(sep)
2813 return [
2814 (unescape_d(x[0]) or unescape_s(x[1]) or x[2])
2815 for x in re_qsplit.findall(s)]
2818g_have_warned_threadpoolctl = False
2821class threadpool_limits_dummy(object):
2823 def __init__(self, *args, **kwargs):
2824 pass
2826 def __enter__(self):
2827 global g_have_warned_threadpoolctl
2829 if not g_have_warned_threadpoolctl:
2830 logger.warning(
2831 'Cannot control number of BLAS threads because '
2832 '`threadpoolctl` module is not available. You may want to '
2833 'install `threadpoolctl`.')
2835 g_have_warned_threadpoolctl = True
2837 return self
2839 def __exit__(self, type, value, traceback):
2840 pass
2843def get_threadpool_limits():
2844 '''
2845 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail.
2846 '''
2848 try:
2849 from threadpoolctl import threadpool_limits
2850 return threadpool_limits
2852 except ImportError:
2853 return threadpool_limits_dummy