Merge branch 'main' into moco_inverse_trajectory

.gitignore

@@ -56,6 +56,7 @@ delete_this_to_stop_optimization*.txt
 .idea/codeStyles/codeStyleConfig.xml
 .idea/codeStyles/Project.xml
 .idea/editor.xml
+.idea/copilot/**
 .cache/clangd/**
 **/.vscode/
 
@@ -71,6 +72,8 @@ delete_this_to_stop_optimization*.txt
 
 .idea/deployment.xml
 .idea/other.xml
+.idea/encodings.xml
+.idea/inspectionProfiles/Project_Default.xml
 
 # Clang-Tidy
 .clang-tidy
@@ -80,4 +83,3 @@ cmake-build-*
 
 ## File-based project format:
 *.iws
-.vscode/settings.json

Bindings/Java/Matlab/examples/Moco/example3DWalking/exampleMocoTrack.m

@@ -324,19 +324,15 @@ function muscleDrivenJointMomentTracking()
 % Ignore coordinates that are locked, prescribed, or coupled to other
 % coordinates via CoordinateCouplerConstraints (true by default).
 jointMomentTracking.setIgnoreConstrainedCoordinates(true);
-coordinateSet = model.getCoordinateSet();
-for i = 0:coordinateSet.getSize()-1
-    coordinate = coordinateSet.get(i);
-    coordName = coordinate.getName();
-    % Don't track generalized forces associated with pelvis residuals.
-    if contains(string(coordName), "pelvis")
-        jointMomentTracking.setWeightForCoordinate(coordName, 0);
-    end
-    % Encourage better tracking of the ankle joint moments.
-    if contains(string(coordName), "ankle")
-        jointMomentTracking.setWeightForCoordinate(coordName, 100);
-    end
-end
+
+% Do not track generalized forces associated with pelvis residuals.
+jointMomentTracking.setWeightForGeneralizedForcePattern(".*pelvis.*", 0);
+
+% Encourage better tracking of the ankle joint moments.
+jointMomentTracking.setWeightForGeneralizedForce("ankle_angle_r_moment", 100);
+jointMomentTracking.setWeightForGeneralizedForce("ankle_angle_l_moment", 100);
+
+% Add the joint moment tracking goal to the problem.
 problem.addGoal(jointMomentTracking);
 
 % Update the solver problem and tolerances.

Bindings/Java/tests/TestMocoSlidingMass.java

