working force plot

rowers/interactiveplots.py

@@ -64,6 +64,7 @@ activate(settings.TIME_ZONE)
 thetimezone = get_current_timezone()
 
 from scipy.stats import linregress,percentileofscore
+from scipy.spatial import ConvexHull,Delaunay
 from scipy import optimize
 from scipy.signal import savgol_filter
 from scipy.interpolate import griddata
@@ -407,172 +408,207 @@ def interactive_forcecurve(theworkouts,workstrokesonly=False):
     # quick linear regression
     # peakforce = slope*peakforceangle + intercept
     slope, intercept, r,p,stderr = linregress(rowdata['peakforceangle'],rowdata['peakforce'])
-    theta = np.arctan(slope)
 
-    a = rowdata['peakforceangle']-rowdata['peakforceangle'].mean()
+    covariancematrix = np.cov(rowdata['peakforceangle'],y=rowdata['peakforce'])
-    F = rowdata['peakforce']-rowdata['peakforce'].mean()
+    eig_vals, eig_vecs = np.linalg.eig(covariancematrix)
 
-    R = np.array([[np.cos(theta),np.sin(theta)],
+    a = rowdata['peakforceangle']-rowdata['peakforceangle'].median()
-                  [-np.sin(theta),np.cos(theta)]])
+    F = rowdata['peakforce']-rowdata['peakforce'].median()
-    Rinv = np.array([[np.cos(theta),np.sin(theta)],
+
-                     [-np.sin(theta),np.cos(theta)]])
+    Rinv = eig_vecs
+    R = np.linalg.inv(Rinv)
 
     x = R[0,0]*a+R[0,1]*F
     y = R[1,0]*a+R[1,1]*F
 
 
-    x05 = x.quantile(q=0.05)
+    x05 = x.quantile(q=0.01)
-    x25 = x.quantile(q=0.25)
+    x25 = x.quantile(q=0.15)
-    x75 = x.quantile(q=0.75)
+    x75 = x.quantile(q=0.85)
-    x95 = x.quantile(q=0.95)
+    x95 = x.quantile(q=0.99)
+    
+    y05 = y.quantile(q=0.01)
+    y25 = y.quantile(q=0.15)
+    y75 = y.quantile(q=0.85)
+    y95 = y.quantile(q=0.99)
+
+    a25 = Rinv[0,0]*x25 + rowdata['peakforceangle'].median()
+    F25 = Rinv[1,0]*x25 + rowdata['peakforce'].median()
+
+    a25b = Rinv[0,1]*y25 + rowdata['peakforceangle'].median()
+    F25b = Rinv[1,1]*y25 + rowdata['peakforce'].median()
+
+    a75 = Rinv[0,0]*x75 + rowdata['peakforceangle'].median()
+    F75 = Rinv[1,0]*x75 + rowdata['peakforce'].median()
+
+    a75b = Rinv[0,1]*y75 + rowdata['peakforceangle'].median()
+    F75b = Rinv[1,1]*y75 + rowdata['peakforce'].median()
+
+    
+    a05 = Rinv[0,0]*x05 + rowdata['peakforceangle'].median()
+    F05 = Rinv[1,0]*x05 + rowdata['peakforce'].median()
+
+    a05b = Rinv[0,1]*y05 + rowdata['peakforceangle'].median()
+    F05b = Rinv[1,1]*y05 + rowdata['peakforce'].median()
+
+    a95 = Rinv[0,0]*x95 + rowdata['peakforceangle'].median()
+    F95 = Rinv[1,0]*x95 + rowdata['peakforce'].median()
+
+    a95b = Rinv[0,1]*y95 + rowdata['peakforceangle'].median()
+    F95b = Rinv[1,1]*y95 + rowdata['peakforce'].median()
+
     
-    y05 = y.quantile(q=0.05)
-    y25 = y.quantile(q=0.25)
-    y75 = y.quantile(q=0.75)
-    y95 = y.quantile(q=0.95)
 
     try:
-        catchav = rowdata['catch'].mean()
+        catchav = rowdata['catch'].median()
         catch25 = rowdata['catch'].quantile(q=0.25)
         catch75 = rowdata['catch'].quantile(q=0.75)
         catch05 = rowdata['catch'].quantile(q=0.05)
-        catch99 = rowdata['catch'].quantile(q=0.95)
+        catch95 = rowdata['catch'].quantile(q=0.95)
     except KeyError:
         catchav = 0
         catch25 = 0
         catch75 = 0
         catch05 = 0
-        catch99 = 0
+        catch95 = 0
 
     try:
