Merge pull request #109 from moorepants/modernize-period-fitting

Update calculating periods from mat files to work with Python 3 and other dependency changes.

.gitignore

@@ -31,7 +31,3 @@ pip-log.txt
 
 #Mr Developer
 .mr.developer.cfg
-
-# matlab
-# mat files
-*.mat

bicycleparameters/com.py

@@ -91,7 +91,7 @@ def center_of_mass(slopes, intercepts):
     # for each line intersection...
     for j, row in enumerate(comb):
         sl = np.array([slopes[row[0]], slopes[row[1]]])
-        a = unumpy.array(np.vstack((-sl, np.ones((2)))).T)
+        a = unumpy.matrix(np.vstack((-sl, np.ones((2)))).T)
         b = np.array([intercepts[row[0]], intercepts[row[1]]])
         lineX[j] = np.dot(a.I, b)
     com = np.mean(lineX, axis=0)

bicycleparameters/geometry.py

@@ -82,11 +82,11 @@ def calculate_abc_geometry(h, d):
     '''
     # extract the values
     h1, h2, h3, h4, h5 = h
-    d1, d2, d3, d4, d = d
+    d1, d2, d3, d4, d_ = d
     # get the perpendicular distances
     a = h1 + h2 - h3 + .5 * d1 - .5 * d2
     b = h4 - .5 * d3 - h5 + .5 * d4
-    c = umath.sqrt(-(a - b)**2 + (d + .5 * (d2 + d3))**2)
+    c = umath.sqrt(-(a - b)**2 + (d_ + .5 * (d2 + d3))**2)
     return a, b, c
 
 

bicycleparameters/inertia.py

@@ -114,7 +114,7 @@ def inertia_components(jay, beta):
     '''
     sb = unumpy.sin(beta)
     cb = unumpy.cos(beta)
-    betaMat = unumpy.array(np.vstack((cb**2, -2 * sb * cb, sb**2)).T)
+    betaMat = unumpy.matrix(np.vstack((cb**2, -2 * sb * cb, sb**2)).T)
     eye = np.squeeze(np.asarray(np.dot(betaMat.I, jay)))
     return eye
 

bicycleparameters/io.py

@@ -79,23 +79,23 @@ def load_pendulum_mat_file(pathToFile):
 
     '''
     pendDat = {}
-    loadmat(pathToFile, mdict=pendDat)
+    loadmat(pathToFile, mdict=pendDat, squeeze_me=True, chars_as_strings=True)
     # clean up the matlab imports
-    del(pendDat['__globals__'], pendDat['__header__'], pendDat['__version__'])
-    for k, v in pendDat.items():
-        try:
-            # change to an ascii string
-            pendDat[k] = v[0].encode('ascii')
-        except:
-            # if an array of a single number
-            if np.shape(v)[0] == 1:
-                pendDat[k] = v[0][0]
-            # else if the notes are empty
-            elif np.shape(v)[0] == 0:
-                pendDat[k] = ''
-            # else it is the data which needs to be a one dimensional array
-            else:
-                pendDat[k] = v.reshape((len(v),))
+    del pendDat['__globals__']
+    del pendDat['__header__']
+    del pendDat['__version__']
+    # If notes is empty it loads like `array([], dtype='<U1')`, so make it an
+    # empty string.
+    if not isinstance(pendDat['notes'], str):
+        if pendDat['notes'].shape == (0,):
+            pendDat['notes'] = ''
+    # If a time array is present, it may be variable sample rate so delete the
+    # sampleRate and duration entries if they are there.
+    if 'time' in pendDat:
+        if 'sampleRate' in pendDat:
+            del pendDat['sampleRate']
+        if 'duration' in pendDat:
+            del pendDat['duration']
     return pendDat
 
 

bicycleparameters/main.py

@@ -1171,7 +1171,7 @@ def eig(self, speeds):
 
         Parameters
         ----------
-        speeds : ndarray, shape (n,) or float
+        speeds : array_like, shape (n,) or float
             The speed at which to calculate the eigenvalues.
 
         Returns
@@ -1188,11 +1188,10 @@ def eig(self, speeds):
         model.
 
         '''
-        # this allows you to enter a float
-        try:
-            speeds.shape
-        except AttributeError:
+        if isinstance(speeds, float):
             speeds = np.array([speeds])
