Merge branch 'main' into opensim_451

Author: aymanhab

Bindings/OpenSimHeaders_moco.h

@@ -7,6 +7,8 @@
 #include <OpenSim/Moco/MocoBounds.h>
 #include <OpenSim/Moco/MocoCasADiSolver/MocoCasADiSolver.h>
 #include <OpenSim/Moco/MocoControlBoundConstraint.h>
+#include <OpenSim/Moco/MocoOutputBoundConstraint.h>
+#include <OpenSim/Moco/MocoStateBoundConstraint.h>
 #include <OpenSim/Moco/MocoFrameDistanceConstraint.h>
 #include <OpenSim/Moco/MocoOutputConstraint.h>
 #include <OpenSim/Moco/MocoGoal/MocoAccelerationTrackingGoal.h>

Bindings/moco.i

@@ -48,6 +48,8 @@ namespace OpenSim {
 %include <OpenSim/Moco/MocoConstraint.h>
 
 %include <OpenSim/Moco/MocoControlBoundConstraint.h>
+%include <OpenSim/Moco/MocoOutputBoundConstraint.h>
+%include <OpenSim/Moco/MocoStateBoundConstraint.h>
 %include <OpenSim/Moco/MocoFrameDistanceConstraint.h>
 %include <OpenSim/Moco/MocoOutputConstraint.h>
 

CHANGELOG.md

@@ -6,12 +6,24 @@ request related to the change, then we may provide the commit.
 
 This is not a comprehensive list of changes but rather a hand-curated collection of the more notable ones. For a comprehensive history, see the [OpenSim Core GitHub repo](https://github.com/opensim-org/opensim-core).
 
+v4.6
+====
+- The performance of `getStateVariableValue`, `getStateVariableDerivativeValue`, and `getModelingOption` was improved in
+  the case where provided string is just the name of the value, rather than a path to it (#3782)
+- Fixed bugs in `MocoStepTimeAsymmetryGoal::printDescriptionImpl()` where there were missing or incorrect values printed. (#3842)
+- Added `ModOpPrescribeCoordinateValues` which can prescribe motion of joints in a model given a table of data. (#3862)
+- Fixed bugs in `convertToMocoTrajectory()` and `MocoTrajectory::resampleWithFrequency()` by updating `interpolate()` to 
+  allow extrapolation using the `extrapolate` flag. Combined with the `ignoreNaNs` flag, this prevents NaNs from 
+  occurring in the output. (#3867)
+- Added `Output`s to `ExpressionBasedCoordinateForce`, `ExpressionBasedPointToPointForce`, and `ExpressionBasedBushingForce` for accessing force values. (#3872)
+
 v4.5.1
 ======
 - Added support for list `Socket`s via the macro `OpenSim_DECLARE_LIST_SOCKET`. The macro-generated method
   `appendSocketConnectee_*` can be used to connect `Object`s to a list `Socket`. In addition, `Component` and Socket have
   new `getConnectee` overloads that take an index to a desired object in the list `Socket` (#3652).
 - Added `ComponentPath::root()`, which returns a `ComponentPath` equivalent to "/"
+- Added `ComponentPath::separator()`, which returns the separator that's placed between elements of the path (i.e. `'/'`)
 - `ComponentPath` is now less-than (`<`) comparable, making it usable in (e.g.) `std::map`
 - `ComponentPath` now has a `std::hash<T>` implementation, making it usable in (e.g.) `std::unordered_map`
 - Added `.clear()` and `.empty()` to `ComponentPath` for more parity with `std::string`'s semantics

CHANGELOG_MOCO.md

@@ -1,6 +1,23 @@
 Moco Change Log
 ===============
 
+1.4.0
+-----
+- 2024-07-26: Added `MocoStateBoundConstraint` and `MocoOutputBoundConstraint` to enable bounding 
+              state variables or output values by one or two `Function`s, similar to 
+              `MocoControlBoundConstraint`.
+
+- 2024-07-22: Added support for `MocoOutputGoal`s and `MocoOutputConstraint`s that are 
+              composed of two `Output`s. This applies to all types of Output goals 
+              (`MocoInitialOutputGoal`, `MocoFinalOutputGoal`, etc.). The two `Output`s 
+              can be combined by addition, subtraction, multiplication, or division.
+
+- 2024-07-08: Fixed a bug in `DeGrooteFregly2016Muscle` where updates to properties 
+              `pennation_angle_at_optimal`, `optimal_fiber_length`, `max_contraction_velocity`, 
+              and `tendon_strain_at_one_norm_force` during parameter optimization did not 
+              affect certain model calculations, and as a result were not changing during
+              optimization.
+
 1.3.1
 -----
 - 2024-07-08: Fixed a bug where deserialization of an OpenSim model with the `Bhargava2004SmoothedMuscleMetabolics`

OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp

@@ -152,16 +152,7 @@ void DeGrooteFregly2016Muscle::extendFinalizeFromProperties() {
             "Pennation angle at optimal fiber length must be in the range [0, "
             "Pi/2).");
 
-    using SimTK::square;
-    const auto normFiberWidth = sin(get_pennation_angle_at_optimal());
-    m_fiberWidth = get_optimal_fiber_length() * normFiberWidth;
-    m_squareFiberWidth = square(m_fiberWidth);
-    m_maxContractionVelocityInMetersPerSecond =
-            get_max_contraction_velocity() * get_optimal_fiber_length();
-    m_kT = log((1.0 + c3) / c1) /
-           (1.0 + get_tendon_strain_at_one_norm_force() - c2);
-    m_isTendonDynamicsExplicit =
-            get_tendon_compliance_dynamics_mode() == "explicit";
+    m_isTendonDynamicsExplicit = get_tendon_compliance_dynamics_mode() == "explicit";
 }
 
 void DeGrooteFregly2016Muscle::extendAddToSystem(
@@ -279,13 +270,13 @@ void DeGrooteFregly2016Muscle::calcMuscleLengthInfoHelper(
     // ------
     mli.fiberLengthAlongTendon = muscleTendonLength - mli.tendonLength;
     mli.fiberLength = sqrt(
-            SimTK::square(mli.fiberLengthAlongTendon) + m_squareFiberWidth);
+            SimTK::square(mli.fiberLengthAlongTendon) + getSquareFiberWidth());
     mli.normFiberLength = mli.fiberLength / get_optimal_fiber_length();
 
     // Pennation.
     // ----------
     mli.cosPennationAngle = mli.fiberLengthAlongTendon / mli.fiberLength;
-    mli.sinPennationAngle = m_fiberWidth / mli.fiberLength;
+    mli.sinPennationAngle = getFiberWidth() / mli.fiberLength;
     mli.pennationAngle = asin(mli.sinPennationAngle);
 
     // Multipliers.
@@ -311,7 +302,7 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper(
         fvi.normFiberVelocity =
                 calcForceVelocityInverseCurve(fvi.fiberForceVelocityMultiplier);
         fvi.fiberVelocity = fvi.normFiberVelocity *
-                            m_maxContractionVelocityInMetersPerSecond;
+                            getMaxContractionVelocityInMetersPerSecond();
         fvi.fiberVelocityAlongTendon =
                 fvi.fiberVelocity / mli.cosPennationAngle;
         fvi.tendonVelocity =
@@ -331,13 +322,13 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper(
         fvi.fiberVelocity =
                 fvi.fiberVelocityAlongTendon * mli.cosPennationAngle;
         fvi.normFiberVelocity =
-                fvi.fiberVelocity / m_maxContractionVelocityInMetersPerSecond;
+                fvi.fiberVelocity / getMaxContractionVelocityInMetersPerSecond();
         fvi.fiberForceVelocityMultiplier =
                 calcForceVelocityMultiplier(fvi.normFiberVelocity);
     }
 
     const SimTK::Real tanPennationAngle =
-            m_fiberWidth / mli.fiberLengthAlongTendon;
+            getFiberWidth() / mli.fiberLengthAlongTendon;
     fvi.pennationAngularVelocity =
             -fvi.fiberVelocity / mli.fiberLength * tanPennationAngle;
 }

OpenSim/Actuators/DeGrooteFregly2016Muscle.h

@@ -80,10 +80,29 @@ that support implicit dynamics (i.e. Moco) and cannot be used to perform a
 time-stepping forward simulation with Manager; use explicit mode for 
 time-stepping.
 
-@note Normalized tendon force is bounded in the range [0, 5] in this class.
-   The methods getMinNormalizedTendonForce() and 
+@section property Property Bounds
+The acceptable bounds for each property are enforced at model initialization.
+These bounds are:
+ - activation_time_constant: (0, inf]
+ - deactivation_time_constant: (0, inf]
+ - active_force_width_scale: [1, inf]
+ - fiber_damping: [0, inf]
+ - passive_fiber_strain_at_one_norm_force: (0, inf]
+ - tendon_strain_at_one_norm_force: (0, inf]
+ - pennation_angle_at_optimal: [0, Pi/2)
+ - default_activation: (0, inf]
+ - default_normalized_tendon_force: [0, 5]
+
+@note The methods getMinNormalizedTendonForce() and
    getMaxNormalizedTendonForce() provide these bounds for use in custom solvers.
 
+@note Muscle properties can be optimized using MocoParameter. The acceptable
+bounds for each property are **not** enforced during parameter optimization, so
+the user must supply these bounds to MocoParameter.
+
+@note The properties `default_activation` and `default_normalized_tendon_force`
+cannot be optimized because they are applied during model initialization only.
+
 @section departures Departures from the Muscle base class
 
 The documentation for Muscle::MuscleLengthInfo states that the
@@ -107,32 +126,32 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
 public:
     OpenSim_DECLARE_PROPERTY(activation_time_constant, double,
             "Smaller value means activation can increase more rapidly. "
-            "Default: 0.015 seconds.");
+            "Default: 0.015 seconds. Bounds: (0, inf]");
     OpenSim_DECLARE_PROPERTY(deactivation_time_constant, double,
             "Smaller value means activation can decrease more rapidly. "
-            "Default: 0.060 seconds.");
+            "Default: 0.060 seconds. Bounds: (0, inf]");
     OpenSim_DECLARE_PROPERTY(default_activation, double,
             "Value of activation in the default state returned by "
-            "initSystem(). Default: 0.5.");
+            "initSystem(). Default: 0.5. Bounds: (0, inf]");
     OpenSim_DECLARE_PROPERTY(default_normalized_tendon_force, double,
             "Value of normalized tendon force in the default state returned by "
-            "initSystem(). Default: 0.5.");
+            "initSystem(). Default: 0.5. Bounds: [0, 5].");
     OpenSim_DECLARE_PROPERTY(active_force_width_scale, double,
             "Scale factor for the width of the active force-length curve. "
-            "Larger values make the curve wider. Default: 1.0.");
+            "Larger values make the curve wider. Default: 1.0. Bounds: [1, inf]");
     OpenSim_DECLARE_PROPERTY(fiber_damping, double,
             "Use this property to define the linear damping force that is "
             "added to the total muscle fiber force. It is computed by "
             "multiplying this damping parameter by the normalized fiber "
-            "velocity and the max isometric force. Default: 0.");
+            "velocity and the max isometric force. Default: 0. Bounds: [0, inf]");
     OpenSim_DECLARE_PROPERTY(ignore_passive_fiber_force, bool,
             "Make the passive fiber force 0. Default: false.");
     OpenSim_DECLARE_PROPERTY(passive_fiber_strain_at_one_norm_force, double,
             "Fiber strain when the passive fiber force is 1 normalized force. "
-            "Default: 0.6.");
+            "Default: 0.6. Bounds: (0, inf]");
     OpenSim_DECLARE_PROPERTY(tendon_strain_at_one_norm_force, double,
             "Tendon strain at a tension of 1 normalized force. "
-            "Default: 0.049.");
+            "Default: 0.049. Bounds: (0, inf]");
     OpenSim_DECLARE_PROPERTY(tendon_compliance_dynamics_mode, std::string,
             "The dynamics method used to enforce tendon compliance dynamics. "
             "Options: 'explicit' or 'implicit'. Default: 'explicit'. ");
