From 8f824a8a4f28c3607de4a6cf8bea750a9693f008 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Fri, 20 Nov 2020 17:23:12 +0100
Subject: [PATCH] python based data analysis
---
data/tools/2d_event_writer.py | 112 ++++++++++++++++++
data/tools/2d_event_writer.py~ | 112 ++++++++++++++++++
data/{ => tools}/config.ini | 2 +-
data/tools/generate-2d.bash | 0
data/tools/generate-3d.bash | 0
data/tools/generate-others.bash | 0
data/tools/generated/events1.csv | 20 ++++
data/tools/support/__init__.py | 1 +
data/tools/support/__init__.pyc | Bin 0 -> 106 bytes
data/tools/support/findQuakes.R | 18 ---
data/tools/support/find_quakes.py | 16 +++
data/tools/support/find_quakes.pyc | Bin 0 -> 741 bytes
data/tools/support/find_quakes.py~ | 17 +++
data/tools/support/maxVelocity.cpp | 29 -----
data/tools/support/max_distance.py | 17 +++
data/tools/support/max_distance.pyc | Bin 0 -> 884 bytes
data/tools/support/max_distance.py~ | 17 +++
data/tools/support/max_velocity.py | 4 +
data/tools/support/max_velocity.pyc | Bin 0 -> 649 bytes
data/tools/support/max_velocity.py~ | 11 ++
data/tools/support/maximum.py | 11 ++
data/tools/support/maximum.pyc | Bin 0 -> 560 bytes
data/tools/support/maximum.py~ | 11 ++
data/tools/support/negativeStarts.cpp | 13 --
data/tools/support/norm.py | 12 ++
data/tools/support/norm.pyc | Bin 0 -> 616 bytes
data/tools/support/norm.py~ | 11 ++
data/tools/support/positiveStarts.cpp | 14 ---
data/tools/support/slip_beginnings.py | 15 +++
data/tools/support/slip_beginnings.pyc | Bin 0 -> 837 bytes
data/tools/support/slip_beginnings.py~ | 11 ++
data/tools/support/slip_endings.py | 16 +++
data/tools/support/slip_endings.pyc | Bin 0 -> 823 bytes
data/tools/support/slip_endings.py~ | 11 ++
data/tools/support/trapezoidal.cpp | 9 --
data/tools/support/trapezoidal.py | 8 ++
data/tools/support/trapezoidal.py~ | 8 ++
.../{writeContours.R => write_contours.py} | 0
debugging.m | 12 ++
39 files changed, 454 insertions(+), 84 deletions(-)
create mode 100644 data/tools/2d_event_writer.py
create mode 100644 data/tools/2d_event_writer.py~
rename data/{ => tools}/config.ini (56%)
mode change 100755 => 100644 data/tools/generate-2d.bash
mode change 100755 => 100644 data/tools/generate-3d.bash
mode change 100755 => 100644 data/tools/generate-others.bash
create mode 100644 data/tools/generated/events1.csv
create mode 100644 data/tools/support/__init__.py
create mode 100644 data/tools/support/__init__.pyc
delete mode 100644 data/tools/support/findQuakes.R
create mode 100644 data/tools/support/find_quakes.py
create mode 100644 data/tools/support/find_quakes.pyc
create mode 100644 data/tools/support/find_quakes.py~
delete mode 100644 data/tools/support/maxVelocity.cpp
create mode 100644 data/tools/support/max_distance.py
create mode 100644 data/tools/support/max_distance.pyc
create mode 100644 data/tools/support/max_distance.py~
create mode 100644 data/tools/support/max_velocity.py
create mode 100644 data/tools/support/max_velocity.pyc
create mode 100644 data/tools/support/max_velocity.py~
create mode 100644 data/tools/support/maximum.py
create mode 100644 data/tools/support/maximum.pyc
create mode 100644 data/tools/support/maximum.py~
delete mode 100644 data/tools/support/negativeStarts.cpp
create mode 100644 data/tools/support/norm.py
create mode 100644 data/tools/support/norm.pyc
create mode 100644 data/tools/support/norm.py~
delete mode 100644 data/tools/support/positiveStarts.cpp
create mode 100644 data/tools/support/slip_beginnings.py
create mode 100644 data/tools/support/slip_beginnings.pyc
create mode 100644 data/tools/support/slip_beginnings.py~
create mode 100644 data/tools/support/slip_endings.py
create mode 100644 data/tools/support/slip_endings.pyc
create mode 100644 data/tools/support/slip_endings.py~
delete mode 100644 data/tools/support/trapezoidal.cpp
create mode 100644 data/tools/support/trapezoidal.py
create mode 100644 data/tools/support/trapezoidal.py~
rename data/tools/support/{writeContours.R => write_contours.py} (100%)
diff --git a/data/tools/2d_event_writer.py b/data/tools/2d_event_writer.py
new file mode 100644
index 00000000..d36e85be
--- /dev/null
+++ b/data/tools/2d_event_writer.py
@@ -0,0 +1,112 @@
+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()
diff --git a/data/tools/2d_event_writer.py~ b/data/tools/2d_event_writer.py~
new file mode 100644
index 00000000..d36e85be
--- /dev/null
+++ b/data/tools/2d_event_writer.py~
@@ -0,0 +1,112 @@
+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()
diff --git a/data/config.ini b/data/tools/config.ini
similarity index 56%
rename from data/config.ini
rename to data/tools/config.ini
index a3a37a49..d8a6165b 100644
--- a/data/config.ini
+++ b/data/tools/config.ini
@@ -1,4 +1,4 @@
[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
diff --git a/data/tools/generate-2d.bash b/data/tools/generate-2d.bash
old mode 100755
new mode 100644
diff --git a/data/tools/generate-3d.bash b/data/tools/generate-3d.bash
old mode 100755
new mode 100644
diff --git a/data/tools/generate-others.bash b/data/tools/generate-others.bash
old mode 100755
new mode 100644
diff --git a/data/tools/generated/events1.csv b/data/tools/generated/events1.csv
new file mode 100644
index 00000000..57582761
--- /dev/null
+++ b/data/tools/generated/events1.csv
@@ -0,0 +1,20 @@
+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
diff --git a/data/tools/support/__init__.py b/data/tools/support/__init__.py
new file mode 100644
index 00000000..8d1c8b69
--- /dev/null
+++ b/data/tools/support/__init__.py
@@ -0,0 +1 @@
+
diff --git a/data/tools/support/__init__.pyc b/data/tools/support/__init__.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..05977bb84ed00be29153c30f2bd03eb66ffde9aa
GIT binary patch
literal 106
zcmZSn%*%Dk{77sv0~9a<X$K%KW&si@3=F{<AQ3+eAi;n}6bl2zic1R$@{3CJ<Kr{)
dGE3s)^$IFWIDq0dx%nxjIjMFa-Niu6003fJ5+VQq
literal 0
HcmV?d00001
diff --git a/data/tools/support/findQuakes.R b/data/tools/support/findQuakes.R
deleted file mode 100644
index 698b6498..00000000
--- a/data/tools/support/findQuakes.R
+++ /dev/null
@@ -1,18 +0,0 @@
-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))
-}
diff --git a/data/tools/support/find_quakes.py b/data/tools/support/find_quakes.py
new file mode 100644
index 00000000..aa4653fd
--- /dev/null
+++ b/data/tools/support/find_quakes.py
@@ -0,0 +1,16 @@
+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]
diff --git a/data/tools/support/find_quakes.pyc b/data/tools/support/find_quakes.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..3ddaf405b860cdef5b12af1266eea56bede793b2
GIT binary patch
literal 741
zcmcIh%}(P$4ED^VfS_FvLgK)QTQ4v#uv#H;XQgthgrblc+7a^;I?h6+_Qam{g?I~I
zh6eyUr3>OtlyN@$bNtz3^)lH0_w;;K!}}xYdqyXpkTLWM@Bl%;Sb$((f|n2!j1~C`
z{sN2(Fa=x=d-R8UrDrGz;VxJ4t|_nW#Kq|1B(I46ny}Guv}Vl>It86`Cb^r1fCa#9
z0ZR!;u7pPcO9e>r6)pgYQm{}=F(?>3N?7(7Q|pT!94`ybSUX)hkrPQBv+#u@14Gsx
zF&cHvD?-6VDhl9j9CCZO&t^Kad78W_|5%@Dhx0BRss|U^unsz#@dst|{0+@mq6-w_
z$?w|g)|MHo8D`x6AE3L)xJfRlmaIq+I+K4B($f^I4o)}8c$?!~=X8Vj)y!(s#+5li
zTVsl@)~3QrV@f{j+%`>`p}uj^l$~}yZRWT`H9v)Tb(tFm(U*PES4UzX2J)iYDIPMr
ss&o2fqurI^!PkWM&31P=;vqlA`M#O9NatrrruFt0FER5@_r#uh1Fgxl=l}o!
literal 0
HcmV?d00001
diff --git a/data/tools/support/find_quakes.py~ b/data/tools/support/find_quakes.py~
new file mode 100644
index 00000000..e5aef040
--- /dev/null
+++ b/data/tools/support/find_quakes.py~
@@ -0,0 +1,17 @@
+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]
diff --git a/data/tools/support/maxVelocity.cpp b/data/tools/support/maxVelocity.cpp
deleted file mode 100644
index 5294e563..00000000
--- a/data/tools/support/maxVelocity.cpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#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;
-}
diff --git a/data/tools/support/max_distance.py b/data/tools/support/max_distance.py
new file mode 100644
index 00000000..95467bad
--- /dev/null
+++ b/data/tools/support/max_distance.py
@@ -0,0 +1,17 @@
+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
diff --git a/data/tools/support/max_distance.pyc b/data/tools/support/max_distance.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..7fc47e8672dbe1ee8d9da3ab51be2c4d3346a3fc
GIT binary patch
literal 884
zcmcIi!EV$r5PeQUp=eiBS_xICXSA1+3rAET!FdJSL#2v>LL9(i6KCUz*sZds_SjG0
zKllrNh7W)@-ik%s>)3Ce$If`>+3BzG!|y+SUX-|P3GGXc@RVTL4HiIcaBg61f@|PB
zflVVxVKbag?{Ta9%{Db;G9r||{^*)QoNq#@F*6R`_lhGtCkXJFXfObS!Bv9S7f<jm
z!G=A>kYGr0$m1vP6D%jVGCWBm&75MSl4`WRj`q^EQ$2+tfmOt{=+>2sX6>e&AVPRU
zmVfbemDetBy<JqHS?9rjl5a&<<+g9Oa4c0RzHw!43n_B(ei8D}x2^9auZz_O>w*+b
zSv_ml@`yphvS)z5ZcQj?5-_Far4ne;QHHl9EO28~PnNwpQB|l1b&h%*s}d?viW^{)
zVPH4-nxQ*V&Of=tKbJ?^&nVhy)L`Zq8atu3Q}spN=E~;xb*c%;9-;NFksuRNt?R!!
z<iV}KEnFvBRU{fBEMv8IEv~+*e+Acs!T|whkKbc6N=CF}Guys`eXcmEx_D)~g6CuU
bcezG(H?Khxd{KMbFRBy0`k<;Er2E+)26@Ja
literal 0
HcmV?d00001
diff --git a/data/tools/support/max_distance.py~ b/data/tools/support/max_distance.py~
new file mode 100644
index 00000000..2bb87571
--- /dev/null
+++ b/data/tools/support/max_distance.py~
@@ -0,0 +1,17 @@
+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
diff --git a/data/tools/support/max_velocity.py b/data/tools/support/max_velocity.py
new file mode 100644
index 00000000..d3ee84eb
--- /dev/null
+++ b/data/tools/support/max_velocity.py
@@ -0,0 +1,4 @@
+import numpy as np
+
+def max_velocity(x) :
+ return maximum(norm(x))
\ No newline at end of file
diff --git a/data/tools/support/max_velocity.pyc b/data/tools/support/max_velocity.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..1e0953b749ce794147ae298b2c25c36082d51a0b
GIT binary patch
literal 649
zcmcIhK~4fO6n*W03elAd!wuMAFF=S1tclSDi3uT-nF1v<(_#A2!N5x0c`~oyF+70&
z7L8uOX8!N*`>(IBo#^ANJ$!$gCiotb=PM4^;gs?xkO6lBsl+oE(j&Lsh`Z*K$2bvO
zHaXfHLk{=KX$P!;y#NWUBzC-Cpz|cKm#7UW`U8!o6YL|@;~!WC=4jf-z5!}r#{v;L
zUaE=Jkr(Int`YnJFA{TZv-lGe1Gm2(VY2RF!pno%v0+y*-fZ`n6M4&S8JC*;m5VX}
zkGirhR8?FoG}{F7xTYhXGgHLbg6%NdR(aP8uqB`lm@|(Dex)6?e9z>1=17-e_r0!2
z*ImXDlDp<gEp4t<Zd7U0OuJ&MoSoA{TxpeVibLv=PRJG}QE5!EB5O0J+@>sTMT))6
z?saApL)*)8OQ$UQ7xVs|2O(PGTwaI^F_7Wsg-aHjyxiXSJKT8Eue*P=c)|zubDM55
MJ@6+x|NN)&3q=!?E&u=k
literal 0
HcmV?d00001
diff --git a/data/tools/support/max_velocity.py~ b/data/tools/support/max_velocity.py~
new file mode 100644
index 00000000..b0ef5d74
--- /dev/null
+++ b/data/tools/support/max_velocity.py~
@@ -0,0 +1,11 @@
+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
diff --git a/data/tools/support/maximum.py b/data/tools/support/maximum.py
new file mode 100644
index 00000000..e7666fbe
--- /dev/null
+++ b/data/tools/support/maximum.py
@@ -0,0 +1,11 @@
+import numpy as np
+
+
+def maximum(x):
+ size = x.shape
+ res = np.zeros(size[0])
+
+ for time_step in range(size[0]):
+ res[time_step] = max(res[time_step], max(x[time_step, :]))
+
+ return res
diff --git a/data/tools/support/maximum.pyc b/data/tools/support/maximum.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..30e653a21ca5c22d270af2091ca051cb8f599f8f
GIT binary patch
literal 560
zcmcIgK}rKL6n#l5wnbMiq@G|lHHh0Hg1QpAXcvM|WSB|ql*uGcGT4@ORlSPG@iHDj
zf6{_pz%c*K%YXTKkBG11?%BuZB*T7j`d?GI9zjkYLkT=G<Q#8Ac<781L2B_$YnU;j
zEvh!v4V9Y_tbsYOVekxdjy26Q^tKEej>R3k_|v+bTw){OTUgUcfnH#Ag^2#ZJI9pA
z;VtR$$l#@q2cqL7<Q*0SNDc^m05?mm@<FL=Y^+pmTx4}xJu7ktbh=1FYzhc0$dcxY
z)lRFG^PQj-T0J?hY(g7EaQDRI%xEQbA+5<v<*KE0X6j$lTFJbrcF^}K^QI~?nWsLL
z-k8!!*H~+6FUh>njXt+ae@OT3w=qLESg{T};iv3^$A^x|t47=9D1;q_Xupa6Csvz0
R2HNCJsjh;l3(t8ZegNHVe7yhw
literal 0
HcmV?d00001
diff --git a/data/tools/support/maximum.py~ b/data/tools/support/maximum.py~
new file mode 100644
index 00000000..3664f3bd
--- /dev/null
+++ b/data/tools/support/maximum.py~
@@ -0,0 +1,11 @@
+import numpy as np
+
+def maximum(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], max(x[time_step, vertex]))
+
+ return res
\ No newline at end of file
diff --git a/data/tools/support/negativeStarts.cpp b/data/tools/support/negativeStarts.cpp
deleted file mode 100644
index 940765c2..00000000
--- a/data/tools/support/negativeStarts.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-#include <Rcpp.h>
-
-// [[Rcpp::export]]
-Rcpp::NumericVector negativeStarts(Rcpp::NumericVector const &x) {
- std::vector<size_t> starts(0);
- bool prev = false;
- for (size_t i = 0; i < x.size(); ++i) {
- if (prev && !x[i])
- starts.push_back(i+1); // R indices start at 1
- prev = x[i];
- }
- return Rcpp::NumericVector(starts.begin(), starts.end());
-}
diff --git a/data/tools/support/norm.py b/data/tools/support/norm.py
new file mode 100644
index 00000000..c0276d80
--- /dev/null
+++ b/data/tools/support/norm.py
@@ -0,0 +1,12 @@
+import numpy as np
+
+
+def norm(x):
+ size = x.shape
+ res = np.zeros((size[0], size[1]))
+
+ for time_step in range(size[0]):
+ for vertex in range(size[1]):
+ res[time_step, vertex] = np.linalg.norm(x[time_step, vertex, :])
+
+ return res
diff --git a/data/tools/support/norm.pyc b/data/tools/support/norm.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..ba10f9d1d00a542cac6c76c38fa465f6a1633e7e
GIT binary patch
literal 616
zcmb_Z!AiqG5PiEzrImUvJ$TPW_oDPtL{QLE5j_-yA~D&vDM>cXZZs{;Nj>`)evUum
z2k7k9;tv=yZ-+N;cHSi3*TMeD$LA!%@5AXlr_g<Zm|lhyxMzqtUMb<0xeGk1Cg1do
zB1UP8qD^r@q2CEMzyfF(48wwBL;VbWdvNrp94`zFNBs!n?TSQQPHxb6sIPW3i5-24
z#&=iwShVcr8*1~NKj4jH>N(+tCb<Gl7%oQ`F9#SCx%E3r$QTmqH`pUE*0r98g*0~S
zyl})?mP)&<4D*>ZR+%Pwn9huERjSMa+enK6wc~WMmb3>@N^Lr}mSpll8(FySCs`U<
z2gJjHzNSpfRVGE2h(g7w)cHziH8sy+DMei6TR1W@GAd6Z5r-xeMyXT_T@{5Yjj%?J
zi`A|5VvXr80a%9}@*eB40eAn%Es|Fm39VnJ{&VvGc1)Y53^Ns1sXVh2wf(eR?*}iy
Bg<${y
literal 0
HcmV?d00001
diff --git a/data/tools/support/norm.py~ b/data/tools/support/norm.py~
new file mode 100644
index 00000000..3e26f65d
--- /dev/null
+++ b/data/tools/support/norm.py~
@@ -0,0 +1,11 @@
+import numpy as np
+
+def norm(x) :
+ size = x.shape
+ res = np.zeros(size[0], size[1])
+
+ for time_step in range(size[0]) :
+ for vertex in range(size[1]) :
+ res[time_step, vertex] = np.linalg.norm(x[time_step, vertex, :])
+
+ return res
diff --git a/data/tools/support/positiveStarts.cpp b/data/tools/support/positiveStarts.cpp
deleted file mode 100644
index 3e9498c0..00000000
--- a/data/tools/support/positiveStarts.cpp
+++ /dev/null
@@ -1,14 +0,0 @@
-#include <Rcpp.h>
-
-// [[Rcpp::export]]
-Rcpp::NumericVector positiveStarts(Rcpp::NumericVector const &x) {
- std::vector<size_t> starts(0);
- bool prev = false;
- for (size_t i = 0; i < x.size(); ++i) {
- if (!prev && x[i])
- starts.push_back(i+1); // R indices start at 1
- prev = x[i];
- }
- return Rcpp::NumericVector(starts.begin(), starts.end());
-}
-
diff --git a/data/tools/support/slip_beginnings.py b/data/tools/support/slip_beginnings.py
new file mode 100644
index 00000000..255adf44
--- /dev/null
+++ b/data/tools/support/slip_beginnings.py
@@ -0,0 +1,15 @@
+import array as ar
+
+
+def slip_beginnings(x):
+ # returns indicies i for which x[i-1]=0 and x[i]>0
+ starts = ar.array('i', (0 for i in range(len(x))))
+ prev = False
+
+ length = 0
+ for i in range(len(x)):
+ if (not prev and x[i]):
+ starts[length] = i
+ length += 1
+ prev = x[i]
+ return starts[:length]
diff --git a/data/tools/support/slip_beginnings.pyc b/data/tools/support/slip_beginnings.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..6ab6756aec374c95a324af74844469346be2e0cf
GIT binary patch
literal 837
zcmcgqO>Yx15PjYa5%hzql{n<aF&FCI$|0(dis%L9kVqgQMQF*It*ec9YddWsB`3J`
zfAr_@XZQi|#!ET$4wh#;JDxY5H?#TA{hg!xAG10Btc>RyN*GX;+zK`DNukW}okk{}
z=`iQ5{lzt&E1FLzeabl{+)+i~ikJsj%n);IfNx6nTRAdp6=Dxtrp&R;v3ifkakJOh
z_OQvY=a&NC^e_;MK7~z=5vA<mlxWA>E`SNIFYvL8csU(PapNdu+ElDi6bT|PsEC{G
zN;TyBL0vPYP?Qi$!FVV-eX~B`_+*kCe^r>$)R)b|=!Mm7Q`RQ9wGPeKcw2a*%cbko
zt7zu2adxiDA{II}O&xStwr%614z+DRe=(JH&blfbx9jkXc}JCTX4U#L>2Hwt`(0U9
z7$tJ4s=`MhT6kZqljRFn8Fn40jZ2<)MIB6s5bP?Zl#$55IRLFUH_1mp6)*o2%Xq-2
zACZ@8!Vwj+q3UOSHBg6YNA0WUD*bvBp*4|aOm?Zd0nzBJe>IQzlrM}HO}VVinRGCX
HV=w;=3RJnG
literal 0
HcmV?d00001
diff --git a/data/tools/support/slip_beginnings.py~ b/data/tools/support/slip_beginnings.py~
new file mode 100644
index 00000000..6c54c9af
--- /dev/null
+++ b/data/tools/support/slip_beginnings.py~
@@ -0,0 +1,11 @@
+import numpy as np
+
+def slip_beginnings(x) :
+ # returns indicies i for which x[i-1]=0 and x[i]>0
+ starts = []
+ prev = False
+ for i in range(len(x)) :
+ if (not prev and x[i]) :
+ starts.append(i)
+ prev = x[i]
+ return starts
diff --git a/data/tools/support/slip_endings.py b/data/tools/support/slip_endings.py
new file mode 100644
index 00000000..9fdcaf88
--- /dev/null
+++ b/data/tools/support/slip_endings.py
@@ -0,0 +1,16 @@
+import array as ar
+
+
+def slip_endings(x):
+ # returns indicies i for which x[i-1]>0 and x[i]=0
+ ends = ar.array('i', (0 for i in range(len(x))))
+ prev = False
+
+ length = 0
+ for i in range(len(x)):
+ if (prev and not x[i]):
+ ends[length] = i
+ length += 1
+ prev = x[i]
+
+ return ends[:length]
diff --git a/data/tools/support/slip_endings.pyc b/data/tools/support/slip_endings.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..89caa52ad2f9e896a68667e641bce5836ddceaeb
GIT binary patch
literal 823
zcmcgqO>Yx15PkNBh^8qLl{n<qV=mO(;1ETXKyX13N~BgoiqI-+(nTBZ#&&2TC8u)Z
zfADkoGyDMh#tRkV4wh#;JD!jC#+&^(>0NyLK40M1(mXFHVL(-ID^$RHg+k*C4K1G8
zFyn3clWVe2G#^m9l$Vt7l_~<CnDaX=YQzkC;8w|gC5OgAA$D+J$qa`K+qd{TZvFy?
z4)z*;BiZ{727NK<Q`l!1QHl<pGu!dHjbOp`8(g&!XS1Pzn?!+GCa^+LB#1nrV&42%
z(2(z?TC=23l#r5w$xw9qCO+W&*)%!sRd`3!T-PgWR;6j`qO!s5OsGG`o7`JdY+S3J
zM!Sf$D;K87V{T$ySHXl$)6_niP?gOG>x$AX!+Em{M=Uy8T4%S7AIorqH1D>UN}-dO
zB}loC!ZP<h-zCfEZfQAkpt3G`UguS?ZHDA-IeCEXLNXd}KSyc8@MXOI9|U8c?fxV!
zfWj#i^iXwmR}Iu-)l(<xkxu{HROn1)2-9P#_El)K)<1cK_>kpCtGd`!b}VB{Z_v-~
E0A$^>+yDRo
literal 0
HcmV?d00001
diff --git a/data/tools/support/slip_endings.py~ b/data/tools/support/slip_endings.py~
new file mode 100644
index 00000000..f05b69c9
--- /dev/null
+++ b/data/tools/support/slip_endings.py~
@@ -0,0 +1,11 @@
+import numpy as np
+
+def slip_endings(x) :
+ # returns indicies i for which x[i-1]>0 and x[i]=0
+ starts = np.array([])
+ prev = False
+ for i in range(len(x)) :
+ if (prev and not x[i]) :
+ starts.append(i)
+ prev = x[i]
+ return starts
\ No newline at end of file
diff --git a/data/tools/support/trapezoidal.cpp b/data/tools/support/trapezoidal.cpp
deleted file mode 100644
index 3b4833d6..00000000
--- a/data/tools/support/trapezoidal.cpp
+++ /dev/null
@@ -1,9 +0,0 @@
-#include <Rcpp.h>
-
-// [[Rcpp::export]]
-double trapezoidal(Rcpp::NumericVector const &x, Rcpp::NumericVector const &y) {
- double ret = 0;
- for (size_t i = 1; i < x.size(); ++i)
- ret += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2;
- return ret;
-}
diff --git a/data/tools/support/trapezoidal.py b/data/tools/support/trapezoidal.py
new file mode 100644
index 00000000..9574b125
--- /dev/null
+++ b/data/tools/support/trapezoidal.py
@@ -0,0 +1,8 @@
+# x: vector of domain values
+# y: vector of function values f(x)
+def trapezoidal(x, y) :
+ # returns integral of (x, y = f(x)) calculated with trapezoidal rule
+ res = 0.0
+ for i in range(1, len(x)) :
+ res += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2.0
+ return res
\ No newline at end of file
diff --git a/data/tools/support/trapezoidal.py~ b/data/tools/support/trapezoidal.py~
new file mode 100644
index 00000000..8d462782
--- /dev/null
+++ b/data/tools/support/trapezoidal.py~
@@ -0,0 +1,8 @@
+# x: vector of domain values
+# y: vector of function values f(x)
+def trapezoidal(x, y) :
+ # returns integral of (x, y = f(x)) calculated with trapezoidal rule
+ res = 0.0
+ for i in range(1, len(x)) :
+ res += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2.0
+ return res
\ No newline at end of file
diff --git a/data/tools/support/writeContours.R b/data/tools/support/write_contours.py
similarity index 100%
rename from data/tools/support/writeContours.R
rename to data/tools/support/write_contours.py
diff --git a/debugging.m b/debugging.m
index 259ae8db..0146058d 100644
--- a/debugging.m
+++ b/debugging.m
@@ -1,3 +1,15 @@
+% householder reflection
+
+v0 = [1; 0];
+v0 = v0./norm(v0);
+
+v1 = [1; -1];
+v1 = v1./norm(v1);
+
+v = 1/2.*(v1-v0);
+vnorm = norm(v);
+
+H = eye(2) - 2/(vnorm^2).*(v*v')
A = [...
];
--
GitLab