Speed up inverse dynamics (Python version) #128
Conversation
gaitanalysis/gait.py
Outdated
| @@ -145,6 +146,29 @@ def __init__(self, data): | |||
| f.close() | |||
| self.load(data) | |||
|
|
|||
| def _leg2d(self, time, marker_values, normalized_force_plate_values, | |||
There was a problem hiding this comment.
This could be a static method or pulled out as a function.
|
So this all seems to work. It gives the same results as the Octave file on the test data but when I run it for the plot in our perturbed data paper I'm seeing some noise in the joint torques in the Python version that I don't see in the Octave version. I must have some slight differences in the filtering/differentiation code. Here is a screenshot: @tvdbogert Can you comment on the filtering/differentiation? You have a custom filterer differentiate. I switched to using a 2nd order Butterworth filter, with filtfilt, and central differencing. |
|
I changed it so I filter the marker positions, compute the derivative, filter the velocities, and compute the derivative instead of filter the marker positions, compute the derivative, and compute the derivative and I get results more similar to the Octave file. FYI: Images are swapped here with respect to the above. |
|
Good to know that Octave performs so poorly on this code. I looked in your speedup inverse dynamics branch but could not locate the python code. Pls give a link if possible. I think I know what is going on. My rtfilter.m processes data with time stamps, one sample at a time. It can handle gaps in the data. If you use standard digital butterworth filters, you need a constant sampling rate. The spikes may be generated if you do the differentiation with the assumption of constant sampling rate. If a frame drops out, you then suddenly get only half of the true derivative. The second derivative then gets spikes. This may happen in a certain part of the gait cycle when the hip marker drops out sometimes. Filtering again after the differentiation will get rid of it, but that just hides the problem and in fact distorts the signal. In fact, you should never filter more than once in this sequence of steps. Filtering and differentiation are both linear operations on the signal, and you can change their order. So filter only once, at the beginning or at the end. My filter has the advantage that the filtering and differentiation is combined into one step, which reduces the distortion from finite differencing. The for loop in myfiltfilt.m is probably the bottleneck for Octave. It may be possible to vectorize that, I would have to take a real good look at rtfilter.m. The problem is that the filter coefficients depend on the time between two samples. It may be possible to code that without loops. This code was written for real-time (I originally wrote it in C, Motek has that code in D-Flow) so it processes one sample at a time. I would like to keep the Matlab code around, with the Python equivalent in a separate folder. In that case, the Python code should produce exactly the same output. Here are the options as I see them:
Whatever we do, the Matlab/Octave and Python code should give identical results. |
inverse dynamics: https://github.com/csu-hmc/GaitAnalysisToolKit/pull/128/files#diff-5e1f4c39e8343cdd4b4dedb9eda0a0efR1084 mimics leg 2d: https://github.com/csu-hmc/GaitAnalysisToolKit/pull/128/files#diff-5e1f4c39e8343cdd4b4dedb9eda0a0efR149 numerical differences: https://github.com/moorepants/DynamicistToolKit/blob/master/dtk/process.py#L575 filter: https://github.com/moorepants/DynamicistToolKit/blob/master/dtk/process.py#L462
That may be it, but the
I thought that too, but oddly filter twice changes the result in this case. I must have a slight error somewhere.
Yes, it definitely is. See the timings here: #30 That loop would ideally be in a low level language. I'm not sure if vectorizing it would help that much.
We can do that. I currently have a unit test that tests my code against yours: The results for the
Yes, this could be done. But I'm pretty sure that code needs to be written in C or something for speed purposes. A direct translation will likely have some speed gains relative to Octave but it want be the 1000x speedup I'm seeing using the constant sample rate filtering here.
This is likely the best solution.
This is what is implemented now, in this PR.
Agreed. |
|
See if my revised Matlab/Octave code helps at all. It eliminates the If you already interpolated before the data gets to leg2d, my Also if you already interpolated we can safely do a standard butterworth In leg2d you then need to do a check: if (max(diff(timestamps)) ~= min(diff(timestamps)) My test data in rawdata.txt has missing samples, so that would not be a One thing you may not have noticed in myfiltfilt.m. If you run a 2nd A C version would be good to have. On 3/30/2015 12:49 PM, Jason K. Moore wrote:
|
Testing now.
I'll investigate more.
For my unit test I pre-interpolate the missing samples and pass that into both the Python and Octave code to see if I get the same results.
I didn't know this. Maybe that's the difference. I'll fix my code.
I can implement the C version if you are fine with including your method in this package. Just let me know. |
moorepants
changed the title
@tvdbogert Where do I find the theory for this? I've never heard of it and would like to implement a generalized version. |
|
I found some reference here for 2nd order filters: Do you know if there is a general formula for any order Butterworth filter? |
|
In the code I cite the Winter book which explains this, but it's not If you apply a filter H(w) twice, the transfer function is H(w)^2. Nth order Butterworth transfer function, magnitude, is: |H(w)| = 1/sqrt(1+(w/w0)^(2n)) w0 is the cutoff frequency, where transfer magnitude is 1/sqrt(2). applied twice: |H(w)| = 1/(1+(w/w0)^(2n)) If you define cutoff frequency as the frequency w0 where transfer is 1/(1+(w0/w0corrected)^(2n)) = 1/sqrt(2) Where w0corrected is the design parameter you need to put in the w0corrected = (sqrt(2) - 1)^(1/2n) * w0 = 0.802 * w0 (when n=2) If you don.t do this, and you apply a 6 Hz lowpass filter twice, you end On 3/30/2015 10:48 PM, Jason K. Moore wrote:
|
|
fThat's a very nice notebook. Good to see that he used Pezzack's He gives a formula for any filter order n, but I get a slightly On 3/30/2015 11:12 PM, Jason K. Moore wrote:
|
He's creating a whole note set (or textbook) here: https://github.com/demotu/BMC
I've adjusted my Note that your equation above is slightly incorrect, it should be: w0corrected = wo / (sqrt(2) - 1)^(1/2n) = w0 / 0.802 (when n=2) |
|
Ok, I now get virtually identical results with these changes:
|
|
That's good enough agreement in results! I assume you're using the conventional digital Butterworth filter, which eliminates the need for a C version. As long as you're interpolating the raw data onto a fixed sampling rate. Which is OK if the gaps are small enough. And yes, I had my correction factor backwards, thanks for the correction. |
Yep, I'm traditional like that :)
I'll implement the C version of your filter if I get an itch. It is a more elegant solution.
Your filter actually had issues with large marker gaps. The data in #132 exposed this. I assume there has got to be a limit to the marker gap that makes either method fail. |
moorepants
force-pushed
the
speedup-inverse-dynamics
branch
from
9c78cf0 to
1363ad8
Compare
The gait/interpolate function is only used in the GaitData.split_at() method and the way I orginally coded it using Pandas paradigms was very slow. This change speeds things up by about 40X. I had to change the test because it no longer supports nans in the DataFrame. Being that it is only used in the split_at method this seemed ok. So far I've never had nans in the data frame by the time it gets to the split at method. GaitData isn't built to handle nans, they should be fixed before creating a GaitData object.
moorepants
force-pushed
the
speedup-inverse-dynamics
branch
from
1363ad8 to
9e80dfd
Compare
Addresses issues #30 and some of #49.
The Octave file
leg2d.mwas quite slow. It takes from 5 to 6 minutes on my machine to compute the inverse dynamics of a 8 minute (100 Hz) chunk of data. This has bothered me a long time because it takes hours to process all of our data when it shouldn't. So I decided to finally fix it.I decided not to spend time improving the performance of the
leg2d.mand its dependent functions because I'd rather purge the Octave stuff from this library (the octave/python bridge has too much overhead and complicates installation especially on Windows).So I've simply ported the
leg2d.mfile to a Python function ingait.pyand also switched to using library filters (scipy.signal) and the finite difference functions in DynamicistToolKit for speed gains. From the simple port I'm seeing close to a 1000 fold performance increase. It now takes less than a second to process 8 minutes of data.lower_extremity_2d_inverse_dynamics()into several smaller functions or convert to a class.