@@ -156,6 +175,7 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
 
     DeGrooteFregly2016Muscle() { constructProperties(); }
 
+
 protected:
     //--------------------------------------------------------------------------
     // COMPONENT INTERFACE
@@ -493,14 +513,14 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
     // TODO: In explicit mode, do not allow negative tendon forces?
     SimTK::Real calcTendonForceMultiplier(
             const SimTK::Real& normTendonLength) const {
-        return c1 * exp(m_kT * (normTendonLength - c2)) - c3;
+        return c1 * exp(getTendonStiffnessParameter() * (normTendonLength - c2)) - c3;
     }
 
     /// This is the derivative of the tendon-force length curve with respect to
     /// normalized tendon length.
     SimTK::Real calcTendonForceMultiplierDerivative(
             const SimTK::Real& normTendonLength) const {
-        return c1 * m_kT * exp(m_kT * (normTendonLength - c2));
+        return c1 * getTendonStiffnessParameter() * exp(getTendonStiffnessParameter() * (normTendonLength - c2));
     }
 
     /// This is the integral of the tendon-force length curve with respect to
@@ -513,8 +533,10 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
         const double minNormTendonLength =
                 calcTendonForceLengthInverseCurve(m_minNormTendonForce);
         const double temp1 =
-                exp(m_kT * normTendonLength) - exp(m_kT * minNormTendonLength);
-        const double temp2 = c1 * exp(-c2 * m_kT) / m_kT;
+                exp(getTendonStiffnessParameter() * normTendonLength)
+                - exp(getTendonStiffnessParameter() * minNormTendonLength);
+        const double temp2 = c1 * exp(-c2 * getTendonStiffnessParameter())
+                                / getTendonStiffnessParameter();
         const double temp3 = c3 * (normTendonLength - minNormTendonLength);
         return temp1 * temp2 - temp3;
     }
