more usable framework

6 files changed, 69 insertions, 40 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..743e88b
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+__pycache__/
+env/
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..5a2b0c0
--- /dev/null
+++ b/README.md
@@ -0,0 +1,5 @@
+Data units are ms and ft. Data derived from Kinovea.
+
+Create a venv and then `pip install -r requirements.txt` and then
+`python run.py` after data is entered in the data folder corresponding
+to `data/robot[robot number]_trial[trial number].csv`
diff --git a/data/robot6_trial2_edited.csv b/data/robot6_trial2_edited.csv
deleted file mode 100644
index 7935553..0000000
--- a/data/robot6_trial2_edited.csv
+++ /dev/null
@@ -1,40 +0,0 @@
-,Time (ms),Trajectory 1
-0,0,-0.979554891586304
-1,4,-0.835712134838104
-2,8,-0.699505805969238
-3,12,-0.548177599906921
-4,16,-0.321406573057175
-5,20,-0.174405112862587
-6,25,-0.0335861444473267
-7,29,0.0799026042222977
-8,33,0.256950944662094
-9,37,0.400715202093124
-10,41,0.544477641582489
-11,45,0.680721402168274
-12,49,0.824513614177704
-13,54,0.968290090560913
-14,58,1.10454130172729
-15,61,1.24079287052155
-16,65,1.39210069179535
-17,70,1.52830684185028
-18,74,1.66456472873688
-19,78,1.81589031219482
-20,82,1.95213425159454
-21,86,2.09591507911682
-22,90,2.23969602584839
-23,94,2.37587356567383
-24,99,2.50452637672424
-25,103,2.64827156066895
-26,106,2.78445959091187
-27,111,2.95850706100464
-28,115,3.0796012878418
-29,119,3.2309718132019
-30,123,3.35961222648621
-31,127,3.49582982063293
-32,131,3.63205075263977
-33,135,3.77579402923584
-34,139,3.91195726394653
-35,144,4.04812145233154
-36,148,4.18429899215698
-37,151,4.32052373886108
-38,156,4.44540596008301
diff --git a/regression.py b/regression.py
new file mode 100644
index 0000000..b6ac7db
--- /dev/null
+++ b/regression.py
@@ -0,0 +1,33 @@
+import numpy as np
+
+# data = [(0,0), (1,1), (2, 4), (3, 9), (4, 16) ] ; t = 3
+def velvar(data, t):
+    n = len(data) # number of entries
+    X = np.array([(datum[0]**2, datum[0], 1) for datum in data])
+    Xt = np.transpose(X) # X^T
+    Y = [datum[1] for datum in data]
+    M = np.matmul( np.linalg.inv(np.matmul(Xt, X)), Xt ) # solution = My
+    est = np.matmul(M, Y);
+    P = np.matmul(X, M) # \hat y = Py
+    R = np.subtract(np.identity(n), P) # residual = Ry
+    bias = sum([sum([entry**2 for entry in row]) for row in R])
+    Ry = np.matmul(R, Y);
+    S2 = np.dot(Ry, Ry); # S^2 estimator
+    variance = S2/bias
+    v = 2*M[0]*t + M[1] # V = v\cdot y
+    velocity = 2*est[0]*t+est[1]
+    variancevel = np.dot(v,v)*variance
+    return (np.abs(velocity), variancevel)
+
+def combine_estimates(velvars):
+    avg = np.average(list(a[0] for a in velvars))
+    var = (np.sum(a[1] + a[0]**2 for a in velvars) - avg**2)/len(velvars)
+    return((avg, var))
+
+if __name__ == '__main__':
+    tot = (velvar([(0,0), (1,1), (2, 4), (3, 9), (4, 17)], 3),
+           velvar([(0, 0), (1,2), (2, 8), (3, 18), (4, 32)], 3))
+    n = 5
+    k = 2
+    print(sum([datum[1]**2+n*datum[0]**2 for datum in tot])
+        - n*k*sum([datum[0]/k for datum in tot])**2)
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..c0358f4
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,5 @@
+numpy==1.21.3
+pandas==1.3.4
+python-dateutil==2.8.2
+pytz==2021.3
+six==1.16.0
diff --git a/run.py b/run.py
new file mode 100644
index 0000000..7ff0b74
--- /dev/null
+++ b/run.py
@@ -0,0 +1,24 @@
+#!/usr/bin/python3
+import os
+import re
+import pandas as pd
+import numpy as np
+from regression import velvar, combine_estimates
+
+robot_velvars = {}
+for entry in os.scandir('data'):
+    if not (entry.is_file() and entry.name.endswith('.csv')):
+        continue
+    match = re.match(r'^robot(\d+)_trial(\d+).csv$', entry.name)
+    if match == None:
+        continue
+    robot = int(match.group(1))
+    if not robot in robot_velvars:
+        robot_velvars[robot] = []
+    trial = int(match.group(2))
+    data = pd.read_csv(entry.path).values
+    robot_velvars[robot].append(velvar(data, data[0][0]))
+
+for robot in robot_velvars:
+    est = [x*1000/3.28 for x in combine_estimates(robot_velvars[robot])]
+    print(f"robot {robot}'s kick velocity is {est[0]} ± {np.sqrt(est[1])*2} m/s (p > .95)")