Merge pull request #4039 from opensim-org/path_fitter_segfault_fix

Fix bugs in `PolynomialPathFitter` related to the number of input coordinate samples

Author: nickbianco

Bindings/Java/Matlab/examples/PolynomialPathFitter/examplePolynomialPathFitter.m

@@ -50,16 +50,35 @@
 %
 % The fitter will randomly sample around the coordinate values provided in the
 % table to generate model configurations for which to compute path lengths and
-% moment arms. This table has many more rows than are needed for the fitter to
-% generate a good fit, so we will remove some of the rows to speed up the
-% fitting process.
+% moment arms. In general, it is up to the user to decide how many sample points
+% are needed to adequately cover the range of motion of the model's coordinates.
+% In this case, the table has many more rows than are needed for the fitter to
+% generate a good fit, so we will remove some of the rows to speed up the fitting
+% process.
 values = TimeSeriesTable('coordinates.sto');
 times = values.getIndependentColumn();
 for i = times.size():-1:1
     if mod(i, 10) ~= 0
         values.removeRowAtIndex(i);
     end
 end
+
+% Add columns for the toe joints. Use a sine function to generate the
+% coordinate values. An amplitude of 0.5 keeps the values within the coordinate
+% range of the MTP joints in the model.
+amplitude = 0.5;
+angular_frequency = 10.0; % rad/s
+phase = 0.0;
+sine = Sine(amplitude, angular_frequency, phase);
+mtp_values = Vector(values.getNumRows(), 0.0);
+time = Vector(1, 0.0);
+for i = 0:values.getNumRows()-1
+    time.set(0, times.get(i));
+    mtp_values.set(i, sine.calcValue(time));
+end
+
+values.appendColumn('/jointset/mtp_r/mtp_angle_r/value', mtp_values);
+values.appendColumn('/jointset/mtp_l/mtp_angle_l/value', mtp_values);
 fitter.setCoordinateValues(TableProcessor(values));
 
 % Configure optional settings

Bindings/Python/examples/PolynomialPathFitter/examplePolynomialPathFitter.py

@@ -50,15 +50,33 @@
 #
 # The fitter will randomly sample around the coordinate values provided in the
 # table to generate model configurations for which to compute path lengths and
-# moment arms. This table has many more rows than are needed for the fitter to
-# generate a good fit, so we will remove some of the rows to speed up the
-# fitting process.
+# moment arms. In general, it is up to the user to decide how many sample points
+# are needed to adequately cover the range of motion of the model's coordinates.
+# In this case, the table has many more rows than are needed for the fitter to
+# generate a good fit, so we will remove some of the rows to speed up the fitting
+# process.
 values = osim.TimeSeriesTable('coordinates.sto')
 times = values.getIndependentColumn()
 for i in range(len(times)):
     if i % 10 != 0:
         values.removeRow(times[i])
 
+# Add columns for the toe joints. Use a sine function to generate the
+# coordinate values. An amplitude of 0.5 keeps the values within the coordinate
+# range of the MTP joints in the model.
+amplitude = 0.5
+angular_frequency = 10.0 # rad/s
+phase = 0.0
+sine = osim.Sine(amplitude, angular_frequency, phase)
+mtp_values = osim.Vector(values.getNumRows(), 0.0)
+times = values.getIndependentColumn()
+time = osim.Vector(1, 0.0)
+for i in range(values.getNumRows()):
+    time.set(0, times[i])
+    mtp_values.set(i, sine.calcValue(time))
+
+values.appendColumn('/jointset/mtp_r/mtp_angle_r/value', mtp_values)
+values.appendColumn('/jointset/mtp_l/mtp_angle_l/value', mtp_values)
 fitter.setCoordinateValues(osim.TableProcessor(values))
 
 # Configure optional settings

CHANGELOG.md