@@ -523,7 +545,8 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
     /// normalized tendon length as a function of the normalized tendon force.
     SimTK::Real calcTendonForceLengthInverseCurve(
             const SimTK::Real& normTendonForce) const {
-        return log((1.0 / c1) * (normTendonForce + c3)) / m_kT + c2;
+        return log((1.0 / c1) * (normTendonForce + c3))
+                / getTendonStiffnessParameter() + c2;
     }
 
     /// This returns normalized tendon velocity given the derivative of 
@@ -534,8 +557,8 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
     SimTK::Real calcTendonForceLengthInverseCurveDerivative(
             const SimTK::Real& derivNormTendonForce,
             const SimTK::Real& normTendonLength) const {
-        return derivNormTendonForce /
-               (c1 * m_kT * exp(m_kT * (normTendonLength - c2)));
+        return derivNormTendonForce / (c1 * getTendonStiffnessParameter()
+            * exp(getTendonStiffnessParameter() * (normTendonLength - c2)));
     }
 
     /// This computes both the total fiber force and the individual components
@@ -632,8 +655,8 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
         // pennationAngle = asin(fiberWidth/fiberLength)
         // d_pennationAngle/d_fiberLength =
         //          d/d_fiberLength (asin(fiberWidth/fiberLength))
-        return (-m_fiberWidth / square(fiberLength)) /
-               sqrt(1.0 - square(m_fiberWidth / fiberLength));
+        return (-getFiberWidth() / square(fiberLength)) /
+               sqrt(1.0 - square(getFiberWidth() / fiberLength));
     }
 
     /// The derivative of the fiber force along the tendon with respect to fiber