-        finishav = rowdata['finish'].mean()
+        finishav = rowdata['finish'].median()
         finish25 = rowdata['finish'].quantile(q=0.25)
         finish75 = rowdata['finish'].quantile(q=0.75)
         finish05 = rowdata['finish'].quantile(q=0.05)
-        finish99 = rowdata['finish'].quantile(q=0.95)
+        finish95 = rowdata['finish'].quantile(q=0.95)
     except KeyError:
         finishav = 0
         finish25 = 0
         finish75 = 0
         finish05 = 0
-        finish99 = 0
+        finish95 = 0
         
     try:
-        washav = (rowdata['finish']-rowdata['wash']).mean()
+        washav = (rowdata['finish']-rowdata['wash']).median()
         wash25 = (rowdata['finish']-rowdata['wash']).quantile(q=0.25)
         wash75 = (rowdata['finish']-rowdata['wash']).quantile(q=0.75)
         wash05 = (rowdata['finish']-rowdata['wash']).quantile(q=0.05)
-        wash99 = (rowdata['finish']-rowdata['wash']).quantile(q=0.95)
+        wash95 = (rowdata['finish']-rowdata['wash']).quantile(q=0.95)
     except KeyError:
         washav = 0
         wash25 = 0
         wash75 = 0
         wash05 = 0
-        wash99 = 0
+        wash95 = 0
 
     try:
-        slipav = (rowdata['slip']+rowdata['catch']).mean()
+        slipav = (rowdata['slip']+rowdata['catch']).median()
         slip25 = (rowdata['slip']+rowdata['catch']).quantile(q=0.25)
         slip75 = (rowdata['slip']+rowdata['catch']).quantile(q=0.75)
         slip05 = (rowdata['slip']+rowdata['catch']).quantile(q=0.05)
-        slip99 = (rowdata['slip']+rowdata['catch']).quantile(q=0.95)
+        slip95 = (rowdata['slip']+rowdata['catch']).quantile(q=0.95)
     except KeyError:
         slipav = 0
         slip25 = 0
         slip75 = 0
         slip05 = 0
-        slip99 = 0
+        slip95 = 0
         
     try:
-        peakforceav = rowdata['peakforce'].mean()
+        peakforceav = rowdata['peakforce'].median()
         peakforce25 = rowdata['peakforce'].quantile(q=0.25)
         peakforce75 = rowdata['peakforce'].quantile(q=0.75)
         peakforce05 = rowdata['peakforce'].quantile(q=0.05)
-        peakforce99 = rowdata['peakforce'].quantile(q=0.95)
+        peakforce95 = rowdata['peakforce'].quantile(q=0.95)
     except KeyError:
         peakforceav = 0
         peakforce25 = 0
         peakforce75 = 0
         peakforce05 = 0
-        peakforce99 = 0
+        peakforce95 = 0
 
     
     try:
-        averageforceav = rowdata['averageforce'].mean()
+        averageforceav = rowdata['averageforce'].median()
     except KeyError:
         averageforceav = 0
         
     try:
-        peakforceangleav = rowdata['peakforceangle'].mean()
+        peakforceangleav = rowdata['peakforceangle'].median()
         peakforceangle05 = rowdata['peakforceangle'].quantile(q=0.05)
         peakforceangle25 = rowdata['peakforceangle'].quantile(q=0.25)
         peakforceangle75 = rowdata['peakforceangle'].quantile(q=0.75)        
-        peakforceangle99 = rowdata['peakforceangle'].quantile(q=0.95)        
+        peakforceangle95 = rowdata['peakforceangle'].quantile(q=0.95)        
     except KeyError:
         peakforceangleav = 0
         peakforceangle25 = 0
         peakforceangle75 = 0
         peakforceangle05 = 0
