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
87try:
88 from urllib.parse import urlencode, quote, unquote # noqa
89 from urllib.request import (
90 Request, build_opener, HTTPDigestAuthHandler, urlopen as _urlopen) # noqa
91 from urllib.error import HTTPError, URLError # noqa
93except ImportError:
94 from urllib import urlencode, quote, unquote # noqa
95 from urllib2 import (Request, build_opener, HTTPDigestAuthHandler, # noqa
96 HTTPError, URLError, urlopen as _urlopen) # noqa
99class URLErrorSSL(URLError):
101 def __init__(self, urlerror):
102 self.__dict__ = urlerror.__dict__.copy()
104 def __str__(self):
105 return (
106 'Requesting web resource failed and the problem could be '
107 'related to SSL. Python standard libraries on some older '
108 'systems (like Ubuntu 14.04) are known to have trouble '
109 'with some SSL setups of today\'s servers: %s'
110 % URLError.__str__(self))
113def urlopen_ssl_check(*args, **kwargs):
114 try:
115 return urlopen(*args, **kwargs)
116 except URLError as e:
117 if str(e).find('SSL') != -1:
118 raise URLErrorSSL(e)
119 else:
120 raise
123def urlopen(*args, **kwargs):
124 if 'cafile' not in kwargs:
125 try:
126 import certifi
127 kwargs['cafile'] = certifi.where()
128 except ImportError:
129 pass
131 return _urlopen(*args, **kwargs)
134try:
135 long
136except NameError:
137 long = int
140force_dummy_progressbar = False
143try:
144 from pyrocko import util_ext
145except ImportError:
146 util_ext = None
149logger = logging.getLogger('pyrocko.util')
151try:
152 import progressbar as progressbar_mod
153except ImportError:
154 from pyrocko import dummy_progressbar as progressbar_mod
157try:
158 num_full = num.full
159except AttributeError:
160 def num_full(shape, fill_value, dtype=None, order='C'):
161 a = num.empty(shape, dtype=dtype, order=order)
162 a.fill(fill_value)
163 return a
165try:
166 num_full_like = num.full_like
167except AttributeError:
168 def num_full_like(arr, fill_value, dtype=None, order='K', subok=True):
169 a = num.empty_like(arr, dtype=dtype, order=order, subok=subok)
170 a.fill(fill_value)
171 return a
174def progressbar_module():
175 return progressbar_mod
178g_setup_logging_args = 'pyrocko', 'warning'
181def setup_logging(programname='pyrocko', levelname='warning'):
182 '''
183 Initialize logging.
185 :param programname: program name to be written in log
186 :param levelname: string indicating the logging level ('debug', 'info',
187 'warning', 'error', 'critical')
189 This function is called at startup by most pyrocko programs to set up a
190 consistent logging format. This is simply a shortcut to a call to
191 :py:func:`logging.basicConfig()`.
192 '''
194 global g_setup_logging_args
195 g_setup_logging_args = (programname, levelname)
197 levels = {'debug': logging.DEBUG,
198 'info': logging.INFO,
199 'warning': logging.WARNING,
200 'error': logging.ERROR,
201 'critical': logging.CRITICAL}
203 logging.basicConfig(
204 level=levels[levelname],
205 format=programname+':%(name)-20s - %(levelname)-8s - %(message)s')
208def subprocess_setup_logging_args():
209 '''
210 Get arguments from previous call to setup_logging.
212 These can be sent down to a worker process so it can setup its logging
213 in the same way as the main process.
214 '''
215 return g_setup_logging_args
218def data_file(fn):
219 return os.path.join(os.path.split(__file__)[0], 'data', fn)
222class DownloadError(Exception):
223 pass
226class PathExists(DownloadError):
227 pass
230def _download(url, fpath, username=None, password=None,
231 force=False, method='download', stats=None,
232 status_callback=None, entries_wanted=None,
233 recursive=False, header=None):
235 import requests
236 from requests.auth import HTTPBasicAuth
237 from requests.exceptions import HTTPError as req_HTTPError
239 requests.adapters.DEFAULT_RETRIES = 5
240 urljoin = requests.compat.urljoin
242 session = requests.Session()
243 session.header = header
244 session.auth = None if username is None\
245 else HTTPBasicAuth(username, password)
247 status = {
248 'ntotal_files': 0,
249 'nread_files': 0,
250 'ntotal_bytes_all_files': 0,
251 'nread_bytes_all_files': 0,
252 'ntotal_bytes_current_file': 0,
253 'nread_bytes_current_file': 0,
254 'finished': False
255 }
257 try:
258 url_to_size = {}
260 if callable(status_callback):
261 status_callback(status)
263 if not recursive and url.endswith('/'):
264 raise DownloadError(
265 'URL: %s appears to be a directory'
266 ' but recurvise download is False' % url)
268 if recursive and not url.endswith('/'):
269 url += '/'
271 def parse_directory_tree(url, subdir=''):
272 r = session.get(urljoin(url, subdir))
273 r.raise_for_status()
275 entries = re.findall(r'href="([a-zA-Z0-9_.-]+/?)"', r.text)
277 files = sorted(set(subdir + fn for fn in entries
278 if not fn.endswith('/')))
280 if entries_wanted is not None:
281 files = [fn for fn in files
282 if (fn in wanted for wanted in entries_wanted)]
284 status['ntotal_files'] += len(files)
286 dirs = sorted(set(subdir + dn for dn in entries
287 if dn.endswith('/')
288 and dn not in ('./', '../')))
290 for dn in dirs:
291 files.extend(parse_directory_tree(
292 url, subdir=dn))
294 return files
296 def get_content_length(url):
297 if url not in url_to_size:
298 r = session.head(url, headers={'Accept-Encoding': ''})
300 content_length = r.headers.get('content-length', None)
301 if content_length is None:
302 logger.debug('Could not get HTTP header '
303 'Content-Length for %s' % url)
305 content_length = None
307 else:
308 content_length = int(content_length)
309 status['ntotal_bytes_all_files'] += content_length
310 if callable(status_callback):
311 status_callback(status)
313 url_to_size[url] = content_length
315 return url_to_size[url]
317 def download_file(url, fn):
318 logger.info('starting download of %s...' % url)
319 ensuredirs(fn)
321 fsize = get_content_length(url)
322 status['ntotal_bytes_current_file'] = fsize
323 status['nread_bytes_current_file'] = 0
324 if callable(status_callback):
325 status_callback(status)
327 r = session.get(url, stream=True, timeout=5)
328 r.raise_for_status()
330 frx = 0
331 fn_tmp = fn + '.%i.temp' % os.getpid()
332 with open(fn_tmp, 'wb') as f:
333 for d in r.iter_content(chunk_size=1024):
334 f.write(d)
335 frx += len(d)
337 status['nread_bytes_all_files'] += len(d)
338 status['nread_bytes_current_file'] += len(d)
339 if callable(status_callback):
340 status_callback(status)
342 os.rename(fn_tmp, fn)
344 if fsize is not None and frx != fsize:
345 logger.warning(
346 'HTTP header Content-Length: %i bytes does not match '
347 'download size %i bytes' % (fsize, frx))
349 logger.info('finished download of %s' % url)
351 status['nread_files'] += 1
352 if callable(status_callback):
353 status_callback(status)
355 if recursive:
356 if op.exists(fpath) and not force:
357 raise PathExists('path %s already exists' % fpath)
359 files = parse_directory_tree(url)
361 dsize = 0
362 for fn in files:
363 file_url = urljoin(url, fn)
364 dsize += get_content_length(file_url) or 0
366 if method == 'calcsize':
367 return dsize
369 else:
370 for fn in files:
371 file_url = urljoin(url, fn)
372 download_file(file_url, op.join(fpath, fn))
374 else:
375 status['ntotal_files'] += 1
376 if callable(status_callback):
377 status_callback(status)
379 fsize = get_content_length(url)
380 if method == 'calcsize':
381 return fsize
382 else:
383 download_file(url, fpath)
385 except req_HTTPError as e:
386 logging.warn("http error: %s" % e)
387 raise DownloadError('could not download file(s) from: %s' % url)
389 finally:
390 status['finished'] = True
391 if callable(status_callback):
392 status_callback(status)
393 session.close()
396def download_file(
397 url, fpath, username=None, password=None, status_callback=None,
398 **kwargs):
399 return _download(
400 url, fpath, username, password,
401 recursive=False,
402 status_callback=status_callback,
403 **kwargs)
406def download_dir(
407 url, fpath, username=None, password=None, status_callback=None,
408 **kwargs):
410 return _download(
411 url, fpath, username, password,
412 recursive=True,
413 status_callback=status_callback,
414 **kwargs)
417class HPFloatUnavailable(Exception):
418 pass
421class dummy_hpfloat(object):
422 def __init__(self, *args, **kwargs):
423 raise HPFloatUnavailable(
424 'NumPy lacks support for float128 or float96 data type on this '
425 'platform.')
428if hasattr(num, 'float128'):
429 hpfloat = num.float128
430elif hasattr(num, 'float96'):
431 hpfloat = num.float96
432else:
433 hpfloat = dummy_hpfloat
436g_time_float = None
437g_time_dtype = None
440class TimeFloatSettingError(Exception):
441 pass
444def use_high_precision_time(enabled):
445 '''
446 Globally force a specific time handling mode.
448 See :ref:`High precision time handling mode <time-handling-mode>`.
450 :param enabled: enable/disable use of high precision time type
451 :type enabled: bool
453 This function should be called before handling/reading any time data.
454 It can only be called once.
456 Special attention is required when using multiprocessing on a platform
457 which does not use fork under the hood. In such cases, the desired setting
458 must be set also in the subprocess.
459 '''
460 _setup_high_precision_time_mode(enabled_app=enabled)
463def _setup_high_precision_time_mode(enabled_app=False):
464 global g_time_float
465 global g_time_dtype
467 if not (g_time_float is None and g_time_dtype is None):
468 raise TimeFloatSettingError(
469 'Cannot set time handling mode: too late, it has already been '
470 'fixed by an earlier call.')
472 from pyrocko import config
474 conf = config.config()
475 enabled_config = conf.use_high_precision_time
477 enabled_env = os.environ.get('PYROCKO_USE_HIGH_PRECISION_TIME', None)
478 if enabled_env is not None:
479 try:
480 enabled_env = int(enabled_env) == 1
481 except ValueError:
482 raise TimeFloatSettingError(
483 'Environment variable PYROCKO_USE_HIGH_PRECISION_TIME '
484 'should be set to 0 or 1.')
486 enabled = enabled_config
487 mode_from = 'config variable `use_high_precision_time`'
488 notify = enabled
490 if enabled_env is not None:
491 if enabled_env != enabled:
492 notify = True
493 enabled = enabled_env
494 mode_from = 'environment variable `PYROCKO_USE_HIGH_PRECISION_TIME`'
496 if enabled_app is not None:
497 if enabled_app != enabled:
498 notify = True
499 enabled = enabled_app
500 mode_from = 'application override'
502 logger.debug('''
503Pyrocko high precision time mode selection (latter override earlier):
504 config: %s
505 env: %s
506 app: %s
507 -> enabled: %s'''.lstrip() % (
508 enabled_config, enabled_env, enabled_app, enabled))
510 if notify:
511 logger.info('Pyrocko high precision time mode %s by %s.' % (
512 'activated' if enabled else 'deactivated',
513 mode_from))
515 if enabled:
516 g_time_float = hpfloat
517 g_time_dtype = hpfloat
518 else:
519 g_time_float = float
520 g_time_dtype = num.float64
523def get_time_float():
524 '''
525 Get the effective float class for timestamps.
527 See :ref:`High precision time handling mode <time-handling-mode>`.
529 :returns: :py:class:`float` or :py:class:`hpfloat`, depending on the
530 current time handling mode
531 '''
532 global g_time_float
534 if g_time_float is None:
535 _setup_high_precision_time_mode()
537 return g_time_float
540def get_time_dtype():
541 '''
542 Get effective NumPy float class to handle timestamps.
544 See :ref:`High precision time handling mode <time-handling-mode>`.
545 '''
547 global g_time_dtype
549 if g_time_dtype is None:
550 _setup_high_precision_time_mode()
552 return g_time_dtype
555def to_time_float(t):
556 '''
557 Convert float to valid time stamp in the current time handling mode.
559 See :ref:`High precision time handling mode <time-handling-mode>`.
560 '''
561 return get_time_float()(t)
564class TimestampTypeError(ValueError):
565 pass
568def check_time_class(t, error='raise'):
569 '''
570 Type-check variable against current time handling mode.
572 See :ref:`High precision time handling mode <time-handling-mode>`.
573 '''
575 if t == 0.0:
576 return
578 if not isinstance(t, get_time_float()):
579 message = (
580 'Timestamp %g is of type %s but should be of type %s with '
581 'Pyrocko\'s currently selected time handling mode.\n\n'
582 'See https://pyrocko.org/docs/current/library/reference/util.html'
583 '#high-precision-time-handling-mode' % (
584 t, type(t), get_time_float()))
586 if error == 'raise':
587 raise TimestampTypeError(message)
588 elif error == 'warn':
589 logger.warn(message)
590 else:
591 assert False
594class Stopwatch(object):
595 '''
596 Simple stopwatch to measure elapsed wall clock time.
598 Usage::
600 s = Stopwatch()
601 time.sleep(1)
602 print s()
603 time.sleep(1)
604 print s()
605 '''
607 def __init__(self):
608 self.start = time.time()
610 def __call__(self):
611 return time.time() - self.start
614def wrap(text, line_length=80):
615 '''
616 Paragraph and list-aware wrapping of text.
617 '''
619 text = text.strip('\n')
620 at_lineend = re.compile(r' *\n')
621 at_para = re.compile(r'((^|(\n\s*)?\n)(\s+[*] )|\n\s*\n)')
623 paragraphs = at_para.split(text)[::5]
624 listindents = at_para.split(text)[4::5]
625 newlist = at_para.split(text)[3::5]
627 listindents[0:0] = [None]
628 listindents.append(True)
629 newlist.append(None)
631 det_indent = re.compile(r'^ *')
633 outlines = []
634 for ip, p in enumerate(paragraphs):
635 if not p:
636 continue
638 if listindents[ip] is None:
639 _indent = det_indent.findall(p)[0]
640 findent = _indent
641 else:
642 findent = listindents[ip]
643 _indent = ' ' * len(findent)
645 ll = line_length - len(_indent)
646 llf = ll
648 oldlines = [s.strip() for s in at_lineend.split(p.rstrip())]
649 p1 = ' '.join(oldlines)
650 possible = re.compile(r'(^.{1,%i}|.{1,%i})( |$)' % (llf, ll))
651 for imatch, match in enumerate(possible.finditer(p1)):
652 parout = match.group(1)
653 if imatch == 0:
654 outlines.append(findent + parout)
655 else:
656 outlines.append(_indent + parout)
658 if ip != len(paragraphs)-1 and (
659 listindents[ip] is None
660 or newlist[ip] is not None
661 or listindents[ip+1] is None):
663 outlines.append('')
665 return outlines
668class BetterHelpFormatter(optparse.IndentedHelpFormatter):
670 def __init__(self, *args, **kwargs):
671 optparse.IndentedHelpFormatter.__init__(self, *args, **kwargs)
673 def format_option(self, option):
674 '''
675 From IndentedHelpFormatter but using a different wrap method.
676 '''
678 help_text_position = 4 + self.current_indent
679 help_text_width = self.width - help_text_position
681 opts = self.option_strings[option]
682 opts = "%*s%s" % (self.current_indent, "", opts)
683 if option.help:
684 help_text = self.expand_default(option)
686 if self.help_position + len(help_text) + 1 <= self.width:
687 lines = [
688 '',
689 '%-*s %s' % (self.help_position, opts, help_text),
690 '']
691 else:
692 lines = ['']
693 lines.append(opts)
694 lines.append('')
695 if option.help:
696 help_lines = wrap(help_text, help_text_width)
697 lines.extend(["%*s%s" % (help_text_position, "", line)
698 for line in help_lines])
699 lines.append('')
701 return "\n".join(lines)
703 def format_description(self, description):
704 if not description:
705 return ''
707 if self.current_indent == 0:
708 lines = []
709 else:
710 lines = ['']
712 lines.extend(wrap(description, self.width - self.current_indent))
713 if self.current_indent == 0:
714 lines.append('\n')
716 return '\n'.join(
717 ['%*s%s' % (self.current_indent, '', line) for line in lines])
720def progressbar(label, maxval):
721 progressbar_mod = progressbar_module()
722 if force_dummy_progressbar:
723 progressbar_mod = dummy_progressbar
725 widgets = [
726 label, ' ',
727 progressbar_mod.Bar(marker='-', left='[', right=']'), ' ',
728 progressbar_mod.Percentage(), ' ',
729 progressbar_mod.ETA()]
731 pbar = progressbar_mod.ProgressBar(widgets=widgets, maxval=maxval).start()
732 return pbar
735def progress_beg(label):
736 '''
737 Notify user that an operation has started.
739 :param label: name of the operation
741 To be used in conjuction with :py:func:`progress_end`.
742 '''
744 sys.stderr.write(label)
745 sys.stderr.flush()
748def progress_end(label=''):
749 '''
750 Notify user that an operation has ended.
752 :param label: name of the operation
754 To be used in conjuction with :py:func:`progress_beg`.
755 '''
757 sys.stderr.write(' done. %s\n' % label)
758 sys.stderr.flush()
761class ArangeError(Exception):
762 pass
765def arange2(start, stop, step, dtype=float, epsilon=1e-6, error='raise'):
766 '''
767 Return evenly spaced numbers over a specified interval.
769 Like :py:func:`numpy.arange` but returning floating point numbers by
770 default and with defined behaviour when stepsize is inconsistent with
771 interval bounds. It is considered inconsistent if the difference between
772 the closest multiple of ``step`` and ``stop`` is larger than ``epsilon *
773 step``. Inconsistencies are handled according to the ``error`` parameter.
774 If it is set to ``'raise'`` an exception of type :py:exc:`ArangeError` is
775 raised. If it is set to ``'round'``, ``'floor'``, or ``'ceil'``, ``stop``
776 is silently changed to the closest, the next smaller, or next larger
777 multiple of ``step``, respectively.
778 '''
780 assert error in ('raise', 'round', 'floor', 'ceil')
782 start = dtype(start)
783 stop = dtype(stop)
784 step = dtype(step)
786 rnd = {'floor': math.floor, 'ceil': math.ceil}.get(error, round)
788 n = int(rnd((stop - start) / step)) + 1
789 stop_check = start + (n-1) * step
791 if error == 'raise' and abs(stop_check - stop) > step * epsilon:
792 raise ArangeError(
793 'inconsistent range specification: start=%g, stop=%g, step=%g'
794 % (start, stop, step))
796 x = num.arange(n, dtype=dtype)
797 x *= step
798 x += start
799 return x
802def polylinefit(x, y, n_or_xnodes):
803 '''
804 Fit piece-wise linear function to data.
806 :param x,y: arrays with coordinates of data
807 :param n_or_xnodes: int, number of segments or x coordinates of polyline
809 :returns: `(xnodes, ynodes, rms_error)` arrays with coordinates of
810 polyline, root-mean-square error
811 '''
813 x = num.asarray(x)
814 y = num.asarray(y)
816 if isinstance(n_or_xnodes, int):
817 n = n_or_xnodes
818 xmin = x.min()
819 xmax = x.max()
820 xnodes = num.linspace(xmin, xmax, n+1)
821 else:
822 xnodes = num.asarray(n_or_xnodes)
823 n = xnodes.size - 1
825 assert len(x) == len(y)
826 assert n > 0
828 ndata = len(x)
829 a = num.zeros((ndata+(n-1), n*2))
830 for i in range(n):
831 xmin_block = xnodes[i]
832 xmax_block = xnodes[i+1]
833 if i == n-1: # don't loose last point
834 indices = num.where(
835 num.logical_and(xmin_block <= x, x <= xmax_block))[0]
836 else:
837 indices = num.where(
838 num.logical_and(xmin_block <= x, x < xmax_block))[0]
840 a[indices, i*2] = x[indices]
841 a[indices, i*2+1] = 1.0
843 w = float(ndata)*100.
844 if i < n-1:
845 a[ndata+i, i*2] = xmax_block*w
846 a[ndata+i, i*2+1] = 1.0*w
847 a[ndata+i, i*2+2] = -xmax_block*w
848 a[ndata+i, i*2+3] = -1.0*w
850 d = num.concatenate((y, num.zeros(n-1)))
851 model = num.linalg.lstsq(a, d, rcond=-1)[0].reshape((n, 2))
853 ynodes = num.zeros(n+1)
854 ynodes[:n] = model[:, 0]*xnodes[:n] + model[:, 1]
855 ynodes[1:] += model[:, 0]*xnodes[1:] + model[:, 1]
856 ynodes[1:n] *= 0.5
858 rms_error = num.sqrt(num.mean((num.interp(x, xnodes, ynodes) - y)**2))
860 return xnodes, ynodes, rms_error
863def plf_integrate_piecewise(x_edges, x, y):
864 '''
865 Calculate definite integral of piece-wise linear function on intervals.
867 Use trapezoidal rule to calculate definite integral of a piece-wise linear
868 function for a series of consecutive intervals. ``x_edges`` and ``x`` must
869 be sorted.
871 :param x_edges: array with edges of the intervals
872 :param x,y: arrays with coordinates of piece-wise linear function's
873 control points
874 '''
876 x_all = num.concatenate((x, x_edges))
877 ii = num.argsort(x_all)
878 y_edges = num.interp(x_edges, x, y)
879 y_all = num.concatenate((y, y_edges))
880 xs = x_all[ii]
881 ys = y_all[ii]
882 y_all[ii[1:]] = num.cumsum((xs[1:] - xs[:-1]) * 0.5 * (ys[1:] + ys[:-1]))
883 return num.diff(y_all[-len(y_edges):])
886def diff_fd_1d_4o(dt, data):
887 '''
888 Approximate first derivative of an array (forth order, central FD).
890 :param dt: sampling interval
891 :param data: NumPy array with data samples
893 :returns: NumPy array with same shape as input
895 Interior points are approximated to fourth order, edge points to first
896 order right- or left-sided respectively, points next to edge to second
897 order central.
898 '''
899 import scipy.signal
901 ddata = num.empty_like(data, dtype=float)
903 if data.size >= 5:
904 ddata[2:-2] = scipy.signal.lfilter(
905 [-1., +8., 0., -8., 1.], [1.], data)[4:] / (12.*dt)
907 if data.size >= 3:
908 ddata[1] = (data[2] - data[0]) / (2. * dt)
909 ddata[-2] = (data[-1] - data[-3]) / (2. * dt)
911 if data.size >= 2:
912 ddata[0] = (data[1] - data[0]) / dt
913 ddata[-1] = (data[-1] - data[-2]) / dt
915 if data.size == 1:
916 ddata[0] = 0.0
918 return ddata
921def diff_fd_1d_2o(dt, data):
922 '''
923 Approximate first derivative of an array (second order, central FD).
925 :param dt: sampling interval
926 :param data: NumPy array with data samples
928 :returns: NumPy array with same shape as input
930 Interior points are approximated to second order, edge points to first
931 order right- or left-sided respectively.
933 Uses :py:func:`numpy.gradient`.
934 '''
936 return num.gradient(data, dt)
939def diff_fd_2d_4o(dt, data):
940 '''
941 Approximate second derivative of an array (forth 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 fourth order, next-to-edge points to
949 second order, edge points repeated.
950 '''
951 import scipy.signal
953 ddata = num.empty_like(data, dtype=float)
955 if data.size >= 5:
956 ddata[2:-2] = scipy.signal.lfilter(
957 [-1., +16., -30., +16., -1.], [1.], data)[4:] / (12.*dt**2)
959 if data.size >= 3:
960 ddata[:2] = (data[2] - 2.0 * data[1] + data[0]) / dt**2
961 ddata[-2:] = (data[-1] - 2.0 * data[-2] + data[-3]) / dt**2
963 if data.size < 3:
964 ddata[:] = 0.0
966 return ddata
969def diff_fd_2d_2o(dt, data):
970 '''
971 Approximate second derivative of an array (second order, central FD).
973 :param dt: sampling interval
974 :param data: NumPy array with data samples
976 :returns: NumPy array with same shape as input
978 Interior points are approximated to second order, edge points repeated.
979 '''
980 import scipy.signal
982 ddata = num.empty_like(data, dtype=float)
984 if data.size >= 3:
985 ddata[1:-1] = scipy.signal.lfilter(
986 [1., -2., 1.], [1.], data)[2:] / (dt**2)
988 ddata[0] = ddata[1]
989 ddata[-1] = ddata[-2]
991 if data.size < 3:
992 ddata[:] = 0.0
994 return ddata
997def diff_fd(n, order, dt, data):
998 '''
999 Approximate 1st or 2nd derivative of an array.
1001 :param n: 1 for first derivative, 2 for second
1002 :param order: order of the approximation 2 and 4 are supported
1003 :param dt: sampling interval
1004 :param data: NumPy array with data samples
1006 :returns: NumPy array with same shape as input
1008 This is a frontend to the functions :py:func:`diff_fd_1d_2o`,
1009 :py:func:`diff_fd_1d_4o`, :py:func:`diff_fd_2d_2o`, and
1010 :py:func:`diff_fd_2d_4o`.
1012 Raises :py:exc:`ValueError` for unsupported `n` or `order`.
1013 '''
1015 funcs = {
1016 1: {2: diff_fd_1d_2o, 4: diff_fd_1d_4o},
1017 2: {2: diff_fd_2d_2o, 4: diff_fd_2d_4o}}
1019 try:
1020 funcs_n = funcs[n]
1021 except KeyError:
1022 raise ValueError(
1023 'pyrocko.util.diff_fd: '
1024 'Only 1st and 2sd derivatives are supported.')
1026 try:
1027 func = funcs_n[order]
1028 except KeyError:
1029 raise ValueError(
1030 'pyrocko.util.diff_fd: '
1031 'Order %i is not supported for %s derivative. Supported: %s' % (
1032 order, ['', '1st', '2nd'][n],
1033 ', '.join('%i' % order for order in sorted(funcs_n.keys()))))
1035 return func(dt, data)
1038class GlobalVars(object):
1039 reuse_store = dict()
1040 decitab_nmax = 0
1041 decitab = {}
1042 decimate_fir_coeffs = {}
1043 decimate_iir_coeffs = {}
1044 re_frac = None
1047def decimate_coeffs(q, n=None, ftype='iir'):
1049 q = int(q)
1051 if n is None:
1052 if ftype == 'fir':
1053 n = 30
1054 else:
1055 n = 8
1057 if ftype == 'fir':
1058 coeffs = GlobalVars.decimate_fir_coeffs
1059 if (n, 1./q) not in coeffs:
1060 coeffs[n, 1./q] = signal.firwin(n+1, 1./q, window='hamming')
1062 b = coeffs[n, 1./q]
1063 return b, [1.], n
1065 else:
1066 coeffs = GlobalVars.decimate_iir_coeffs
1067 if (n, 0.05, 0.8/q) not in coeffs:
1068 coeffs[n, 0.05, 0.8/q] = signal.cheby1(n, 0.05, 0.8/q)
1070 b, a = coeffs[n, 0.05, 0.8/q]
1071 return b, a, n
1074def decimate(x, q, n=None, ftype='iir', zi=None, ioff=0):
1075 '''
1076 Downsample the signal x by an integer factor q, using an order n filter
1078 By default, an order 8 Chebyshev type I filter is used or a 30 point FIR
1079 filter with hamming window if ftype is 'fir'.
1081 :param x: the signal to be downsampled (1D NumPy array)
1082 :param q: the downsampling factor
1083 :param n: order of the filter (1 less than the length of the filter for a
1084 'fir' filter)
1085 :param ftype: type of the filter; can be 'iir' or 'fir'
1087 :returns: the downsampled signal (1D NumPy array)
1089 '''
1091 b, a, n = decimate_coeffs(q, n, ftype)
1093 if zi is None or zi is True:
1094 zi_ = num.zeros(max(len(a), len(b))-1, dtype=float)
1095 else:
1096 zi_ = zi
1098 y, zf = signal.lfilter(b, a, x, zi=zi_)
1100 if zi is not None:
1101 return y[n//2+ioff::q].copy(), zf
1102 else:
1103 return y[n//2+ioff::q].copy()
1106class UnavailableDecimation(Exception):
1107 '''
1108 Exception raised by :py:func:`decitab` for unavailable decimation factors.
1109 '''
1111 pass
1114def gcd(a, b, epsilon=1e-7):
1115 '''
1116 Greatest common divisor.
1117 '''
1119 while b > epsilon*a:
1120 a, b = b, a % b
1122 return a
1125def lcm(a, b):
1126 '''
1127 Least common multiple.
1128 '''
1130 return a*b // gcd(a, b)
1133def mk_decitab(nmax=100):
1134 '''
1135 Make table with decimation sequences.
1137 Decimation from one sampling rate to a lower one is achieved by a
1138 successive application of :py:func:`decimation` with small integer
1139 downsampling factors (because using large downampling factors can make the
1140 decimation unstable or slow). This function sets up a table with downsample
1141 sequences for factors up to ``nmax``.
1142 '''
1144 tab = GlobalVars.decitab
1145 for i in range(1, 10):
1146 for j in range(1, i+1):
1147 for k in range(1, j+1):
1148 for l_ in range(1, k+1):
1149 for m in range(1, l_+1):
1150 p = i*j*k*l_*m
1151 if p > nmax:
1152 break
1153 if p not in tab:
1154 tab[p] = (i, j, k, l_, m)
1155 if i*j*k*l_ > nmax:
1156 break
1157 if i*j*k > nmax:
1158 break
1159 if i*j > nmax:
1160 break
1161 if i > nmax:
1162 break
1164 GlobalVars.decitab_nmax = nmax
1167def zfmt(n):
1168 return '%%0%ii' % (int(math.log10(n - 1)) + 1)
1171def _year_to_time(year):
1172 tt = (year, 1, 1, 0, 0, 0)
1173 return to_time_float(calendar.timegm(tt))
1176def _working_year(year):
1177 try:
1178 tt = (year, 1, 1, 0, 0, 0)
1179 t = calendar.timegm(tt)
1180 tt2_ = time.gmtime(t)
1181 tt2 = tuple(tt2_)[:6]
1182 if tt != tt2:
1183 return False
1185 s = '%i-01-01 00:00:00' % year
1186 s2 = time.strftime('%Y-%m-%d %H:%M:%S', tt2_)
1187 if s != s2:
1188 return False
1190 t3 = str_to_time(s2, format='%Y-%m-%d %H:%M:%S')
1191 s3 = time_to_str(t3, format='%Y-%m-%d %H:%M:%S')
1192 if s3 != s2:
1193 return False
1195 except Exception:
1196 return False
1198 return True
1201def working_system_time_range(year_min_lim=None, year_max_lim=None):
1202 '''
1203 Check time range supported by the systems's time conversion functions.
1205 Returns system time stamps of start of year of first/last fully supported
1206 year span. If this is before 1900 or after 2100, return first/last century
1207 which is fully supported.
1209 :returns: ``(tmin, tmax, year_min, year_max)``
1210 '''
1212 year0 = 2000
1213 year_min = year0
1214 year_max = year0
1216 itests = list(range(101))
1217 for i in range(19):
1218 itests.append(200 + i*100)
1220 for i in itests:
1221 year = year0 - i
1222 if year_min_lim is not None and year < year_min_lim:
1223 year_min = year_min_lim
1224 break
1225 elif not _working_year(year):
1226 break
1227 else:
1228 year_min = year
1230 for i in itests:
1231 year = year0 + i
1232 if year_max_lim is not None and year > year_max_lim:
1233 year_max = year_max_lim
1234 break
1235 elif not _working_year(year + 1):
1236 break
1237 else:
1238 year_max = year
1240 return (
1241 _year_to_time(year_min),
1242 _year_to_time(year_max),
1243 year_min, year_max)
1246g_working_system_time_range = None
1249def get_working_system_time_range():
1250 '''
1251 Caching variant of :py:func:`working_system_time_range`.
1252 '''
1254 global g_working_system_time_range
1255 if g_working_system_time_range is None:
1256 g_working_system_time_range = working_system_time_range()
1258 return g_working_system_time_range
1261def is_working_time(t):
1262 tmin, tmax, _, _ = get_working_system_time_range()
1263 return tmin <= t <= tmax
1266def julian_day_of_year(timestamp):
1267 '''
1268 Get the day number after the 1st of January of year in ``timestamp``.
1270 :returns: day number as int
1271 '''
1273 return time.gmtime(int(timestamp)).tm_yday
1276def hour_start(timestamp):
1277 '''
1278 Get beginning of hour for any point in time.
1280 :param timestamp: time instant as system timestamp (in seconds)
1282 :returns: instant of hour start as system timestamp
1283 '''
1285 tt = time.gmtime(int(timestamp))
1286 tts = tt[0:4] + (0, 0)
1287 return to_time_float(calendar.timegm(tts))
1290def day_start(timestamp):
1291 '''
1292 Get beginning of day for any point in time.
1294 :param timestamp: time instant as system timestamp (in seconds)
1296 :returns: instant of day start as system timestamp
1297 '''
1299 tt = time.gmtime(int(timestamp))
1300 tts = tt[0:3] + (0, 0, 0)
1301 return to_time_float(calendar.timegm(tts))
1304def month_start(timestamp):
1305 '''
1306 Get beginning of month for any point in time.
1308 :param timestamp: time instant as system timestamp (in seconds)
1310 :returns: instant of month start as system timestamp
1311 '''
1313 tt = time.gmtime(int(timestamp))
1314 tts = tt[0:2] + (1, 0, 0, 0)
1315 return to_time_float(calendar.timegm(tts))
1318def year_start(timestamp):
1319 '''
1320 Get beginning of year for any point in time.
1322 :param timestamp: time instant as system timestamp (in seconds)
1324 :returns: instant of year start as system timestamp
1325 '''
1327 tt = time.gmtime(int(timestamp))
1328 tts = tt[0:1] + (1, 1, 0, 0, 0)
1329 return to_time_float(calendar.timegm(tts))
1332def iter_days(tmin, tmax):
1333 '''
1334 Yields begin and end of days until given time span is covered.
1336 :param tmin,tmax: input time span
1338 :yields: tuples with (begin, end) of days as system timestamps
1339 '''
1341 t = day_start(tmin)
1342 while t < tmax:
1343 tend = day_start(t + 26*60*60)
1344 yield t, tend
1345 t = tend
1348def iter_months(tmin, tmax):
1349 '''
1350 Yields begin and end of months until given time span is covered.
1352 :param tmin,tmax: input time span
1354 :yields: tuples with (begin, end) of months as system timestamps
1355 '''
1357 t = month_start(tmin)
1358 while t < tmax:
1359 tend = month_start(t + 24*60*60*33)
1360 yield t, tend
1361 t = tend
1364def iter_years(tmin, tmax):
1365 '''
1366 Yields begin and end of years until given time span is covered.
1368 :param tmin,tmax: input time span
1370 :yields: tuples with (begin, end) of years as system timestamps
1371 '''
1373 t = year_start(tmin)
1374 while t < tmax:
1375 tend = year_start(t + 24*60*60*369)
1376 yield t, tend
1377 t = tend
1380def decitab(n):
1381 '''
1382 Get integer decimation sequence for given downampling factor.
1384 :param n: target decimation factor
1386 :returns: tuple with downsampling sequence
1387 '''
1389 if n > GlobalVars.decitab_nmax:
1390 mk_decitab(n*2)
1391 if n not in GlobalVars.decitab:
1392 raise UnavailableDecimation('ratio = %g' % n)
1393 return GlobalVars.decitab[n]
1396def ctimegm(s, format="%Y-%m-%d %H:%M:%S"):
1397 '''
1398 Convert string representing UTC time to system time.
1400 :param s: string to be interpreted
1401 :param format: format string passed to :py:func:`strptime`
1403 :returns: system time stamp
1405 Interpretes string with format ``'%Y-%m-%d %H:%M:%S'``, using strptime.
1407 .. note::
1408 This function is to be replaced by :py:func:`str_to_time`.
1409 '''
1411 return calendar.timegm(time.strptime(s, format))
1414def gmctime(t, format="%Y-%m-%d %H:%M:%S"):
1415 '''
1416 Get string representation from system time, UTC.
1418 Produces string with format ``'%Y-%m-%d %H:%M:%S'``, using strftime.
1420 .. note::
1421 This function is to be repaced by :py:func:`time_to_str`.
1422 '''
1424 return time.strftime(format, time.gmtime(t))
1427def gmctime_v(t, format="%a, %d %b %Y %H:%M:%S"):
1428 '''
1429 Get string representation from system time, UTC. Same as
1430 :py:func:`gmctime` but with a more verbose default format.
1432 .. note::
1433 This function is to be replaced by :py:func:`time_to_str`.
1434 '''
1436 return time.strftime(format, time.gmtime(t))
1439def gmctime_fn(t, format="%Y-%m-%d_%H-%M-%S"):
1440 '''
1441 Get string representation from system time, UTC. Same as
1442 :py:func:`gmctime` but with a default usable in filenames.
1444 .. note::
1445 This function is to be replaced by :py:func:`time_to_str`.
1446 '''
1448 return time.strftime(format, time.gmtime(t))
1451class TimeStrError(Exception):
1452 pass
1455class FractionalSecondsMissing(TimeStrError):
1456 '''
1457 Exception raised by :py:func:`str_to_time` when the given string lacks
1458 fractional seconds.
1459 '''
1461 pass
1464class FractionalSecondsWrongNumberOfDigits(TimeStrError):
1465 '''
1466 Exception raised by :py:func:`str_to_time` when the given string has an
1467 incorrect number of digits in the fractional seconds part.
1468 '''
1470 pass
1473def _endswith_n(s, endings):
1474 for ix, x in enumerate(endings):
1475 if s.endswith(x):
1476 return ix
1477 return -1
1480def str_to_time(s, format='%Y-%m-%d %H:%M:%S.OPTFRAC'):
1481 '''
1482 Convert string representing UTC time to floating point system time.
1484 :param s: string representing UTC time
1485 :param format: time string format
1486 :returns: system time stamp as floating point value
1488 Uses the semantics of :py:func:`time.strptime` but allows for fractional
1489 seconds. If the format ends with ``'.FRAC'``, anything after a dot is
1490 interpreted as fractional seconds. If the format ends with ``'.OPTFRAC'``,
1491 the fractional part, including the dot is made optional. The latter has the
1492 consequence, that the time strings and the format may not contain any other
1493 dots. If the format ends with ``'.xFRAC'`` where x is 1, 2, or 3, it is
1494 ensured, that exactly that number of digits are present in the fractional
1495 seconds.
1496 '''
1498 time_float = get_time_float()
1500 if util_ext is not None:
1501 try:
1502 t, tfrac = util_ext.stt(s, format)
1503 except util_ext.UtilExtError as e:
1504 raise TimeStrError(
1505 '%s, string=%s, format=%s' % (str(e), s, format))
1507 return time_float(t)+tfrac
1509 fracsec = 0.
1510 fixed_endings = '.FRAC', '.1FRAC', '.2FRAC', '.3FRAC'
1512 iend = _endswith_n(format, fixed_endings)
1513 if iend != -1:
1514 dotpos = s.rfind('.')
1515 if dotpos == -1:
1516 raise FractionalSecondsMissing(
1517 'string=%s, format=%s' % (s, format))
1519 if iend > 0 and iend != (len(s)-dotpos-1):
1520 raise FractionalSecondsWrongNumberOfDigits(
1521 'string=%s, format=%s' % (s, format))
1523 format = format[:-len(fixed_endings[iend])]
1524 fracsec = float(s[dotpos:])
1525 s = s[:dotpos]
1527 elif format.endswith('.OPTFRAC'):
1528 dotpos = s.rfind('.')
1529 format = format[:-8]
1530 if dotpos != -1 and len(s[dotpos:]) > 1:
1531 fracsec = float(s[dotpos:])
1533 if dotpos != -1:
1534 s = s[:dotpos]
1536 try:
1537 return time_float(calendar.timegm(time.strptime(s, format))) \
1538 + fracsec
1539 except ValueError as e:
1540 raise TimeStrError('%s, string=%s, format=%s' % (str(e), s, format))
1543stt = str_to_time
1546def time_to_str(t, format='%Y-%m-%d %H:%M:%S.3FRAC'):
1547 '''
1548 Get string representation for floating point system time.
1550 :param t: floating point system time
1551 :param format: time string format
1552 :returns: string representing UTC time
1554 Uses the semantics of :py:func:`time.strftime` but additionally allows for
1555 fractional seconds. If ``format`` contains ``'.xFRAC'``, where ``x`` is a
1556 digit between 1 and 9, this is replaced with the fractional part of ``t``
1557 with ``x`` digits precision.
1558 '''
1560 if pyrocko.grumpy > 0:
1561 check_time_class(t, 'warn' if pyrocko.grumpy == 1 else 'raise')
1563 if isinstance(format, int):
1564 if format > 0:
1565 format = '%Y-%m-%d %H:%M:%S.' + '%iFRAC' % format
1566 else:
1567 format = '%Y-%m-%d %H:%M:%S'
1569 if util_ext is not None:
1570 t0 = num.floor(t)
1571 try:
1572 return util_ext.tts(int(t0), float(t - t0), format)
1573 except util_ext.UtilExtError as e:
1574 raise TimeStrError(
1575 '%s, timestamp=%f, format=%s' % (str(e), t, format))
1577 if not GlobalVars.re_frac:
1578 GlobalVars.re_frac = re.compile(r'\.[1-9]FRAC')
1579 GlobalVars.frac_formats = dict(
1580 [('.%sFRAC' % x, '%.'+x+'f') for x in '123456789'])
1582 ts = float(num.floor(t))
1583 tfrac = t-ts
1585 m = GlobalVars.re_frac.search(format)
1586 if m:
1587 sfrac = (GlobalVars.frac_formats[m.group(0)] % tfrac)
1588 if sfrac[0] == '1':
1589 ts += 1.
1591 format, nsub = GlobalVars.re_frac.subn(sfrac[1:], format, 1)
1593 return time.strftime(format, time.gmtime(ts))
1596tts = time_to_str
1599def mystrftime(fmt=None, tt=None, milliseconds=0):
1601 if fmt is None:
1602 fmt = '%Y-%m-%d %H:%M:%S .%r'
1604 if tt is None:
1605 tt = time.time()
1607 fmt = fmt.replace('%r', '%03i' % int(round(milliseconds)))
1608 fmt = fmt.replace('%u', '%06i' % int(round(milliseconds*1000)))
1609 fmt = fmt.replace('%n', '%09i' % int(round(milliseconds*1000000)))
1610 return time.strftime(fmt, tt)
1613def gmtime_x(timestamp):
1614 etimestamp = float(num.floor(timestamp))
1615 tt = time.gmtime(etimestamp)
1616 ms = (timestamp-etimestamp)*1000
1617 return tt, ms
1620def plural_s(n):
1621 if n == 1:
1622 return ''
1623 else:
1624 return 's'
1627def ensuredirs(dst):
1628 '''
1629 Create all intermediate path components for a target path.
1631 :param dst: target path
1633 The leaf part of the target path is not created (use :py:func:`ensuredir`
1634 if a the target path is a directory to be created).
1635 '''
1637 d, x = os.path.split(dst.rstrip(os.sep))
1638 dirs = []
1639 while d and not os.path.exists(d):
1640 dirs.append(d)
1641 d, x = os.path.split(d)
1643 dirs.reverse()
1645 for d in dirs:
1646 try:
1647 os.mkdir(d)
1648 except OSError as e:
1649 if not e.errno == errno.EEXIST:
1650 raise
1653def ensuredir(dst):
1654 '''
1655 Create directory and all intermediate path components to it as needed.
1657 :param dst: directory name
1659 Nothing is done if the given target already exists.
1660 '''
1662 if os.path.exists(dst):
1663 return
1665 dst.rstrip(os.sep)
1667 ensuredirs(dst)
1668 try:
1669 os.mkdir(dst)
1670 except OSError as e:
1671 if not e.errno == errno.EEXIST:
1672 raise
1675def reuse(x):
1676 '''
1677 Get unique instance of an object.
1679 :param x: hashable object
1680 :returns: reference to x or an equivalent object
1682 Cache object ``x`` in a global dict for reuse, or if x already
1683 is in that dict, return a reference to it.
1684 '''
1686 grs = GlobalVars.reuse_store
1687 if x not in grs:
1688 grs[x] = x
1689 return grs[x]
1692def deuse(x):
1693 grs = GlobalVars.reuse_store
1694 if x in grs:
1695 del grs[x]
1698class Anon(object):
1699 '''
1700 Dict-to-object utility.
1702 Any given arguments are stored as attributes.
1704 Example::
1706 a = Anon(x=1, y=2)
1707 print a.x, a.y
1708 '''
1710 def __init__(self, **dict):
1711 for k in dict:
1712 self.__dict__[k] = dict[k]
1715def select_files(paths, selector=None, regex=None, show_progress=True):
1716 '''
1717 Recursively select files.
1719 :param paths: entry path names
1720 :param selector: callback for conditional inclusion
1721 :param regex: pattern for conditional inclusion
1722 :param show_progress: if True, indicate start and stop of processing
1723 :returns: list of path names
1725 Recursively finds all files under given entry points ``paths``. If
1726 parameter ``regex`` is a regular expression, only files with matching path
1727 names are included. If additionally parameter ``selector`` is given a
1728 callback function, only files for which the callback returns ``True`` are
1729 included. The callback should take a single argument. The callback is
1730 called with a single argument, an object, having as attributes, any named
1731 groups given in ``regex``.
1733 Examples
1735 To find all files ending in ``'.mseed'`` or ``'.msd'``::
1737 select_files(paths,
1738 regex=r'\\.(mseed|msd)$')
1740 To find all files ending with ``'$Year.$DayOfYear'``, having set 2009 for
1741 the year::
1743 select_files(paths,
1744 regex=r'(?P<year>\\d\\d\\d\\d)\\.(?P<doy>\\d\\d\\d)$',
1745 selector=(lambda x: int(x.year) == 2009))
1746 '''
1748 if show_progress:
1749 progress_beg('selecting files...')
1750 if logger.isEnabledFor(logging.DEBUG):
1751 sys.stderr.write('\n')
1753 good = []
1754 if regex:
1755 rselector = re.compile(regex)
1757 def addfile(path):
1758 if regex:
1759 logger.debug("looking at filename: '%s'" % path)
1760 m = rselector.search(path)
1761 if m:
1762 infos = Anon(**m.groupdict())
1763 logger.debug(" regex '%s' matches." % regex)
1764 for k, v in m.groupdict().items():
1765 logger.debug(
1766 " attribute '%s' has value '%s'" % (k, v))
1767 if selector is None or selector(infos):
1768 good.append(os.path.abspath(path))
1770 else:
1771 logger.debug(" regex '%s' does not match." % regex)
1772 else:
1773 good.append(os.path.abspath(path))
1775 if isinstance(paths, str):
1776 paths = [paths]
1778 for path in paths:
1779 if os.path.isdir(path):
1780 for (dirpath, dirnames, filenames) in os.walk(path):
1781 for filename in filenames:
1782 addfile(op.join(dirpath, filename))
1783 else:
1784 addfile(path)
1786 if show_progress:
1787 progress_end('%i file%s selected.' % (len(good), plural_s(len(good))))
1789 return good
1792def base36encode(number, alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
1793 '''
1794 Convert positive integer to a base36 string.
1795 '''
1797 if not isinstance(number, (int, long)):
1798 raise TypeError('number must be an integer')
1799 if number < 0:
1800 raise ValueError('number must be positive')
1802 # Special case for small numbers
1803 if number < 36:
1804 return alphabet[number]
1806 base36 = ''
1807 while number != 0:
1808 number, i = divmod(number, 36)
1809 base36 = alphabet[i] + base36
1811 return base36
1814def base36decode(number):
1815 '''
1816 Decode base36 endcoded positive integer.
1817 '''
1819 return int(number, 36)
1822class UnpackError(Exception):
1823 '''
1824 Exception raised when :py:func:`unpack_fixed` encounters an error.
1825 '''
1827 pass
1830ruler = ''.join(['%-10i' % i for i in range(8)]) \
1831 + '\n' + '0123456789' * 8 + '\n'
1834def unpack_fixed(format, line, *callargs):
1835 '''
1836 Unpack fixed format string, as produced by many fortran codes.
1838 :param format: format specification
1839 :param line: string to be processed
1840 :param callargs: callbacks for callback fields in the format
1842 The format is described by a string of comma-separated fields. Each field
1843 is defined by a character for the field type followed by the field width. A
1844 questionmark may be appended to the field description to allow the argument
1845 to be optional (The data string is then allowed to be filled with blanks
1846 and ``None`` is returned in this case).
1848 The following field types are available:
1850 ==== ================================================================
1851 Type Description
1852 ==== ================================================================
1853 A string (full field width is extracted)
1854 a string (whitespace at the beginning and the end is removed)
1855 i integer value
1856 f floating point value
1857 @ special type, a callback must be given for the conversion
1858 x special field type to skip parts of the string
1859 ==== ================================================================
1860 '''
1862 ipos = 0
1863 values = []
1864 icall = 0
1865 for form in format.split(','):
1866 form = form.strip()
1867 optional = form[-1] == '?'
1868 form = form.rstrip('?')
1869 typ = form[0]
1870 ln = int(form[1:])
1871 s = line[ipos:ipos+ln]
1872 cast = {
1873 'x': None,
1874 'A': str,
1875 'a': lambda x: x.strip(),
1876 'i': int,
1877 'f': float,
1878 '@': 'extra'}[typ]
1880 if cast == 'extra':
1881 cast = callargs[icall]
1882 icall += 1
1884 if cast is not None:
1885 if optional and s.strip() == '':
1886 values.append(None)
1887 else:
1888 try:
1889 values.append(cast(s))
1890 except Exception:
1891 mark = [' '] * 80
1892 mark[ipos:ipos+ln] = ['^'] * ln
1893 mark = ''.join(mark)
1894 raise UnpackError(
1895 'Invalid cast to type "%s" at position [%i:%i] of '
1896 'line: \n%s%s\n%s' % (
1897 typ, ipos, ipos+ln, ruler, line.rstrip(), mark))
1899 ipos += ln
1901 return values
1904_pattern_cache = {}
1907def _nslc_pattern(pattern):
1908 if pattern not in _pattern_cache:
1909 rpattern = re.compile(fnmatch.translate(pattern), re.I)
1910 _pattern_cache[pattern] = rpattern
1911 else:
1912 rpattern = _pattern_cache[pattern]
1914 return rpattern
1917def match_nslc(patterns, nslc):
1918 '''
1919 Match network-station-location-channel code against pattern or list of
1920 patterns.
1922 :param patterns: pattern or list of patterns
1923 :param nslc: tuple with (network, station, location, channel) as strings
1925 :returns: ``True`` if the pattern matches or if any of the given patterns
1926 match; or ``False``.
1928 The patterns may contain shell-style wildcards: \\*, ?, [seq], [!seq].
1930 Example::
1932 match_nslc('*.HAM3.*.BH?', ('GR', 'HAM3', '', 'BHZ')) # -> True
1933 '''
1935 if isinstance(patterns, str):
1936 patterns = [patterns]
1938 if not isinstance(nslc, str):
1939 s = '.'.join(nslc)
1940 else:
1941 s = nslc
1943 for pattern in patterns:
1944 if _nslc_pattern(pattern).match(s):
1945 return True
1947 return False
1950def match_nslcs(patterns, nslcs):
1951 '''
1952 Get network-station-location-channel codes that match given pattern or any
1953 of several given patterns.
1955 :param patterns: pattern or list of patterns
1956 :param nslcs: list of (network, station, location, channel) tuples
1958 See also :py:func:`match_nslc`
1959 '''
1961 matching = []
1962 for nslc in nslcs:
1963 if match_nslc(patterns, nslc):
1964 matching.append(nslc)
1966 return matching
1969class Timeout(Exception):
1970 pass
1973def create_lockfile(fn, timeout=None, timewarn=10.):
1974 t0 = time.time()
1976 while True:
1977 try:
1978 f = os.open(fn, os.O_CREAT | os.O_WRONLY | os.O_EXCL)
1979 os.close(f)
1980 return
1982 except OSError as e:
1983 if e.errno == errno.EEXIST:
1984 tnow = time.time()
1986 if timeout is not None and tnow - t0 > timeout:
1987 raise Timeout(
1988 'Timeout (%gs) occured while waiting to get exclusive '
1989 'access to: %s' % (timeout, fn))
1991 if timewarn is not None and tnow - t0 > timewarn:
1992 logger.warn(
1993 'Waiting since %gs to get exclusive access to: %s' % (
1994 timewarn, fn))
1996 timewarn *= 2
1998 time.sleep(0.01)
1999 else:
2000 raise
2003def delete_lockfile(fn):
2004 os.unlink(fn)
2007class Lockfile(Exception):
2009 def __init__(self, path, timeout=5, timewarn=10.):
2010 self._path = path
2011 self._timeout = timeout
2012 self._timewarn = timewarn
2014 def __enter__(self):
2015 create_lockfile(
2016 self._path, timeout=self._timeout, timewarn=self._timewarn)
2017 return None
2019 def __exit__(self, type, value, traceback):
2020 delete_lockfile(self._path)
2023class SoleError(Exception):
2024 '''
2025 Exception raised by objects of type :py:class:`Sole`, when an concurrent
2026 instance is running.
2027 '''
2029 pass
2032class Sole(object):
2033 '''
2034 Use POSIX advisory file locking to ensure that only a single instance of a
2035 program is running.
2037 :param pid_path: path to lockfile to be used
2039 Usage::
2041 from pyrocko.util import Sole, SoleError, setup_logging
2042 import os
2044 setup_logging('my_program')
2046 pid_path = os.path.join(os.environ['HOME'], '.my_program_lock')
2047 try:
2048 sole = Sole(pid_path)
2050 except SoleError, e:
2051 logger.fatal( str(e) )
2052 sys.exit(1)
2053 '''
2055 def __init__(self, pid_path):
2056 self._pid_path = pid_path
2057 self._other_running = False
2058 ensuredirs(self._pid_path)
2059 self._lockfile = None
2060 self._os = os
2062 if not fcntl:
2063 raise SoleError(
2064 'Python standard library module "fcntl" not available on '
2065 'this platform.')
2067 self._fcntl = fcntl
2069 try:
2070 self._lockfile = os.open(self._pid_path, os.O_CREAT | os.O_WRONLY)
2071 except Exception:
2072 raise SoleError(
2073 'Cannot open lockfile (path = %s)' % self._pid_path)
2075 try:
2076 fcntl.lockf(self._lockfile, fcntl.LOCK_EX | fcntl.LOCK_NB)
2078 except IOError:
2079 self._other_running = True
2080 try:
2081 f = open(self._pid_path, 'r')
2082 pid = f.read().strip()
2083 f.close()
2084 except Exception:
2085 pid = '?'
2087 raise SoleError('Other instance is running (pid = %s)' % pid)
2089 try:
2090 os.ftruncate(self._lockfile, 0)
2091 os.write(self._lockfile, '%i\n' % os.getpid())
2092 os.fsync(self._lockfile)
2094 except Exception:
2095 # the pid is only stored for user information, so this is allowed
2096 # to fail
2097 pass
2099 def __del__(self):
2100 if not self._other_running:
2101 if self._lockfile is not None:
2102 self._fcntl.lockf(self._lockfile, self._fcntl.LOCK_UN)
2103 self._os.close(self._lockfile)
2104 try:
2105 self._os.unlink(self._pid_path)
2106 except Exception:
2107 pass
2110re_escapequotes = re.compile(r"(['\\])")
2113def escapequotes(s):
2114 return re_escapequotes.sub(r"\\\1", s)
2117class TableWriter(object):
2118 '''
2119 Write table of space separated values to a file.
2121 :param f: file like object
2123 Strings containing spaces are quoted on output.
2124 '''
2126 def __init__(self, f):
2127 self._f = f
2129 def writerow(self, row, minfieldwidths=None):
2131 '''
2132 Write one row of values to underlying file.
2134 :param row: iterable of values
2135 :param minfieldwidths: minimum field widths for the values
2137 Each value in in ``row`` is converted to a string and optionally padded
2138 with blanks. The resulting strings are output separated with blanks. If
2139 any values given are strings and if they contain whitespace, they are
2140 quoted with single quotes, and any internal single quotes are
2141 backslash-escaped.
2142 '''
2144 out = []
2146 for i, x in enumerate(row):
2147 w = 0
2148 if minfieldwidths and i < len(minfieldwidths):
2149 w = minfieldwidths[i]
2151 if isinstance(x, str):
2152 if re.search(r"\s|'", x):
2153 x = "'%s'" % escapequotes(x)
2155 x = x.ljust(w)
2156 else:
2157 x = str(x).rjust(w)
2159 out.append(x)
2161 self._f.write(' '.join(out).rstrip() + '\n')
2164class TableReader(object):
2166 '''
2167 Read table of space separated values from a file.
2169 :param f: file-like object
2171 This uses Pythons shlex module to tokenize lines. Should deal correctly
2172 with quoted strings.
2173 '''
2175 def __init__(self, f):
2176 self._f = f
2177 self.eof = False
2179 def readrow(self):
2180 '''
2181 Read one row from the underlying file, tokenize it with shlex.
2183 :returns: tokenized line as a list of strings.
2184 '''
2186 line = self._f.readline()
2187 if not line:
2188 self.eof = True
2189 return []
2190 line.strip()
2191 if line.startswith('#'):
2192 return []
2194 return qsplit(line)
2197def gform(number, significant_digits=3):
2198 '''
2199 Pretty print floating point numbers.
2201 Align floating point numbers at the decimal dot.
2203 ::
2205 | -d.dde+xxx|
2206 | -d.dde+xx |
2207 |-ddd. |
2208 | -dd.d |
2209 | -d.dd |
2210 | -0.ddd |
2211 | -0.0ddd |
2212 | -0.00ddd |
2213 | -d.dde-xx |
2214 | -d.dde-xxx|
2215 | nan|
2218 The formatted string has length ``significant_digits * 2 + 6``.
2219 '''
2221 no_exp_range = (pow(10., -1),
2222 pow(10., significant_digits))
2223 width = significant_digits+significant_digits-1+1+1+5
2225 if (no_exp_range[0] <= abs(number) < no_exp_range[1]) or number == 0.:
2226 s = ('%#.*g' % (significant_digits, number)).rstrip('0')
2227 else:
2228 s = '%.*E' % (significant_digits-1, number)
2229 s = (' '*(-s.find('.')+(significant_digits+1))+s).ljust(width)
2230 if s.strip().lower() == 'nan':
2231 s = 'nan'.rjust(width)
2232 return s
2235def human_bytesize(value):
2237 exts = 'Bytes kB MB GB TB PB EB ZB YB'.split()
2239 if value == 1:
2240 return '1 Byte'
2242 for i, ext in enumerate(exts):
2243 x = float(value) / 1000**i
2244 if round(x) < 10. and not value < 1000:
2245 return '%.1f %s' % (x, ext)
2246 if round(x) < 1000.:
2247 return '%.0f %s' % (x, ext)
2249 return '%i Bytes' % value
2252re_compatibility = re.compile(
2253 r'!pyrocko\.(trace|gf\.(meta|seismosizer)|fomosto\.' +
2254 r'(dummy|poel|qseis|qssp))\.'
2255)
2258def pf_is_old(fn):
2259 oldstyle = False
2260 with open(fn, 'r') as f:
2261 for line in f:
2262 if re_compatibility.search(line):
2263 oldstyle = True
2265 return oldstyle
2268def pf_upgrade(fn):
2269 need = pf_is_old(fn)
2270 if need:
2271 fn_temp = fn + '.temp'
2273 with open(fn, 'r') as fin:
2274 with open(fn_temp, 'w') as fout:
2275 for line in fin:
2276 line = re_compatibility.sub('!pf.', line)
2277 fout.write(line)
2279 os.rename(fn_temp, fn)
2281 return need
2284def read_leap_seconds(tzfile='/usr/share/zoneinfo/right/UTC'):
2285 '''
2286 Extract leap second information from tzdata.
2288 Based on example at http://stackoverflow.com/questions/19332902/\
2289 extract-historic-leap-seconds-from-tzdata
2291 See also 'man 5 tzfile'.
2292 '''
2294 from struct import unpack, calcsize
2295 out = []
2296 with open(tzfile, 'rb') as f:
2297 # read header
2298 fmt = '>4s c 15x 6l'
2299 (magic, format, ttisgmtcnt, ttisstdcnt, leapcnt, timecnt,
2300 typecnt, charcnt) = unpack(fmt, f.read(calcsize(fmt)))
2301 assert magic == 'TZif'.encode('US-ASCII'), 'Not a timezone file'
2303 # skip over some uninteresting data
2304 fmt = '>%(timecnt)dl %(timecnt)dB %(ttinfo)s %(charcnt)ds' % dict(
2305 timecnt=timecnt, ttinfo='lBB'*typecnt, charcnt=charcnt)
2306 f.read(calcsize(fmt))
2308 # read leap-seconds
2309 fmt = '>2l'
2310 for i in range(leapcnt):
2311 tleap, nleap = unpack(fmt, f.read(calcsize(fmt)))
2312 out.append((tleap-nleap+1, nleap))
2314 return out
2317class LeapSecondsError(Exception):
2318 pass
2321class LeapSecondsOutdated(LeapSecondsError):
2322 pass
2325class InvalidLeapSecondsFile(LeapSecondsOutdated):
2326 pass
2329def parse_leap_seconds_list(fn):
2330 data = []
2331 texpires = None
2332 try:
2333 t0 = int(round(str_to_time('1900-01-01 00:00:00')))
2334 except TimeStrError:
2335 t0 = int(round(str_to_time('1970-01-01 00:00:00'))) - 2208988800
2337 tnow = int(round(time.time()))
2339 if not op.exists(fn):
2340 raise LeapSecondsOutdated('no leap seconds file found')
2342 try:
2343 with open(fn, 'rb') as f:
2344 for line in f:
2345 if line.strip().startswith(b'<!DOCTYPE'):
2346 raise InvalidLeapSecondsFile('invalid leap seconds file')
2348 if line.startswith(b'#@'):
2349 texpires = int(line.split()[1]) + t0
2350 elif line.startswith(b'#') or len(line) < 5:
2351 pass
2352 else:
2353 toks = line.split()
2354 t = int(toks[0]) + t0
2355 nleap = int(toks[1]) - 10
2356 data.append((t, nleap))
2358 except IOError:
2359 raise LeapSecondsError('cannot read leap seconds file %s' % fn)
2361 if texpires is None or tnow > texpires:
2362 raise LeapSecondsOutdated('leap seconds list is outdated')
2364 return data
2367def read_leap_seconds2():
2368 from pyrocko import config
2369 conf = config.config()
2370 fn = conf.leapseconds_path
2371 url = conf.leapseconds_url
2372 # check for outdated default URL
2373 if url == 'http://www.ietf.org/timezones/data/leap-seconds.list':
2374 url = 'https://www.ietf.org/timezones/data/leap-seconds.list'
2375 logger.info(
2376 'Leap seconds default URL is now: %s\nUsing new default.' % url)
2378 for i in range(3):
2379 try:
2380 return parse_leap_seconds_list(fn)
2382 except LeapSecondsOutdated:
2383 try:
2384 logger.info('updating leap seconds list...')
2385 download_file(url, fn)
2387 except Exception as e:
2388 raise LeapSecondsError(
2389 'cannot download leap seconds list from %s to %s (%s)'
2390 % (url, fn, e))
2392 raise LeapSecondsError('Could not retrieve/read leap seconds file.')
2395def gps_utc_offset(t_utc):
2396 '''
2397 Time offset t_gps - t_utc for a given t_utc.
2398 '''
2399 ls = read_leap_seconds2()
2400 i = 0
2401 if t_utc < ls[0][0]:
2402 return ls[0][1] - 1 - 9
2404 while i < len(ls) - 1:
2405 if ls[i][0] <= t_utc and t_utc < ls[i+1][0]:
2406 return ls[i][1] - 9
2407 i += 1
2409 return ls[-1][1] - 9
2412def utc_gps_offset(t_gps):
2413 '''
2414 Time offset t_utc - t_gps for a given t_gps.
2415 '''
2416 ls = read_leap_seconds2()
2418 if t_gps < ls[0][0] + ls[0][1] - 9:
2419 return - (ls[0][1] - 1 - 9)
2421 i = 0
2422 while i < len(ls) - 1:
2423 if ls[i][0] + ls[i][1] - 9 <= t_gps \
2424 and t_gps < ls[i+1][0] + ls[i+1][1] - 9:
2425 return - (ls[i][1] - 9)
2426 i += 1
2428 return - (ls[-1][1] - 9)
2431def make_iload_family(iload_fh, doc_fmt='FMT', doc_yielded_objects='FMT'):
2432 import itertools
2433 import glob
2434 from pyrocko.io.io_common import FileLoadError
2436 def iload_filename(filename, **kwargs):
2437 try:
2438 with open(filename, 'rb') as f:
2439 for cr in iload_fh(f, **kwargs):
2440 yield cr
2442 except FileLoadError as e:
2443 e.set_context('filename', filename)
2444 raise
2446 def iload_dirname(dirname, **kwargs):
2447 for entry in os.listdir(dirname):
2448 fpath = op.join(dirname, entry)
2449 if op.isfile(fpath):
2450 for cr in iload_filename(fpath, **kwargs):
2451 yield cr
2453 def iload_glob(pattern, **kwargs):
2455 for fn in glob.iglob(pattern):
2456 for cr in iload_filename(fn, **kwargs):
2457 yield cr
2459 def iload(source, **kwargs):
2460 if isinstance(source, str):
2461 if op.isdir(source):
2462 return iload_dirname(source, **kwargs)
2463 elif op.isfile(source):
2464 return iload_filename(source, **kwargs)
2465 else:
2466 return iload_glob(source, **kwargs)
2468 elif hasattr(source, 'read'):
2469 return iload_fh(source, **kwargs)
2470 else:
2471 return itertools.chain.from_iterable(
2472 iload(subsource, **kwargs) for subsource in source)
2474 iload_filename.__doc__ = '''
2475 Read %s information from named file.
2476 ''' % doc_fmt
2478 iload_dirname.__doc__ = '''
2479 Read %s information from directory of %s files.
2480 ''' % (doc_fmt, doc_fmt)
2482 iload_glob.__doc__ = '''
2483 Read %s information from files matching a glob pattern.
2484 ''' % doc_fmt
2486 iload.__doc__ = '''
2487 Load %s information from given source(s)
2489 The ``source`` can be specified as the name of a %s file, the name of a
2490 directory containing %s files, a glob pattern of %s files, an open
2491 filehandle or an iterator yielding any of the forementioned sources.
2493 This function behaves as a generator yielding %s objects.
2494 ''' % (doc_fmt, doc_fmt, doc_fmt, doc_fmt, doc_yielded_objects)
2496 for f in iload_filename, iload_dirname, iload_glob, iload:
2497 f.__module__ = iload_fh.__module__
2499 return iload_filename, iload_dirname, iload_glob, iload
2502class Inconsistency(Exception):
2503 pass
2506def consistency_check(list_of_tuples, message='values differ:'):
2507 '''
2508 Check for inconsistencies.
2510 Given a list of tuples, check that all tuple elements except for first one
2511 match. E.g. ``[('STA.N', 55.3, 103.2), ('STA.E', 55.3, 103.2)]`` would be
2512 valid because the coordinates at the two channels are the same.
2513 '''
2515 if len(list_of_tuples) >= 2:
2516 if any(t[1:] != list_of_tuples[0][1:] for t in list_of_tuples[1:]):
2517 raise Inconsistency('%s\n' % message + '\n'.join(
2518 ' %s: %s' % (t[0], ', '.join('%g' % x for x in t[1:]))
2519 for t in list_of_tuples))
2522class defaultzerodict(dict):
2523 def __missing__(self, k):
2524 return 0
2527def mostfrequent(x):
2528 c = defaultzerodict()
2529 for e in x:
2530 c[e] += 1
2532 return sorted(list(c.keys()), key=lambda k: c[k])[-1]
2535def consistency_merge(list_of_tuples,
2536 message='values differ:',
2537 error='raise',
2538 merge=mostfrequent):
2540 assert error in ('raise', 'warn', 'ignore')
2542 if len(list_of_tuples) == 0:
2543 raise Exception('cannot merge empty sequence')
2545 try:
2546 consistency_check(list_of_tuples, message)
2547 return list_of_tuples[0][1:]
2548 except Inconsistency as e:
2549 if error == 'raise':
2550 raise
2552 elif error == 'warn':
2553 logger.warning(str(e))
2555 return tuple([merge(x) for x in list(zip(*list_of_tuples))[1:]])
2558def parse_md(f):
2559 try:
2560 with open(op.join(
2561 op.dirname(op.abspath(f)),
2562 'README.md'), 'r') as readme:
2563 mdstr = readme.read()
2564 except IOError as e:
2565 return 'Failed to get README.md: %s' % e
2567 # Remve the title
2568 mdstr = re.sub(r'^# .*\n?', '', mdstr)
2569 # Append sphinx reference to `pyrocko.` modules
2570 mdstr = re.sub(r'`pyrocko\.(.*)`', r':py:mod:`pyrocko.\1`', mdstr)
2571 # Convert Subsections to toc-less rubrics
2572 mdstr = re.sub(r'## (.*)\n', r'.. rubric:: \1\n', mdstr)
2573 return mdstr
2576def mpl_show(plt):
2577 import matplotlib
2578 if matplotlib.get_backend().lower() == 'agg':
2579 logger.warning('Cannot show() when using matplotlib "agg" backend')
2580 else:
2581 plt.show()
2584g_re_qsplit = re.compile(
2585 r'"([^"\\]*(?:\\.[^"\\]*)*)"|\'([^\'\\]*(?:\\.[^\'\\]*)*)\'|(\S+)')
2588g_re_trivial = re.compile(r'\A[^\'"\s]+\Z')
2589g_re_escape_s = re.compile(r'([\\\'])')
2590g_re_unescape_s = re.compile(r'\\([\\\'])')
2591g_re_escape_d = re.compile(r'([\\"])')
2592g_re_unescape_d = re.compile(r'\\([\\"])')
2595def escape_s(s):
2596 '''
2597 Backslash-escape single-quotes and backslashes.
2599 Example: ``Jack's`` => ``Jack\\'s``
2601 '''
2602 return g_re_escape_s.sub(r'\\\1', s)
2605def unescape_s(s):
2606 '''
2607 Unescape backslash-escaped single-quotes and backslashes.
2609 Example: ``Jack\\'s`` => ``Jack's``
2610 '''
2611 return g_re_unescape_s.sub(r'\1', s)
2614def escape_d(s):
2615 '''
2616 Backslash-escape double-quotes and backslashes.
2618 Example: ``"Hello \\O/"`` => ``\\"Hello \\\\O/\\"``
2619 '''
2620 return g_re_escape_d.sub(r'\\\1', s)
2623def unescape_d(s):
2624 '''
2625 Unescape backslash-escaped double-quotes and backslashes.
2627 Example: ``\\"Hello \\\\O/\\"`` => ``"Hello \\O/"``
2628 '''
2629 return g_re_unescape_d.sub(r'\1', s)
2632def qjoin_s(it):
2633 '''
2634 Join sequence of strings into a line, single-quoting non-trivial strings.
2636 Example: ``["55", "Sparrow's Island"]`` => ``"55 'Sparrow\\\\'s Island'"``
2637 '''
2638 return ' '.join(
2639 w if g_re_trivial.search(w) else "'%s'" % escape_s(w) for w in it)
2642def qjoin_d(it):
2643 '''
2644 Join sequence of strings into a line, double-quoting non-trivial strings.
2646 Example: ``['55', 'Pete "The Robot" Smith']`` =>
2647 ``'55' "Pete \\\\"The Robot\\\\" Smith"'``
2648 '''
2649 return ' '.join(
2650 w if g_re_trivial.search(w) else '"%s"' % escape_d(w) for w in it)
2653def qsplit(s):
2654 '''
2655 Split line into list of strings, allowing for quoted strings.
2657 Example: ``"55 'Sparrow\\\\'s Island'"`` =>
2658 ``["55", "Sparrow's Island"]``,
2659 ``'55' "Pete \\\\"The Robot\\\\" Smith"'`` =>
2660 ``['55', 'Pete "The Robot" Smith']``
2661 '''
2662 return [
2663 (unescape_d(x[0]) or unescape_s(x[1]) or x[2])
2664 for x in g_re_qsplit.findall(s)]
2667g_have_warned_threadpoolctl = False
2670class threadpool_limits_dummy(object):
2672 def __init__(self, *args, **kwargs):
2673 pass
2675 def __enter__(self):
2676 global g_have_warned_threadpoolctl
2678 if not g_have_warned_threadpoolctl:
2679 logger.warning(
2680 'Cannot control number of BLAS threads because '
2681 '`threadpoolctl` module is not available. You may want to '
2682 'install `threadpoolctl`.')
2684 g_have_warned_threadpoolctl = True
2686 return self
2688 def __exit__(self, type, value, traceback):
2689 pass
2692def get_threadpool_limits():
2693 '''
2694 Try to import threadpoolctl.threadpool_limits, provide dummy if not avail.
2695 '''
2697 try:
2698 from threadpoolctl import threadpool_limits
2699 return threadpool_limits
2701 except ImportError:
2702 return threadpool_limits_dummy