@@ -838,6 +861,23 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
     void calcMusclePotentialEnergyInfoHelper(const bool& ignoreTendonCompliance,
             const MuscleLengthInfo& mli, MusclePotentialEnergyInfo& mpei) const;
 
+    SimTK::Real getFiberWidth() const {
+        const auto normFiberWidth = sin(get_pennation_angle_at_optimal());
+        return get_optimal_fiber_length() * normFiberWidth;
+    }
+    SimTK::Real getSquareFiberWidth() const {
+        return SimTK::square(getFiberWidth());
+    }
+    SimTK::Real getMaxContractionVelocityInMetersPerSecond() const {
+        return get_max_contraction_velocity() * get_optimal_fiber_length();
+    }
+    /// Tendon stiffness parameter kT from De Groote et al., 2016. Instead of
+    /// kT, users specify tendon strain at 1 norm force, which is more intuitive.
+    SimTK::Real getTendonStiffnessParameter() const {
+        return log((1.0 + c3) / c1) /
+           (1.0 + get_tendon_strain_at_one_norm_force() - c2);
+    }
+
     /// This is a Gaussian-like function used in the active force-length curve.
     /// A proper Gaussian function does not have the variable in the denominator
     /// of the exponent.
@@ -949,14 +989,6 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle {
 
     // Computed from properties.
     // -------------------------
-
-    // The square of (fiber_width / optimal_fiber_length).
-    SimTK::Real m_fiberWidth = SimTK::NaN;
-    SimTK::Real m_squareFiberWidth = SimTK::NaN;
-    SimTK::Real m_maxContractionVelocityInMetersPerSecond = SimTK::NaN;
-    // Tendon stiffness parameter from De Groote et al., 2016. Instead of
-    // kT, users specify tendon strain at 1 norm force, which is more intuitive.
-    SimTK::Real m_kT = SimTK::NaN;
     bool m_isTendonDynamicsExplicit = true;
 
     // Indices for MuscleDynamicsInfo::userDefinedDynamicsExtras.

OpenSim/Actuators/ModelOperators.h

@@ -19,9 +19,10 @@
  * -------------------------------------------------------------------------- */
 
 #include "ModelProcessor.h"
-
 #include <OpenSim/Actuators/ModelFactory.h>
+#include <OpenSim/Common/GCVSplineSet.h>
 #include <OpenSim/Simulation/Model/ExternalLoads.h>
+#include "OpenSim/Simulation/TableProcessor.h"
 
 namespace OpenSim {
 
@@ -208,6 +209,47 @@ class OSIMACTUATORS_API ModOpReplacePathsWithFunctionBasedPaths
     }
 };
 