+        else:
+            speeds = np.asarray(speeds)
 
         par = io.remove_uncertainties(self.parameters['Benchmark'])
 
@@ -1209,6 +1208,7 @@ def eig(self, speeds):
 
         return evals, evecs
 
+
     def plot_eigenvalues_vs_speed(self, speeds, fig=None, generic=False,
                                   color='black', show=False, largest=False,
                                   linestyle='-', grid=False, show_legend=True):

bicycleparameters/period.py

@@ -87,22 +87,9 @@ def calc_periods_for_files(directory, filenames, forkIsSplit):
 
     periods = {}
 
-    def pathParts(path):
-        '''Splits a path into a list of its parts.'''
-        components = []
-        while True:
-            (path, tail) = os.path.split(path)
-            if tail == "":
-                components.reverse()
-                return components
-            components.append(tail)
-
-    pathToRawDataParts = pathParts(directory)
-    pathToRawDataParts.pop()
-    pathToBicycleDir = os.path.join(pathToRawDataParts[0],
-                                    pathToRawDataParts[1],
-                                    pathToRawDataParts[2])
-    pathToPlotDir = os.path.join(pathToBicycleDir, 'Plots', 'PendulumFit')
+    # directory is /path/to/data/bicycles/BikeName/RawData
+    path_to_bicycle_dir = os.sep.join(directory.split(os.sep)[:-1])
+    pathToPlotDir = os.path.join(path_to_bicycle_dir, 'Plots', 'PendulumFit')
 
     # make sure there is a place to save the plots
     if not os.path.exists(pathToPlotDir):
@@ -116,12 +103,17 @@ def pathParts(path):
         # generate a variable name for this period
         periodKey = get_period_key(matData, forkIsSplit)
         # calculate the period
-        sampleRate = get_sample_rate(matData)
         pathToPlotFile = os.path.join(pathToPlotDir,
                                       os.path.splitext(f)[0] + '.png')
-        period = get_period_from_truncated(matData['data'],
-                                           sampleRate,
-                                           pathToPlotFile)
+        if 'time' in matData:
+            period = get_period_from_truncated(matData['data'],
+                                               matData['time'],
+                                               pathToPlotFile)
+        else:
+            sampleRate = get_sample_rate(matData)
+            period = get_period_from_truncated(matData['data'],
+                                               sampleRate,
+                                               pathToPlotFile)
         print("The period is:", period, "\n")
         # either append the the period or if it isn't there yet, then
         # make a new list
@@ -212,16 +204,17 @@ def fit_goodness(ym, yp):
     return rsq, SSE, SST, SSR
 
 
-def get_period(data, sampleRate, pathToPlotFile):
+def get_period(data, sample_rate_or_time, pathToPlotFile):
     '''Returns the period and uncertainty for data resembling a decaying
     oscillation.
 
     Parameters
     ----------
-    data : ndarray, shape(n,)
+    data : array_like, shape(n,)
         A time series that resembles a decaying oscillation.
-    sampleRate : int
-        The frequency that data was sampled at.
+    sample_rate_or_time : int or array_like, shape(n,)
+        Either the frequency in Hertz that data was sampled at or a time array
+        that corresponds to ``data``.
     pathToPlotFile : string
         A path to the file to print the plots.
 
@@ -233,7 +226,12 @@ def get_period(data, sampleRate, pathToPlotFile):
     '''
 
     y = data
-    x = np.linspace(0., (len(y) - 1) / float(sampleRate), num=len(y))
+    if isinstance(sample_rate_or_time, int):
+        sample_rate = sample_rate_or_time
+        x = np.linspace(0.0, (len(y) - 1)/float(sample_rate), num=len(y))
+    else:
+        x = sample_rate_or_time
+        sample_rate = int(1.0/np.mean(np.diff(x)))  # approximate
 
     def fitfunc(p, t):
         '''Decaying oscillation function.'''
@@ -246,7 +244,7 @@ def fitfunc(p, t):
     # initial guesses
     # p0 = np.array([1.35, -.5, -.75, 0.01, 3.93]) # guess from delft
     # p0 = np.array([2.5, -.75, -.75, 0.001, 4.3]) # guess from ucd