@@ -82,7 +82,7 @@ public static void testMocoSlidingMass() throws Exception {
     // =====================
     if (MocoCasADiSolver.isAvailable()) {
       MocoCasADiSolver ms = study.initCasADiSolver();
-      ms.set_num_mesh_intervals(100);
+      ms.set_num_mesh_intervals(50);
 
       // Now that we've finished setting up the tool, print it to a file.
       study.print("sliding_mass.omoco");

Bindings/OpenSimHeaders_common.h

@@ -22,6 +22,7 @@
 #include <OpenSim/Common/Event.h>
 #include <OpenSim/Common/Exception.h>
 #include <OpenSim/Common/ExperimentalSensor.h>
+#include <OpenSim/Common/ExpressionBasedFunction.h>
 #include <OpenSim/Common/FileAdapter.h>
 #include <OpenSim/Common/Function.h>
 #include <OpenSim/Common/FunctionSet.h>

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/Python/examples/Moco/example3DWalking/exampleMocoTrack.py

@@ -279,18 +279,15 @@ def muscleDrivenJointMomentTracking():
     # Ignore coordinates that are locked, prescribed, or coupled to other
     # coordinates via CoordinateCouplerConstraints (true by default).
     jointMomentTracking.setIgnoreConstrainedCoordinates(True)
-    coordinateSet = model.getCoordinateSet()
-    for i in range(coordinateSet.getSize()):
-        coordinate = coordinateSet.get(i)
-        coordName = coordinate.getName()
-        # Don't track generalized forces associated with pelvis residuals.
-        if 'pelvis' in coordName:
-            jointMomentTracking.setWeightForCoordinate(coordName, 0)
-
-        # Encourage better tracking of the ankle joint moments.
-        if 'ankle' in coordName:
-            jointMomentTracking.setWeightForCoordinate(coordName, 100)
-
+
+    # Do not track generalized forces associated with pelvis residuals.
+    jointMomentTracking.setWeightForGeneralizedForcePattern('.*pelvis.*', 0)
+
+    # Encourage better tracking of the ankle joint moments.
+    jointMomentTracking.setWeightForGeneralizedForce('ankle_angle_r_moment', 100)
+    jointMomentTracking.setWeightForGeneralizedForce('ankle_angle_l_moment', 100)
+
+    # Add the joint moment tracking goal to the problem.
     problem.addGoal(jointMomentTracking)
 
     # Update the solver tolerances.

Bindings/Python/tests/test_moco.py

@@ -35,6 +35,27 @@ def createSlidingMassModel():
 
     return model
 
+def createDoubleSlidingMassModel():
+    model = createSlidingMassModel()
+    body = osim.Body("body2", 10.0, osim.Vec3(0), osim.Inertia(0))
+    model.addComponent(body)
+
+    joint = osim.SliderJoint("slider2", model.getGround(), body)
+    coord = joint.updCoordinate(osim.SliderJoint.Coord_TranslationX)
+    coord.setName("position");
+    model.addComponent(joint);
+
+    actu = osim.CoordinateActuator()
+    actu.setCoordinate(coord)
+    actu.setName("actuator2")
+    actu.setOptimalForce(1)
+    model.addComponent(actu)
+
+    model.finalizeConnections()
+    model.initSystem()
+    return model
+
+
 class TestSwigAddtlInterface(unittest.TestCase):
     def test_bounds(self):
         model = osim.Model()
@@ -391,3 +412,33 @@ def test_changing_costs(self):
             # Change the weights of the costs.
             effort.setWeight(0.1)
             assert(study.solve().getFinalTime() < 0.8 * finalTime0)
+
+    def test_expression_based_parameter_goal(self):
+        study = osim.MocoStudy()
+        mp = study.updProblem()
+        mp.setModel(createDoubleSlidingMassModel())
+        mp.setTimeBounds(0, 1)
+        mp.setStateInfo("/slider/position/value", [-5, 5], 0, [0.2, 0.3])
+        mp.setStateInfo("/slider/position/speed", [-20, 20])
+        mp.setStateInfo("/slider2/position/value", [-5, 5], 1, [1.2, 1.3])
+        mp.setStateInfo("/slider2/position/speed", [-20, 20])
+
+        parameter = osim.MocoParameter("sphere_mass", "body", "mass",
+                                    osim.MocoBounds(0, 10))
+        mp.addParameter(parameter)
+        parameter2 = osim.MocoParameter("sphere2_mass", "body2", "mass",
+                                     osim.MocoBounds(0, 10))
+        mp.addParameter(parameter2)
+        total_weight = 7
+        mass_goal = osim.MocoExpressionBasedParameterGoal()
+        mp.addGoal(mass_goal)
+        mass_goal.setExpression(f"(p+q-{total_weight})^2")
+        mass_goal.addParameter(parameter, "p")
+        mass_goal.addParameter(parameter2, "q")
+
+        ms = study.initTropterSolver()
+        ms.set_num_mesh_intervals(25)
+        sol = study.solve()
+
+        self.assertAlmostEqual(sol.getParameter("sphere_mass") + sol.getParameter("sphere2_mass"),
+                               total_weight)

Bindings/common.i

@@ -60,6 +60,7 @@
 %include <OpenSim/Common/Sine.h>
 %include <OpenSim/Common/PolynomialFunction.h>
 %include <OpenSim/Common/MultivariatePolynomialFunction.h>
+%include <OpenSim/Common/ExpressionBasedFunction.h>
 
 %include <OpenSim/Common/SmoothSegmentedFunctionFactory.h>
 %include <OpenSim/Common/SmoothSegmentedFunction.h>

Bindings/moco.i

@@ -24,6 +24,7 @@ namespace OpenSim {
 %include <OpenSim/Moco/MocoGoal/MocoControlTrackingGoal.h>
 %include <OpenSim/Moco/MocoGoal/MocoContactTrackingGoal.h>
 %include <OpenSim/Moco/MocoGoal/MocoContactImpulseTrackingGoal.h>
+%include <OpenSim/Moco/MocoGoal/MocoExpressionBasedParameterGoal.h>
 %include <OpenSim/Moco/MocoGoal/MocoInitialActivationGoal.h>
 %include <OpenSim/Moco/MocoGoal/MocoJointReactionGoal.h>
 %include <OpenSim/Moco/MocoGoal/MocoSumSquaredStateGoal.h>
@@ -48,6 +49,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

@@ -8,6 +8,20 @@ This is not a comprehensive list of changes but rather a hand-curated collection
 
 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)
+- `PointForceDirection` no longer has a virtual destructor, is `final`, and its `scale` functionality
+  has been marked as `[[deprecated]]` (#3890)
+- Added `ExpressionBasedFunction` for creating `Function`s based on user-defined mathematical expressions. (#3892)
+
+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).
@@ -68,10 +82,7 @@ pointer to avoid crashes in scripting due to invalid pointer ownership (#3781).
 - Fixed bugs in `MocoCasOCProblem` and `CasOC::Problem` with incorrect string formatting. (#3828)
 - Fixed `MocoOrientationTrackingGoal::initializeOnModelImpl` to check for missing kinematic states, but allow other missing columns. (#3830)
 - Improved exception handling for internal errors in `MocoCasADiSolver`. Problems will now abort and print a descriptive error message (rather than fail due to an empty trajectory). (#3834)
-- 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)
-
+- Upgraded the Ipopt dependency Metis to version 5.1.0 on Unix and macOS to enable building on `osx-arm64` (#3874).
 
 v4.5
 ====

CHANGELOG_MOCO.md

@@ -1,17 +1,52 @@
 Moco Change Log
 ===============
 
-1.3.1
+1.4.0
 -----
 - 2024-08-30: Added `MocoInverse::initializeKinematics()` to allow users to retrieve kinematics after
               converting the `MocoInverse` to a `MocoStudy`. 
 
+- 2024-08-26: Changed all `printDescription()` and `printDescriptionImpl()` methods to log at level info 
+              instead of cout.
+
+- 2024-08-21: Added support for the 'projection' method for enforcing kinematic
+              constraints from Bordalba et al. (2023) to `MocoCasADiSolver`. 
+              This method allows enforcing kinematic constraints with any
+              transcription scheme and can be enabled using by setting the
+              property `kinematic_constraint_method` to `'Bordalba2023'`. The 
+              existing default method from Posa et al. (2016) can also be
+              specified using `'Posa2016'`.
+
+- 2024-08-21: The `MocoSolver` properties `interpolate_control_midpoints` and
+              `enforce_path_constraint_midpoints` have been renamed to 
+              `interpolate_control_mesh_interior_points` and 
+              `enforce_path_constraint_mesh_interior_points`, respectively. This
+              change was made to clarify that these properties apply to the new
+              Legendre-Gauss and Legendre-Gauss-Radau transcription schemes which 
+              have collocation points in the mesh interval interior that do not 
+              necessarily coincide with the mesh interval midpoints (unlike 
+              Hermite-Simpson transcription).
+
+- 2024-08-15: Added `MocoExpressionBasedParameterGoal` to enable minimizing any arithmetic expression
+              of parameter values.
+
+- 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`
               component would not properly set the muscle masses listed, resulting in incorrect
               metabolics values being computed.
@@ -22,9 +57,9 @@ Moco Change Log
               (e.g., getInputControlsTrajectory(), generateControlsFromModelControllers()).
               `MocoGoal`s and `MocoPathConstraint`s have been updated to support
               Input controls. See the Moco User Guide for details.
-
+              
 - 2024-04-01: Added `MocoGeneralizedForceTrackingGoal` to enable joint moment tracking
-              in `MocoProblem`s, and added the utility `calcGeneralizeForces()` to 
+              in `MocoProblem`s, and added the utility `calcGeneralizedForces()` to 
               `MocoStudy` for computing joint moments from a `MocoTrajectory`. 
               Added a sub-example to exampleMocoTrack (C++, Python, and Matlab) to 
               feature this new functionality.

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