+/// Prescribe motion to Coordinate%s in a model by providing a table containing
+/// time series data of Coordinate values. Any columns in the provided table
+/// (e.g., "/jointset/ankle_r/ankle_angle_r/value") that do not match a valid
+/// path to a Joint Coordinate value in the model will be ignored. A GCVSpline
+/// function is created for each column of Coordinate values and this function
+/// is assigned to the `prescribed_function` property for the matching Coordinate.
+/// In addition, the `prescribed` property for each matching Coordinate is set
+/// to "true".
+class OSIMACTUATORS_API ModOpPrescribeCoordinateValues : public ModelOperator {
+    OpenSim_DECLARE_CONCRETE_OBJECT(
+            ModOpPrescribeCoordinateValues, ModelOperator);
+    OpenSim_DECLARE_PROPERTY(coordinate_values, TableProcessor,
+            "The table of coordinate value data to prescribe to the model")
+
+public:
+    ModOpPrescribeCoordinateValues(TableProcessor table) {
+        constructProperty_coordinate_values(table);
+    }
+    void operate(Model& model, const std::string&) const override {
+        model.finalizeFromProperties();
+        TimeSeriesTable table = get_coordinate_values().process();
+        GCVSplineSet statesSpline(table);
+
+        for (const std::string& pathString: table.getColumnLabels()) {
+            ComponentPath path = ComponentPath(pathString);
+            if (path.getNumPathLevels() < 3) {
+                continue;
+            }
+            std::string jointPath = path.getParentPath().getParentPath().toString();
+            if (!model.hasComponent<Joint>(jointPath)) {
+                log_warn("Found column label '{}', but it does not match a "
+                         "joint coordinate value in the model.", pathString);
+                continue;
+            }
+            Coordinate& q = model.updComponent<Joint>(jointPath).updCoordinate();
+            q.setPrescribedFunction(statesSpline.get(pathString));
+            q.setDefaultIsPrescribed(true);
+        }
+    }
+};
+
 } // namespace OpenSim
 
 #endif // OPENSIM_MODELOPERATORS_H