-    p0 = make_guess(data, sampleRate)  # tries to make a good guess
+    p0 = make_guess(data, sample_rate)  # tries to make a good guess
 
     # create the error function
     errfunc = lambda p, t, y: fitfunc(p, t) - y
@@ -284,9 +282,9 @@ def fitfunc(p, t):
     T = 1. / fd
 
     # plot the data and save it to file
-    fig = plt.figure()
+    fig, ax = plt.subplots(layout='constrained')
     plot_osfit(x, y, lscurve, p1, rsq, T, m=np.max(x), fig=fig)
-    plt.savefig(pathToPlotFile)
+    fig.savefig(pathToPlotFile)
     plt.close()
 
     # return the period
@@ -296,7 +294,10 @@ def fitfunc(p, t):
 def get_period_from_truncated(data, sampleRate, pathToPlotFile):
     # dataRec = average_rectified_sections(data)
     dataRec = data
-    dataGood = select_good_data(dataRec, 0.1)
+    if isinstance(sampleRate, int):
+        dataGood = select_good_data(dataRec, 0.1)
+    else:  # don't truncate if a time array is given in the file
+        dataGood = dataRec
     return get_period(dataGood, sampleRate, pathToPlotFile)
 
 
@@ -456,7 +457,7 @@ def plot_osfit(t, ym, yf, p, rsq, T, m=None, fig=None):
         The measured voltage
     yf : ndarray (n,)
     p : ndarray (5,)
-        The fit parameters for the decaying osicallation fucntion
+        The fit parameters for the decaying oscillation function.
     rsq : float
         The r squared value of y (the fit)
     T : float
