Commit 63396a6c authored by Marius Kriegerowski's avatar Marius Kriegerowski
Browse files

Merge branch 'master' of gitext.gfz-potsdam.de:heimann/lassie

parents 0fd6549e 91a98533
import logging import logging
import os.path as op import os.path as op
from pyrocko.guts import String, Float, Timestamp, List, Bool from pyrocko.guts import String, Float, Timestamp, List, Bool, Int
from pyrocko import model, guts from pyrocko import model, guts
from pyrocko.io import stationxml from pyrocko.io import stationxml
from pyrocko.gf import TPDef from pyrocko.gf import TPDef
...@@ -120,6 +120,18 @@ class Config(HasPaths): ...@@ -120,6 +120,18 @@ class Config(HasPaths):
default='lassie_phases.cache', default='lassie_phases.cache',
help='directory where lassie stores tabulated phases etc.') help='directory where lassie stores tabulated phases etc.')
stacking_blocksize = Int.T(
optional=True,
help='enable chunked stacking to reduce memory usage. Setting this to '
'e.g. 64 will use ngridpoints * 64 * 8 bytes of memory to hold '
'the stacking results, instead of computing the whole processing '
'time window in one shot. Setting this to a very small number '
'may lead to bad performance. If this is enabled together with '
'plotting, the cutout of the image function seen in the map '
'image must be stacked again just for plotting (redundantly and '
'memory greedy) because it may intersect more than one '
'processing chunk.')
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
HasPaths.__init__(self, *args, **kwargs) HasPaths.__init__(self, *args, **kwargs)
self._receivers = None self._receivers = None
...@@ -131,8 +143,8 @@ class Config(HasPaths): ...@@ -131,8 +143,8 @@ class Config(HasPaths):
''' '''
Post-init setup of image function contributors. Post-init setup of image function contributors.
''' '''
for ifc in self.image_function_contributions: for ifc_ in self.image_function_contributions:
ifc.setup(self) ifc_.setup(self)
def set_config_name(self, config_name): def set_config_name(self, config_name):
self._config_name = config_name self._config_name = config_name
......
...@@ -2,6 +2,7 @@ import math ...@@ -2,6 +2,7 @@ import math
import logging import logging
import os.path as op import os.path as op
import shutil import shutil
import time
from collections import defaultdict from collections import defaultdict
...@@ -309,10 +310,14 @@ def search( ...@@ -309,10 +310,14 @@ def search(
if config.fill_incomplete_with_zeros: if config.fill_incomplete_with_zeros:
trs = zero_fill(trs, wmin - tpad, wmax + tpad) trs = zero_fill(trs, wmin - tpad, wmax + tpad)
frames = None t0 = math.floor(wmin / deltat_cf) * deltat_cf
iwmin = int(round((wmin-tpeaksearch-t0) / deltat_cf))
iwmax = int(round((wmax+tpeaksearch-t0) / deltat_cf))
lengthout = iwmax - iwmin + 1
pdata = [] pdata = []
trs_debug = [] trs_debug = []
shift_maxs = [] parstack_params = []
for iifc, ifc in enumerate(ifcs): for iifc, ifc in enumerate(ifcs):
dataset = ifc.preprocess( dataset = ifc.preprocess(
trs, wmin-tpeaksearch, wmax+tpeaksearch, trs, wmin-tpeaksearch, wmax+tpeaksearch,
...@@ -329,8 +334,6 @@ def search( ...@@ -329,8 +334,6 @@ def search(
trs_debug.extend(trs + list(trs_selected)) trs_debug.extend(trs + list(trs_selected))
t0 = (wmin / deltat_cf) * deltat_cf
istations_selected = num.array( istations_selected = num.array(
[station_index[nsl] for nsl in nsls_selected], [station_index[nsl] for nsl in nsls_selected],
dtype=num.int) dtype=num.int)
...@@ -362,35 +365,58 @@ def search( ...@@ -362,35 +365,58 @@ def search(
pdata.append((list(trs_selected), shift_table, ifc)) pdata.append((list(trs_selected), shift_table, ifc))
iwmin = int(round((wmin-tpeaksearch-t0) / deltat_cf)) parstack_params.append((arrays, offsets, shifts, weights))
iwmax = int(round((wmax+tpeaksearch-t0) / deltat_cf))
lengthout = iwmax - iwmin + 1 if config.stacking_blocksize is not None:
ipstep = config.stacking_blocksize
shift_maxs.append(num.max(-shifts) * deltat_cf) frames = None
else:
ipstep = lengthout
frames = num.zeros((ngridpoints, lengthout))
twall_start = time.time()
frame_maxs = num.zeros(lengthout)
frame_argmaxs = num.zeros(lengthout, dtype=num.int)
ipmin = iwmin
while ipmin < iwmin+lengthout:
ipsize = min(ipstep, iwmin + lengthout - ipmin)
if ipstep == lengthout:
frames_p = frames
else:
frames_p = num.zeros((ngridpoints, ipsize))
frames, ioff = parstack( for (arrays, offsets, shifts, weights) in parstack_params:
frames_p, _ = parstack(
arrays, offsets, shifts, weights, 0, arrays, offsets, shifts, weights, 0,
offsetout=iwmin, offsetout=ipmin,
lengthout=lengthout, lengthout=ipsize,
result=frames, result=frames_p,
nparallel=nparallel, nparallel=nparallel,
impl='openmp') impl='openmp')
shift_max = max(shift_maxs)
if config.sharpness_normalization: if config.sharpness_normalization:
frame_maxs = frames.max(axis=0) frame_p_maxs = frames_p.max(axis=0)
frame_means = num.abs(frames).mean(axis=0) frame_p_means = num.abs(frames_p).mean(axis=0)
frames *= (frame_maxs / frame_means)[num.newaxis, :] frames_p *= (frame_p_maxs / frame_p_means)[num.newaxis, :]
frames *= norm_map[:, num.newaxis] frames_p *= norm_map[:, num.newaxis]
if config.ifc_count_normalization: if config.ifc_count_normalization:
frames *= 1.0 / len(ifcs) frames_p *= 1.0 / len(ifcs)
frame_maxs[ipmin-iwmin:ipmin-iwmin+ipsize] = \
frames_p.max(axis=0)
frame_argmaxs[ipmin-iwmin:ipmin-iwmin+ipsize] = \
pargmax(frames_p)
frame_maxs = frames.max(axis=0) ipmin += ipstep
del frames_p
tmin_frames = t0 + ioff * deltat_cf twall_end = time.time()
logger.info('wallclock time for stacking: %g s' % (
twall_end - twall_start))
tmin_frames = t0 + iwmin * deltat_cf
tr_stackmax = trace.Trace( tr_stackmax = trace.Trace(
'', 'SMAX', '', '', '', 'SMAX', '', '',
...@@ -418,17 +444,13 @@ def search( ...@@ -418,17 +444,13 @@ def search(
wmin <= tpeak and tpeak < wmax])) or ([], []) wmin <= tpeak and tpeak < wmax])) or ([], [])
tr_stackmax_indx = tr_stackmax.copy(data=False) tr_stackmax_indx = tr_stackmax.copy(data=False)
tr_stackmax_indx.set_ydata(frame_argmaxs.astype(num.int32))
imaxs = pargmax(frames)
tr_stackmax_indx.set_ydata(imaxs.astype(num.int32))
tr_stackmax_indx.set_location('i') tr_stackmax_indx.set_location('i')
for (tpeak, apeak) in zip(tpeaks, apeaks): for (tpeak, apeak) in zip(tpeaks, apeaks):
iframe = int(round(((tpeak-t0) - ioff*deltat_cf) / deltat_cf)) iframe = int(round((tpeak - tmin_frames) / deltat_cf))
frame = frames[:, iframe] imax = frame_argmaxs[iframe]
imax = imaxs[iframe]
latpeak, lonpeak, xpeak, ypeak, zpeak = \ latpeak, lonpeak, xpeak, ypeak, zpeak = \
grid.index_to_location(imax) grid.index_to_location(imax)
...@@ -477,23 +499,60 @@ def search( ...@@ -477,23 +499,60 @@ def search(
util.ensuredirs(fn) util.ensuredirs(fn)
try: if frames is not None:
frames_p = frames
tmin_frames_p = tmin_frames
iframe_p = iframe
else:
iframe_min = max(
0,
int(round(iframe - tpeaksearch/deltat_cf)))
iframe_max = min(
lengthout-1,
int(round(iframe + tpeaksearch/deltat_cf)))
ipsize = iframe_max - iframe_min + 1
frames_p = num.zeros((ngridpoints, ipsize))
tmin_frames_p = tmin_frames + iframe_min*deltat_cf
iframe_p = iframe - iframe_min
for (arrays, offsets, shifts, weights) \
in parstack_params:
frames_p, _ = parstack(
arrays, offsets, shifts, weights, 0,
offsetout=iwmin+iframe_min,
lengthout=ipsize,
result=frames_p,
nparallel=nparallel,
impl='openmp')
if config.sharpness_normalization:
frame_p_maxs = frames_p.max(axis=0)
frame_p_means = num.abs(frames_p).mean(axis=0)
frames_p *= (
frame_p_maxs/frame_p_means)[num.newaxis, :]
frames_p *= norm_map[:, num.newaxis]
if config.ifc_count_normalization:
frames_p *= 1.0 / len(ifcs)
plot.plot_detection( plot.plot_detection(
grid, receivers, frames, tmin_frames, grid, receivers, frames_p, tmin_frames_p,
deltat_cf, imax, iframe, fsmooth_min, xpeak, ypeak, deltat_cf, imax, iframe_p, xpeak, ypeak,
zpeak, zpeak,
tr_stackmax, tpeaks, apeaks, tr_stackmax, tpeaks, apeaks,
config.detector_threshold, config.detector_threshold,
wmin, wmax, wmin, wmax,
pdata, trs, fmin, fmax, idetection, pdata, trs, fmin, fmax, idetection,
grid_station_shift_max=shift_max, tpeaksearch,
movie=show_movie, movie=show_movie,
show=show_detections, show=show_detections,
save_filename=fn) save_filename=fn)
except AttributeError as e:
logger.warn(e)
del frame del frames_p
if stop_after_first: if stop_after_first:
return return
......
...@@ -108,9 +108,9 @@ def plot_geometry_carthesian(grid, receivers): ...@@ -108,9 +108,9 @@ def plot_geometry_carthesian(grid, receivers):
def plot_detection( def plot_detection(
grid, receivers, frames, tmin_frames, deltat_cf, imax, iframe, grid, receivers, frames, tmin_frames, deltat_cf, imax, iframe,
fsmooth_min, xpeak, ypeak, zpeak, tr_stackmax, tpeaks, apeaks, xpeak, ypeak, zpeak, tr_stackmax, tpeaks, apeaks,
detector_threshold, wmin, wmax, pdata, trs_raw, fmin, fmax, detector_threshold, wmin, wmax, pdata, trs_raw, fmin, fmax,
idetection, grid_station_shift_max, idetection, tpeaksearch,
movie=False, save_filename=None, show=True): movie=False, save_filename=None, show=True):
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
...@@ -123,11 +123,10 @@ def plot_detection( ...@@ -123,11 +123,10 @@ def plot_detection(
distances = grid.distances(receivers) distances = grid.distances(receivers)
plot.mpl_init() plot.mpl_init(fontsize=9)
fig = plt.figure(figsize=plot.mpl_papersize('a4', 'landscape')) fig = plt.figure(figsize=plot.mpl_papersize('a4', 'landscape'))
axes = plt.subplot2grid((2, 3), (0, 2), aspect=1.0) axes = plt.subplot2grid((2, 3), (0, 2), aspect=1.0)
plot.mpl_labelspace(axes) plot.mpl_labelspace(axes)
...@@ -161,7 +160,7 @@ def plot_detection( ...@@ -161,7 +160,7 @@ def plot_detection(
tpeak_current = tmin_frames + deltat_cf * iframe tpeak_current = tmin_frames + deltat_cf * iframe
t0 = tpeak_current t0 = tpeak_current
tduration = 2.0*grid_station_shift_max + 1./fsmooth_min tduration = 2.0*tpeaksearch
axes2.axvspan( axes2.axvspan(
tr_stackmax.tmin - t0, wmin - t0, tr_stackmax.tmin - t0, wmin - t0,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment