python based data analysis

data/tools/2d_event_writer.py

+import ConfigParser as cp
+import os
+import numpy as np
+import csv
+import h5py
+
+from support.maximum import maximum
+from support.norm import norm
+from support.find_quakes import find_quakes
+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+from support.max_distance import max_distance
+
+
+NBODIES = 2
+FINAL_TIME = 1000  # s
+FINAL_VELOCITY = 1e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 1
+
+# read config ini
+config = cp.ConfigParser()
+config_path = os.path.join('config.ini')
+config.read(config_path)
+sim_path = config.get('directories', 'simulation')
+exp_path = config.get('directories', 'experiment')
+out_path = config.get('directories', 'output')
+
+# create output directory
+out_dir = os.path.join(out_path)
+if not os.path.exists(out_dir):
+    os.mkdir(out_dir)
+
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+# read time
+relative_time = np.array(h5file['relativeTime'])
+real_time = relative_time * FINAL_TIME
+
+for body_ID in range(NBODIES):
+    body = 'body' + str(body_ID)
+
+    if body not in h5file:
+        continue
+
+    # read data
+    coordinates = np.array(h5file[body + '/coordinates'])
+    displacement = np.array(h5file[body + '/displacement'])
+    velocity = np.array(h5file[body + '/velocity'])
+
+    num_vertices = displacement.shape[1]
+    velocity_norm = norm(velocity[:, :, TANGENTIAL_COORDS:])
+    maximum_velocity = maximum(velocity_norm)
+
+    [quake_starts, quake_ends] = find_quakes(THRESHOLD_VELOCITY,
+                                             maximum_velocity)
+
+    quakes = {}
+    for quake_ID in range(len(quake_starts)):
+        quake_start = int(quake_starts[quake_ID])
+        quake_end = int(quake_ends[quake_ID])
+
+        quake = {}
+        quake['start'] = real_time[quake_start]
+        quake['end'] = real_time[quake_end]
+
+        local_slipping_times = velocity_norm[quake_start:quake_end, :] \
+            > THRESHOLD_VELOCITY
+        quake_displacement = displacement[quake_start:quake_end, :, :]
+        slip = np.zeros(num_vertices)
+
+        for i in range(num_vertices):
+            if any(local_slipping_times[:, i]):
+                starts = slip_beginnings(local_slipping_times[:, i])
+                ends = slip_endings(local_slipping_times[:, i])
+                slip[i] = np.linalg.norm(quake_displacement[ends, i, :]
+                                         - quake_displacement[starts, i, :])
+
+        quake['peakSlip'] = max(slip)
+        quake['peakSlipRate'] = max(maximum_velocity[quake_start:quake_end])
+
+        max_rupture_width = 0
+        for time_step in range(local_slipping_times.shape[0]):
+            slipping_time = local_slipping_times[time_step, :]
+            if not any(slipping_time):
+                continue
+
+            slipping_coords = coordinates[slipping_time] \
+                + displacement[quake_start + time_step, slipping_time]
+            if len(slipping_coords) > 1:
+                pair = np.array(max_distance(slipping_coords))
+                max_rupture_width = max(max_rupture_width,
+                                        np.linalg.norm(pair[0] - pair[1]))
+
+        quake['ruptureWidth'] = max_rupture_width
+        quakes[quake_ID] = quake
+
+    # output quake data to csv file
+    csv_columns = quakes[0].keys()
+    csv_file_path = os.path.join(out_path, 'events' + str(body_ID) + '.csv')
+    try:
+        with open(csv_file_path, 'w') as csv_file:
+            writer = csv.DictWriter(csv_file, fieldnames=csv_columns)
+            writer.writeheader()
+            for quake_ID in quakes:
+                writer.writerow(quakes[quake_ID])
+    except IOError:
+        print('I/O error')
+
+h5file.close()

data/tools/2d_event_writer.py~

+import ConfigParser as cp
+import os
+import numpy as np
+import csv
+import h5py
+
+from support.maximum import maximum
+from support.norm import norm
+from support.find_quakes import find_quakes
+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+from support.max_distance import max_distance
+
+
+NBODIES = 2
+FINAL_TIME = 1000  # s
+FINAL_VELOCITY = 1e-5  # m/s
+THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY  # 1000e-6 + FINAL_VELOCITY
+
+TANGENTIAL_COORDS = 1
+
+# read config ini
+config = cp.ConfigParser()
+config_path = os.path.join('config.ini')
+config.read(config_path)
+sim_path = config.get('directories', 'simulation')
+exp_path = config.get('directories', 'experiment')
+out_path = config.get('directories', 'output')
+
+# create output directory
+out_dir = os.path.join(out_path)
+if not os.path.exists(out_dir):
+    os.mkdir(out_dir)
+
+h5path = os.path.join(sim_path)
+h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')
+
+# read time
+relative_time = np.array(h5file['relativeTime'])
+real_time = relative_time * FINAL_TIME
+
+for body_ID in range(NBODIES):
+    body = 'body' + str(body_ID)
+
+    if body not in h5file:
+        continue
+
+    # read data
+    coordinates = np.array(h5file[body + '/coordinates'])
+    displacement = np.array(h5file[body + '/displacement'])
+    velocity = np.array(h5file[body + '/velocity'])
+
+    num_vertices = displacement.shape[1]
+    velocity_norm = norm(velocity[:, :, TANGENTIAL_COORDS:])
+    maximum_velocity = maximum(velocity_norm)
+
+    [quake_starts, quake_ends] = find_quakes(THRESHOLD_VELOCITY,
+                                             maximum_velocity)
+
+    quakes = {}
+    for quake_ID in range(len(quake_starts)):
+        quake_start = int(quake_starts[quake_ID])
+        quake_end = int(quake_ends[quake_ID])
+
+        quake = {}
+        quake['start'] = real_time[quake_start]
+        quake['end'] = real_time[quake_end]
+
+        local_slipping_times = velocity_norm[quake_start:quake_end, :] \
+            > THRESHOLD_VELOCITY
+        quake_displacement = displacement[quake_start:quake_end, :, :]
+        slip = np.zeros(num_vertices)
+
+        for i in range(num_vertices):
+            if any(local_slipping_times[:, i]):
+                starts = slip_beginnings(local_slipping_times[:, i])
+                ends = slip_endings(local_slipping_times[:, i])
+                slip[i] = np.linalg.norm(quake_displacement[ends, i, :]
+                                         - quake_displacement[starts, i, :])
+
+        quake['peakSlip'] = max(slip)
+        quake['peakSlipRate'] = max(maximum_velocity[quake_start:quake_end])
+
+        max_rupture_width = 0
+        for time_step in range(local_slipping_times.shape[0]):
+            slipping_time = local_slipping_times[time_step, :]
+            if not any(slipping_time):
+                continue
+
+            slipping_coords = coordinates[slipping_time] \
+                + displacement[quake_start + time_step, slipping_time]
+            if len(slipping_coords) > 1:
+                pair = np.array(max_distance(slipping_coords))
+                max_rupture_width = max(max_rupture_width,
+                                        np.linalg.norm(pair[0] - pair[1]))
+
+        quake['ruptureWidth'] = max_rupture_width
+        quakes[quake_ID] = quake
+
+    # output quake data to csv file
+    csv_columns = quakes[0].keys()
+    csv_file_path = os.path.join(out_path, 'events' + str(body_ID) + '.csv')
+    try:
+        with open(csv_file_path, 'w') as csv_file:
+            writer = csv.DictWriter(csv_file, fieldnames=csv_columns)
+            writer.writeheader()
+            for quake_ID in quakes:
+                writer.writerow(quakes[quake_ID])
+    except IOError:
+        print('I/O error')
+
+h5file.close()

data/config.ini → data/tools/config.ini

 [directories]
-simulation = ~/group/publications/2016-PippingKornhuberRosenauOncken
+simulation = /storage/mi/podlesny/software/dune/build-debug/dune-tectonic/src/strikeslip
 experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
 output     = generated

data/tools/generated/events1.csv

+start,peakSlipRate,end,ruptureWidth,peakSlip
+20.086059570312422,5.7045853762396602e-06,20.086669921874922,0,0.0
+20.090332031249918,5.6219949973608878e-06,20.092163085937418,0,0.0
+20.094604492187415,5.8412336903395239e-06,20.095825195312415,0,0.0
+20.099487304687411,5.3352579051826403e-06,20.100097656249908,0,0.0
+20.143737792968619,5.4285740142335498e-06,20.144958496093619,0,0.0
+20.148620605468615,5.6370561043864526e-06,20.149536132812365,0,0.0
+20.353851318359059,6.0525078256123743e-06,20.354461669921559,0,0.0
+20.358734130859055,5.1806538408711459e-06,20.359344482421555,0,0.0
+20.371856689452795,5.6547173619088325e-06,20.373077392577795,0,0.0
+20.425262451171495,6.4545930360194781e-06,20.428314208983991,0.088385894124722503,1.4298690605481336e-08
+20.479583740233945,5.8310594282664788e-06,20.481414794921445,0,0.0
+20.483245849608942,5.2537138171449348e-06,20.485076904296442,0,6.2000101383992609e-09
+20.531158447265149,5.2316818926703734e-06,20.531768798827649,0,0.0
+20.544586181640138,5.4092184353134343e-06,20.545196533202635,0,0.0
+20.593109130858842,6.190190600109647e-06,20.595550537108842,0,0.0
+20.608673095702581,5.5618674305029011e-06,20.609283447265081,0,0.0
+20.625152587890067,5.4920875934131527e-06,20.625762939452567,0,0.0
+20.630035400390064,5.1740512343599518e-06,20.630645751952564,0,0.0
+20.63369750976506,5.7392806837644885e-06,20.63491821289006,0,0.0