-        peakforceangle99 = 0
+        peakforceangle95 = 0
 
 
+    #thresholdforce /= 4.45 # N to lbs
+    thresholdforce = 100 if 'x' in boattype else 200
+    points2575 = [
+        (catch25,0),  #0
+        (slip25,thresholdforce), #1
+        (a75,F75),#4
+        (a25b,F25b), #9
+        (a25,F25), #2
+        (wash75,thresholdforce), #5
+        (finish75,0), #6
+        (finish25,0), #7
+        (wash25,thresholdforce), #8
+        (a75b,F75b), #3
+        (slip75,thresholdforce), #10
+        (catch75,0), #11
+        ]
+
+    points0595 = [
+        (catch05,0),  #0
+        (slip05,thresholdforce), #1
+        (a95,F95),#4
+        (a05b,F05b), #9
+        (a05,F05), #2
+        (wash95,thresholdforce), #5
+        (finish95,0), #6
+        (finish05,0), #7
+        (wash05,thresholdforce), #8
+        (a95b,F95b), #3
+        (slip95,thresholdforce), #10
+        (catch95,0), #11
+        ]
+
+
+    hull2575 = ConvexHull(points2575)
+
+    
+    angles2575 = []
+    forces2575 = []
+
+    for x,y in points2575:
+        angles2575.append(x)
+        forces2575.append(y)
+
+    
+    angles0595 = []
+    forces0595 = []
+
+    for x,y in points0595:
+        angles0595.append(x)
+        forces0595.append(y)
+
+
+
+    
+        
     x = [catchav,
          slipav,
          peakforceangleav,
          washav,
          finishav]
 
-    thresholdforce = 100 if 'x' in boattype else 200
-    #thresholdforce /= 4.45 # N to lbs
     y = [0,thresholdforce,
          peakforceav,
          thresholdforce,0]
 
-    x2575 = [catch25,
-             slip25,
-             peakforceangle25,
-#             peakforceangleav,
-             peakforceangle75,
-           wash75,
-           finish75,
-           finish25,
-           wash25,
-           peakforceangleav,
-           slip75]
-
-    y2575 = [0,
-           thresholdforce,
-             peakforce25,
-#             peakforceav,
-             peakforce75,
-           thresholdforce,
-           0,
-           0,
-           thresholdforce,
-           peakforce25,
-           0]
-
-    
-
-    x0599 = [catch05,
-             slip05,
-             peakforceangleav,
-             wash99,
-             finish99,
-             finish05,
-             wash05,
-             peakforceangleav,
-             slip99]
-
-    y0599 = [0,
-           thresholdforce,
-           peakforce99,
-           thresholdforce,
-           0,
-           0,
-           thresholdforce,
-           peakforce05,
-           0]
 
     
     source = ColumnDataSource(
@@ -581,21 +617,13 @@ def interactive_forcecurve(theworkouts,workstrokesonly=False):
             y = y,
         ))
 
-    sourcerange = ColumnDataSource(
-        data = dict(
-            x2575=x2575,
-            y2575=y2575,
-            x0599=x0599,
-            y0599=y0599,
-        ))
-
     sourceslipwash = ColumnDataSource(
         data = dict(
             xslip = [slipav,washav],
             yslip = [thresholdforce,thresholdforce]
             )
         )
-
+    
     sourcetrend = ColumnDataSource(
         data = dict(
             x = [peakforceangle25,peakforceangle75],
@@ -605,8 +633,8 @@ def interactive_forcecurve(theworkouts,workstrokesonly=False):
 
     sourcefit = ColumnDataSource(
         data = dict(
-            x = rowdata['peakforceangle'],
+            x = np.array([peakforceangle25,peakforceangle75]),
-            y = slope*rowdata['peakforceangle']+intercept
+            y = slope*np.array([peakforceangle25,peakforceangle75])+intercept
             )
         )
     
@@ -621,6 +649,15 @@ def interactive_forcecurve(theworkouts,workstrokesonly=False):
             )
         )
 
+    sourcerange = ColumnDataSource(
+        data = dict(
+            x2575 = angles2575,
+            y2575 = forces2575,
+            x0595 = angles0595,
+            y0595 = forces0595,
+            )
+        )
+
     plot = Figure(tools=TOOLS,
                   toolbar_sticky=False,toolbar_location="above")
 
@@ -654,15 +691,13 @@ def interactive_forcecurve(theworkouts,workstrokesonly=False):
     avf = Span(location=averageforceav,dimension='width',line_color='blue',
                line_dash=[6,6],line_width=2)
     
-    plot.patch('x0199','y0199',source=sourcerange,color="red",alpha=0.05)
+    plot.patch('x0595','y0595',source=sourcerange,color="red",alpha=0.05)
     plot.patch('x2575','y2575',source=sourcerange,color="red",alpha=0.2)
     plot.line('x','y',source=source,color="red")
     plot.circle('xslip','yslip',source=sourceslipwash,color="red")
 
-    plot.circle('peakforceangle','peakforce',source=sourcepoints,color='cyan')
+    plot.circle('peakforceangle','peakforce',source=sourcepoints,color='black',alpha=0.1)
     plot.line('x','y',source=source,color="red")
-    plot.line('x','y',source=sourcetrend,color="blue")
-    plot.line('x','y',source=sourcefit,color="green")
 
     plot.add_layout(avf)