@@ -469,38 +470,20 @@ def plot_osfit(t, ym, yf, p, rsq, T, m=None, fig=None):
     fig : the figure
 
     '''
-    # figure properties
-    figwidth = 4.  # in inches
-    goldenMean = (np.sqrt(5) - 1.0) / 2.0
-    figsize = [figwidth, figwidth * goldenMean]
-    params = {#'backend': 'ps',
-        'axes.labelsize': 8,
-        'axes.titlesize': 8,
-        'text.fontsize': 8,
-        'legend.fontsize': 8,
-        'xtick.labelsize': 6,
-        'ytick.labelsize': 6,
-        'text.usetex': True,
-        #'figure.figsize': figsize
-        }
     if fig:
-        fig = fig
+        ax1 = fig.axes[0]
     else:
-        fig = plt.figure(2)
-    fig.set_size_inches(figsize)
-    plt.rcParams.update(params)
-    ax1 = plt.axes([0.125, 0.125, 0.9-0.125, 0.65])
-    #if m == None:
-        #end = len(t)
-    #else:
-        #end = t[round(m/t[-1]*len(t))]
-    ax1.plot(t, ym, '.', markersize=2)
+        fig, ax1 = plt.subplots(layout='constrained')
+    ax1.plot(t, ym, '.', markersize=4)
     plt.plot(t, yf, 'k-')
     plt.xlabel('Time [s]')
     plt.ylabel('Amplitude [V]')
-    equation = r'$f(t)={0:1.2f}+e^{{-({3:1.3f})({4:1.1f})t}}\left[{1:1.2f}\sin{{\sqrt{{1-{3:1.3f}^2}}{4:1.1f}t}}+{2:1.2f}\cos{{\sqrt{{1-{3:1.3f}^2}}{4:1.1f}t}}\right]$'.format(p[0], p[1], p[2], p[3], p[4])
-    rsquare = '$r^2={0:1.3f}$'.format(rsq)
-    period = '$T={0} s$'.format(T)
+    equation = (r'$f(t)={0:1.2f}+e^{{-({3:1.3f})({4:1.1f})t}}\left[{1:1.2f}'
+                r'\sin{{\sqrt{{1-{3:1.3f}^2}}{4:1.1f}t}}+{2:1.2f}'
+                r'\cos{{\sqrt{{1-{3:1.3f}^2}}{4:1.1f}t}}\right]$')
+    equation = equation.format(p[0], p[1], p[2], p[3], p[4])
+    rsquare = r'$r^2={0:1.3f}$'.format(rsq)
+    period = r'$T={0} s$'.format(T)
     plt.title(equation + '\n' + rsquare + ', ' + period)
     plt.legend(['Measured', 'Fit'])
     if m is not None:

bicycleparameters/tests/test_io.py

@@ -0,0 +1,26 @@
+import os
+import numpy as np
+from bicycleparameters.io import load_pendulum_mat_file
+
+TESTS_DIR = os.path.dirname(__file__)
+
+
+def test_load_pendulum_mat_file():
+    path = os.path.join(TESTS_DIR, 'sample_data',
+                        'BrowserForkCompoundFirst1.mat')
+
+    d = load_pendulum_mat_file(path)
+
+    assert d['angle'] == 'First'
+    assert d['ActualRate'] == 1000
+    assert d['bicycle'] == 'BrowserIns'
+    assert d['data'].shape == (30000,)
+    np.testing.assert_allclose(
+        d['data'][0:5],
+        np.array([0.33034183, 0.32525861, 0.31509219, 0.31509219, 0.33034183]))
+    assert d['duration'] == 30
+    assert d['filename'] == 'BrowserInsForkCompoundFirst1'
+    assert d['notes'] == ''
+    assert d['part'] == 'Fork'
+    assert d['pendulum'] == 'Compound'
+    assert d['trial'] == '1'

conda/bicycleparameters-dev.yml

@@ -17,6 +17,7 @@ dependencies:
   - pandas >=1.3.5
   - plotly
   # dev
+  - ipython
   - nose
   - numpydoc >=1.2
   - pytest

docs/data.rst

@@ -211,14 +211,31 @@ Notes
 
 Pendulum Data Files
 ===================
+
 If you have raw signal data that the periods can be estimated from, then these
 should be included in the ``RawData`` directory. There should be at least one
 file for every period typically found in ``<bicycle name>Measured.txt`` file. The
 signals collected should exhibit very typical decayed oscillations. Currently
 the only supported file is a Matlab mat file with these variables:
 
-- ``data`` : signal vector of a decaying oscillation
-- ``sampleRate`` : sample rate of data in hertz
+- ``bicycle``, char : short name of the bicycle, e.g. ``'Browser'``
+- ``part``, char : examples are ``'Fork'``, ``'Frame'``, ``'Rwheel'``
+- ``pendulum``, char : either ``'Compound'`` or ``'Torsional'``
+- ``angle`` or ``angleOrder``, char : the ``'First'``, ``'Second'``, ...,
+  ``'Sixth'`` orientation
+- ``trial``, char : the integer number of the repetition ``'1'``, ``'2'``, etc.
+- ``filename``, char : filename constructed from above variables (see below for
+  explanation)
+- ``notes``, char : any notes recorded for that measurement
+- ``data``, double size n x 1 : signal vector of a decaying oscillation (from
+  rate gyro)
+- Time information, either ``time`` OR both ``duration`` and ``sampleRate``:
+
+   - ``duration``, double size 1 x 1 : time in seconds for the measurement
+   - ``sampleRate`` or ``ActualRate``, double size 1 x 1 (integer) : sample
+     rate of data in Hertz
+   - ``time``, double size n x 1 : time values in seconds that correspond to
+     ``data``
 
 The files should be named in this manner ``<short
 name><part><pendulum><orientation><trial>.mat`` where:
@@ -230,6 +247,11 @@ name><part><pendulum><orientation><trial>.mat`` where:
   ``Fifth``, or ``Sixth``
 - ``<trial>`` is an integer greater than or equal to 1
 
+These data files were originally collected with Matlab and a National
+Instruments USB-6218 using a Silicon Sensing CRS03-04S single axis angular rate
+gyro and a `data collection script
+<https://github.com/moorepants/PhysicalParameters/blob/master/PhysicalParameters/matlab/acquire_data.m>`_.
+
 Notes
 -----
 
@@ -273,7 +295,7 @@ files.
 riders/<rider name>/RawData/
 ============================
 
-**These files must follow the YAML format used in the `yeadon package`_.**
+**These files must follow the YAML format used in the yeadon package.**
 
 <rider name><bicycle name>YeadonCFG.txt
 ---------------------------------------
@@ -291,4 +313,4 @@ all of the geometric measurements of the rider. See the `yeadon documentation`_
 for more details.
 
 .. _yeadon package: http://pypi.python.org/pypi/yeadon
-.. _yeadon documentation : http://packages.python.org/yeadon/
+.. _yeadon documentation : https://yeadon.readthedocs.io