data/tools/support/__init__.py

+ 

data/tools/support/findQuakes.R

-Rcpp::sourceCpp('tools/support/maxVelocity.cpp')
-Rcpp::sourceCpp('tools/support/negativeStarts.cpp')
-Rcpp::sourceCpp('tools/support/positiveStarts.cpp')
-
-findQuakes <- function (criticalVelocity, velocity, indices, dimensions) {
-  maximumVelocities<- maxVelocity(velocity[indices,,dimensions])
-  slippingTimes  <- maximumVelocities > criticalVelocity
-
-  ns <- negativeStarts(slippingTimes) + indices[[1]] - 1
-  ps <- positiveStarts(slippingTimes) + indices[[1]] - 1
-
-  ## NOTE: we drop incomplete quakes here
-  mlen <- min(length(ns), length(ps))
-  ns <- ns[1:mlen]
-  ps <- ps[1:mlen]
-
-  return(data.frame(endingIndex = ns, beginningIndex = ps))
-}

data/tools/support/find_quakes.py

+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+
+
+def find_quakes(threshold_velocity, maximum_velocities):
+    slipping_times = maximum_velocities > threshold_velocity
+
+    quake_starts = slip_beginnings(slipping_times)
+    quake_ends = slip_endings(slipping_times)
+
+    # remove incomplete quakes
+    min_len = min(len(quake_starts), len(quake_ends))
+    quake_ends = quake_ends[0:min_len]
+    quake_starts = quake_starts[0:min_len]
+
+    return [quake_starts, quake_ends]

data/tools/support/find_quakes.py~

+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+
+
+def find_quakes(threshold_velocity, maximum_velocities):
+    slipping_times = maximum_velocities > threshold_velocity
+    print(slipping_times)
+
+    quake_starts = slip_beginnings(slipping_times)
+    quake_ends = slip_endings(slipping_times)
+
+    # remove incomplete quakes
+    min_len = min(len(quake_starts), len(quake_ends))
+    quake_ends = quake_ends[0:min_len]
+    quake_starts = quake_starts[0:min_len]
+
+    return [quake_starts, quake_ends]

data/tools/support/maxVelocity.cpp

-#include <Rcpp.h>
-
-// For std::hypot(x,y)
-// [[Rcpp::plugins(cpp11)]]
-
-// Layout of vectors is such that loops over inner indices are traversed first
-// [[Rcpp::export]]
-Rcpp::NumericVector maxVelocity(Rcpp::NumericVector const &x) {
-  Rcpp::IntegerVector size = x.attr("dim");
-  Rcpp::NumericVector ret(size[0]);
-  switch (size[2]) {
-  case 1:
-    for (size_t ts = 0; ts < size[0]; ++ts)
-      for (size_t coord = 0; coord < size[1]; ++coord)
-        ret[ts] = std::max(ret[ts],
-                           std::abs(x[0 * size[1] * size[0] + coord * size[0] + ts]));
-    break;
-  case 2:
-    for (size_t ts = 0; ts < size[0]; ++ts)
-      for (size_t coord = 0; coord < size[1]; ++coord)
-        ret[ts] = std::max(ret[ts],
-                           std::hypot(x[0 * size[1] * size[0] + coord * size[0] + ts],
-                                      x[1 * size[1] * size[0] + coord * size[0] + ts]));
-    break;
-  default:
-    throw std::range_error("Inadmissible value");
-  }
-  return ret;
-}

data/tools/support/max_distance.py

+from itertools import combinations
+
+
+def square_distance(x, y):
+    return sum([(xi-yi)**2 for xi, yi in zip(x, y)])
+
+
+def max_distance(points):
+    max_square_distance = 0
+    max_pair = (0, 0)
+    for pair in combinations(points, 2):
+        sq_dist = square_distance(*pair)
+
+        if sq_dist > max_square_distance:
+            max_square_distance = sq_dist
+            max_pair = pair
+    return max_pair

data/tools/support/max_distance.py~

+from itertools import combinations
+
+
+def square_distance(x, y):
+    return sum([(xi-yi)**2 for xi, yi in zip(x, y)])
+
+
+def max_distance(points):
+    max_square_distance = 0
+    max_pair = (0, 0)
+    for pair in combinations(points, 2):
+        sq_dist = square_distance(*pair)
+        print(sq_dist)
+        if sq_dist > max_square_distance:
+            max_square_distance = sq_dist
+            max_pair = pair
+    return max_pair

data/tools/support/max_velocity.py

+import numpy as np
+
+def max_velocity(x) :
+    return maximum(norm(x))
\ No newline at end of file

data/tools/support/max_velocity.py~

+import numpy as np
+
+def max_velocity(x) :
+    size = x.shape 
+    res = np.zeros(size[0])
+
+    for time_step in range(size[0]) :
+        for vertex in range(size[1]) :
+            res[time_step] = max(res[time_step], np.linalg.norm(x[time_step, vertex, :]))
+
+    return res
\ No newline at end of file