@@ -88,6 +88,7 @@ v4.6
   other path-based force elements dependent on a user-provided expression. (#4035)
 - Fixed a bug where `DeGrooteFregly2016Muscle::getBoundsNormalizedFiberLength()` was returning
   tendon force bounds rather than fiber length bounds. (#4040)
+- Fixed bugs in `PolynomialPathFitter` when too few coordinate samples were provided. (#4039)
 
 
 v4.5.1

OpenSim/Actuators/PolynomialPathFitter.cpp

@@ -170,6 +170,12 @@ void PolynomialPathFitter::run() {
             "Expected 'num_parallel_threads' to be between 1 and {}, but "
             "received {}.", std::thread::hardware_concurrency(),
             get_num_parallel_threads())
+    OPENSIM_THROW_IF_FRMOBJ(
+            static_cast<int>(values.getNumRows()) < get_num_parallel_threads(), 
+            Exception, "Expected the number of time points in the coordinate "
+            "values table to be greater than 'num_parallel_threads', but "
+            "received {} and {}, respectively.", 
+            values.getNumRows(), get_num_parallel_threads())
     log_info("Number of parallel threads = {}", get_num_parallel_threads());
 
     // Number of samples per frame.
@@ -662,7 +668,8 @@ TimeSeriesTable PolynomialPathFitter::sampleCoordinateValues(
     int timeIdx = 0;
     const auto& times = values.getIndependentColumn();
     TimeSeriesTable valuesSampled;
-    double dt = (times[1] - times[0]) / (get_num_samples_per_frame() + 2);
+    double dt = (times.size() < 2) ? 0.01 : 
+            (times[1] - times[0]) / (get_num_samples_per_frame() + 2);
     for (int i = 0; i < get_num_parallel_threads(); ++i) {
         int numTimeIndexes = outputs[i].nrow() / get_num_samples_per_frame();
         for (int j = 0; j < numTimeIndexes; ++j) {

OpenSim/Actuators/PolynomialPathFitter.h

@@ -151,6 +151,13 @@ class OSIMACTUATORS_API PolynomialPathFitterBounds : public Object {
  * for the fitting process, or if the fit for a particular path did not meet the
  * specified tolerances.
  *
+ * In general, it is up to the user to decide how many sample points are needed
+ * to adequately cover the range of motion of the model's coordinates. As the
+ * complexity of a muscle path increases, more sample points over a larger
+ * dimension of coordinate values are needed to achieve a good fit. Users may
+ * consider manually creating the coordinates value table to ensure that the
+ * sampling covers the full range of motion for the model.
+ *
  * It is highly recommended to use the files printed to the output directory to
  * evaluate the quality of the fitted paths (see `setOutputDirectory()` for more
  * details). Depending on the quality of the original model, it may not be

OpenSim/Actuators/Test/testPolynomialPathFitter.cpp

@@ -57,8 +57,9 @@ namespace {
         return processor;
     }
 
-    TableProcessor createCoordinatesTable(bool correctLabel, bool addMetaData) {
-        SimTK::Vector column = SimTK::Test::randVector(100);
+    TableProcessor createCoordinatesTable(bool correctLabel, bool addMetaData,
+            int numRows = 100) {
+        SimTK::Vector column = SimTK::Test::randVector(numRows);
         std::vector<double> times;
         times.reserve(column.size());
         for (int i = 0; i < column.size(); ++i) {
@@ -130,3 +131,17 @@ TEST_CASE("Invalid configurations") {
             "Expected 'maximum_polynomial_order' to be at most 9"));
     }
 }
+
+TEST_CASE("Number of rows less than the number of threads") {    
+    int numAvailableThreads = std::thread::hardware_concurrency();
+    int numRows = (numAvailableThreads == 1) ? 1 : numAvailableThreads - 1;
+    CAPTURE(numAvailableThreads);
+    CAPTURE(numRows);
+   
+    PolynomialPathFitter fitter;  
+    fitter.setModel(createHangingMuscleModel());
+    fitter.setCoordinateValues(createCoordinatesTable(true, true, numRows));
+    fitter.setNumParallelThreads(numAvailableThreads);
+    REQUIRE_THROWS_WITH(fitter.run(), ContainsSubstring(
+        "Expected the number of time points in the coordinate values table"));
+}