csu-hmc / gaitanalysistoolkit Goto Github PK
View Code? Open in Web Editor NEWTools for the Cleveland State Human Motion and Control Lab
Home Page: http://hmc.csuohio.edu
License: Other
Tools for the Cleveland State Human Motion and Control Lab
Home Page: http://hmc.csuohio.edu
License: Other
If you try to load in a compensation trial there are some issues.
This takes 23 seconds on my machine for 500 steps.
I'm running the numpy 1.8, pandas 0.13 and this test now fails. Seems to have something to do with the changes in min_time and max_time.
result = solver.solve(sparse_a=True)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-44-c96f1dfa9907> in <module>()
----> 1 result = solver.solve(sparse_a=True)
/home/moorepants/src/Gait-Analysis-Toolkit/gaitanalysis/controlid.pyc in solve(self, sparse_a, gain_omission_matrix)
749 x, variance, covariance = self.least_squares(A, b)
750
--> 751 deconstructed_solution = self.deconstruct_solution(x, covariance)
752
753 gain_matrices = deconstructed_solution[0]
/home/moorepants/src/Gait-Analysis-Toolkit/gaitanalysis/controlid.pyc in deconstruct_solution(self, x, covariance)
299 control_vectors_variance = np.zeros((self.n, self.q))
300
--> 301 parameter_variance = np.diag(covariance)
302
303 for i in range(self.n):
/home/moorepants/envs/walk/local/lib/python2.7/site-packages/numpy/lib/twodim_base.pyc in diag(v, k)
287 return v.diagonal(k)
288 else:
--> 289 raise ValueError("Input must be 1- or 2-d.")
290
291 def diagflat(v, k=0):
ValueError: Input must be 1- or 2-d.
We can wrap Ton's filter that can handle varying sample rates and NaN values and use that as the universal filter.
Note that it is typical to compute inverse dynamics before filtering and filter the results, but we are not doing that currently.
You can filter before or after other linear operations (like differentiation) and get the same effect. But the order matters for non-linear operations, like non-linear inverse dynamics.
We record the forces at the beginning and end of a trial with no load on the treadmill and not treadmill motion. These force measurements can be used to account for any drift. We need to check into that.
If event 'A' is recorded as firing multiple times, we should store it. Either as:
events['A'] = [(start1, stop1), (start2, stop2)]
or
events['A1'] = (start1, stop1)
events['A2'] = (start2, stop2)
The label methods in DFlowData should be moved to a class dedicated to that.
It would be nice to have a method in DFlowData that output all the data expressed in a different coordinate frame. You could pass in a rotation matrix and all of the data that has the xyz coordinates from the mocap system could be output into a different global reference frame.
Sandy,
I think it will be best if we don't use spaces in directories or file names, there are cross platform issues sometimes, and life will likely be easier if we avoid it.
You can change your directories and have git track it using the git mv
command:
git mv Inertial\ Compensation inertial_compensation
Here is a profile of essentially running DFlowData.clean_data() + all the methods in WalkingData (including inverse dynamics compuations).
7641430 function calls (7640347 primitive calls) in 181.271 seconds
Ordered by: internal time
ncalls tottime percall cumtime percall filename:lineno(function)
320218 165.994 0.001 165.994 0.001 {method 'readline' of 'file' objects}
3 2.709 0.903 2.709 0.903 {method 'read' of 'pandas.parser.TextReader' objects}
36844 1.017 0.000 1.017 0.000 {pandas.lib.maybe_convert_objects}
19 0.879 0.046 168.193 8.852 session.py:494(evaluate)
36839 0.816 0.000 4.204 0.000 series.py:3129(interpolate)
37612 0.574 0.000 0.643 0.000 common.py:134(_isnull_ndarraylike)
120 0.530 0.004 0.532 0.004 fitpack2.py:441(__init__)
41185 0.385 0.000 0.385 0.000 {method 'reduce' of 'numpy.ufunc' objects}
320229 0.331 0.000 0.331 0.000 {_codecs.utf_8_decode}
320226 0.284 0.000 0.777 0.000 {method 'decode' of 'str' objects}
268534 0.267 0.000 0.485 0.000 {method 'view' of 'numpy.ndarray' objects}
37998 0.258 0.000 0.851 0.000 series.py:430(__new__)
37612 0.238 0.000 0.988 0.000 common.py:64(_isnull_new)
951420/951418 0.226 0.000 0.226 0.000 {isinstance}
38160 0.208 0.000 0.209 0.000 {method 'copy' of 'numpy.ndarray' objects}
111799 0.208 0.000 0.426 0.000 index.py:332(__getitem__)
36642 0.203 0.000 1.341 0.000 frame.py:1928(_ixs)
320110 0.185 0.000 0.297 0.000 __init__.py:1343(isEnabledFor)
98213 0.169 0.000 0.169 0.000 {numpy.core.multiarray.array}
320223 0.162 0.000 0.493 0.000 utf_8.py:15(decode)
320110 0.154 0.000 0.451 0.000 __init__.py:1128(debug)
151206 0.133 0.000 0.193 0.000 numeric.py:1810(isscalar)
38175 0.132 0.000 0.298 0.000 series.py:3299(_sanitize_array)
801 0.130 0.000 0.130 0.000 {pandas.algos.take_2d_axis0_float64_float64}
164141/164139 0.127 0.000 0.129 0.000 {getattr}
1 0.119 0.119 0.277 0.277 __init__.py:20(<module>)
1246 0.117 0.000 0.117 0.000 {numpy.core.multiarray.concatenate}
36839 0.112 0.000 4.317 0.000 frame.py:4389(<lambda>)
320110 0.112 0.000 0.112 0.000 __init__.py:1329(getEffectiveLevel)
197 0.110 0.001 6.391 0.032 frame.py:4433(_apply_standard)
36642 0.098 0.000 0.647 0.000 internals.py:1681(iget)
76780 0.095 0.000 0.108 0.000 index.py:318(__contains__)
75404 0.092 0.000 0.213 0.000 series.py:532(__array_finalize__)
I think using readline to pull the events from the record files is a major slow down. Also all the calls to session.py are oct2py calls, which are slow. I believe the slow parts are the inverse dynamics real time filter and the soder.m file.
# TODO : The time stamp is removed from dflow_data.data here, so the
# following command fails because there is no timestamp column. i.e.
# the class is mutable. This problematic if you want to extract more
# than one event.
Try extracting two events in a row from DFlowData and setting the index_col flag to the same column on both.
@campanawanna Is this something you added in the recent merge? If so, we need to avoid showing plots when running tests, or at least make them optional.
We need to use lil_matrix to preallocate a sparse matrix and then fill it up.
We should put a license with the code so it is clear to people how they can use it. I've been using http://unlicense.org/ lately which makes it public domain (no restrictions). If you all desire some more restrictive (require citation, non-commerical, etc) we should discuss that. Right now the UNLICENSE is in the directory as that is what I had in DynamicisitToolKit.
octave:8> profile on
octave:9> leg2d(rand(5000, 1), rand(5000, 12), rand(5000, 3), options);
Marker 1: 0 samples are missing, longest gap is 0 samples.
Marker 2: 0 samples are missing, longest gap is 0 samples.
Marker 3: 0 samples are missing, longest gap is 0 samples.
Marker 4: 0 samples are missing, longest gap is 0 samples.
Marker 5: 0 samples are missing, longest gap is 0 samples.
Marker 6: 0 samples are missing, longest gap is 0 samples.
octave:10> profile off
octave:11> T = profile('info');
octave:12> profshow(T)
# Function Attr Time (s) Calls
-------------------------------------------------------------
26 rtfilter 3.560 70000
25 myfiltfilt>rtfilter_batch 2.287 14
3 binary * 0.306 947201
12 binary + 0.138 315762
10 zeros 0.112 69898
24 binary / 0.104 315715
29 pi 0.074 70162
5 size 0.072 70001
11 binary - 0.065 175555
30 binary ^ 0.060 175388
27 isempty 0.046 70000
28 binary <= 0.042 70018
31 sqrt 0.022 35081
2 leg2d 0.005 1
33 lookup 0.004 18
1 rand 0.004 3
32 interp1q 0.003 18
18 myfiltfilt 0.001 7
19 flipud 0.001 21
44 unwrap 0.001 4
octave:13>
The do_plot
flag in grf_landmarks no longer calls show() and it now plots two windows instead of subplots. The number of steps is missing from the axes titles.
@campanawanna Let's agree on some of these changes before merging. You must have added this after I reviewed the PR for your acceleration stuff.
In test_markers.py
np.linalg.det() doesn't work on three-dimensional arrays in numpy 1.7.1
Now that we have a stationary V-gait we do not need the inertial compensation with the accelerometers (although it would still help, once we get better accelerometers).
But we do need to to the compensation for the belt/roller acceleration. Sandy has found a linear relationship between Mx (Nm) and the acceleration (m s^-2). Theoreticaly, no other force plate signals should be affected by belt acceleration. Sandy can probably confirm that with her data.
I would suggest to do a 6 Hz zero-lag filter on the belt speed signal (make sure it's exactly the same filter that is used for all other signals), then do a 3-point zero lag derivative estimation (xdot_i = (x_i+1 - x_i-1)/(t_i+1 - t_i-1)) to get the acceleration. Multiply by the correction coefficient, subtract result from Mx and that's it.
If you use filtered velocity, you should correct after filtering the force plate signals, otherwise the correction signal gets filtered twice. Since all these processing steps are linear, I suppose you can use unfiltered velocity and correct before filtering the force plate signals. The order does not matter. The only thing is, if you look at a corrected but unfiltered Mx signal, it may look much noisier than the original Mx signal because it has correction related to a possibly noisy belt acceleration signal. At the end it's all the same, it only matters if you want to inspect these intermediate results.
In test_motek all sample data in self.mocap_data_frame is essentially the same sine wave. It would be nice for each column to be different.
Seems more intuitive.
solver = SimpleControlSolver(data)
solver.solve(['sensor1', 'sensor2'], ['control1'], gain_omission_matrix=np.array(2 * [True]))
I think it was a poor choice to use a time vector of floats for the index in walking classes. It may be possible to use Panda's datetime index classes, but accuracy issues arise because it can only resolve to nanoseconds (based on integers from Epoch). I got the float indices to work for the WalkingData class and the SimpleController class, but I should probably revert to an integer index and move the time stamp's into the data columns. This doesn't work as cleanly as one would hope but there are fundamental issues with using floating point indices.
If num_samples=None
in the split_at method, then the shortest step (fewest number of samples) is used for the output number of samples. This works fine if all steps are similar length but if one step has incorrect heelstrike times then you can get steps with too few data points.
There have been measurable delays when comparing accelerometer readings to corresponding force plate readings and also comparing marker readings to force plate readings. We need to ensure that these are taken care of in the D-Flow/Cortex software pipeline and if not incorporate these delays into our processing pipleline.
My logic is messed up in the inverse 2d dynamics function in WalkingData.
This is an issue with the D-Flow data, so ideally this would be a method on DFlowData but it will likely require the identification of steps so that you can isolate the stance phase of each step to find a mean offset for each.
Should this happen before or after inertial compensation? I guess it has to happen after inertial compensation because the forces and moments will not be zero.
Needs preallocation here:
and here:
Stacking arrays in loops might not be the best either:
I think these must be HBM outputs that give the segment rotations relative to something.
If we plot the gains from 0 to 100% then first and last sample can't be both at 0% and 100% because that is the same time, i.e. heel strike.
If the delay is specified we can solve a linear system such as:
m(t) = m0(t) + K(t) * [s0(t)- s(t)] + Kd * [s0(t-tau) - s(t-tau)]
I think the gain_omission_matrix should probably be named gain_inclusion_matrix. I may have goofed up the name with respect to the definition.
In Octave-Matlab-Codes/2D-Inverse-Dynamics/test/rawdata.txt
the force and moments are labeled FP1.*
where as the markers are R*
. The right force plate should be called FP2. I don't think there is an issue with the data, i.e. all the data is from the same side, but the labels are just not consistent. Still haven't found the real issue yet...
Probably should keep track of that.
@tvdbogert @spinningplates @campanawanna
I'm going to shoot for submitting a paper on what's been done so far around July 1st, 2014. Both the data and the source code (including this repository) for processing it will be released alongside the paper publication. This means that a bit of work is needed to get this repository in shape for it's initial release. So I'm going to introduce some standards for code quality. These are open for discussion now if you have any concerns. This article, http://dx.doi.org/10.1371/journal.pbio.1001745 is a good starting point to read about reproducible scientific software.
So here are the standards I'd like us to meet for the publication:
Optional, but nice:
Notes:
I've created a milestone to track the issues that must be fixed before the first release.
https://github.com/csu-hmc/Gait-Analysis-Toolkit/issues/milestones
I had this test fail once saying that numpy couldn't assign 160000 values to something like 15692 values produced by a boolean mask. This could possible be a case that depends on the random test values provided. Not sure. And not sure how to reproduce unless I keep running the tests and hit it again.
I'm not sure the WalkingData class is that useful in it self. Could do:
data_frame = gait.inverse_dynamics_2D(data_frame)
data_frame = gait.grf_landmarks(data_frame, left, right)
data_panel = gait.split_at(data_frame, ...)
We basically have that now, but don't need the first arg because that is self
in the class.
This can be useful for plot scaling. In particular in plot_gains
Should list all contributors and their email.
Right now the classes have a weird combination of modifying attributes of the object in-place and methods that don't modifying attributes but return the results, of which other methods use to modify attributes.
This inconsistency is tied around thinking about object immutability. Some objects work nicely as immutable, i.e. their methods always return a new modified object. But mutable objects let you modify the object in place.
pandas data frames are not immutable, you can change the underyling data, but it seems like most of their methods return a new data frame as opposed to modifying it in place. So it is a combination. I'm really not sure what the best OO design paradigm is for this stuff.
Octave seems to have some unit test framework: http://www.gnu.org/software/octave/doc/interpreter/Test-Functions.html
FYI, the release candidate version of D-Flow that we are required to use due to the new stabilization struts has backward incompatible outputs. Missing markers used to be indicated by holding the previous marker value constant. See: http://gait-analysis-toolkit.readthedocs.org/en/latest/motek.html#missing-values. After failing to process the data from the most recent trials I tracked down the reason our software no longer gives correct answers. It seems that the latest DFlow (3.16.2rc4) now outputs a string in the text file representing zero, '0.000000', for missing markers. This does not work with the current version of Gait-Analysis-Toolkit and some refactoring is now needed.
Right now the event detection method in DFlowData looks at the end of the file and counts up to th -6 line to find the last time value. This needs to support any number of events instead of just a fixed number.
When running inertial_compensation.m
on the compensation data needed for run T006 the backslash operator fails to solve the linear least squares here:
This failure gives this error:
octave:7> new_forces = inertial_compensation(calib_forces, calib_accel, markers, forces, accel);
** On entry to DLASD4 parameter number -1 had an illegal value
error: inertial_compensation: exception encountered in Fortran subroutine dgelsd_
error: called from:
error: /home/moorepants/src/Gait-Analysis-Toolkit/Octave-Matlab-Codes/Inertial-Compensation/inertial_compensation.m at line 88, column 3
This error only occurs if the ATLAS LAPACK version is enabled and does not occur when the normal LAPACK is enabled. It would be nice to output the matrices at line 88 and then just try the backslash with ATLAS LAPACK and normal LAPACK.
Loading compensation data is very similar to loading a regular trial. It may make sense to subclass DFlowData and put all the associated compensation computations with that class instead of having them buried in DFlowData. This would also allow the missing marker identification issue in #69 to be addressed for the compensation files.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.