Add NMSM Pipeline contact elements for use in Moco

[SpencerTWilliams]

Brief summary of changes

This fork includes the new MeyerFregly2016Muscle and completes the implementation of the MeyerFregly2016Force included in the StationPlaneContactForce class.

The MeyerFregly2016Muscle is intended to be used without tendon compliance.

Testing I've completed

I've successfully compiled OpenSim on Ubuntu with these changes. However, I have not been able to test using the new elements.

Looking for feedback on...

I ran into an error while trying to compile a test for the MeyerFregly2016Muscle by modifying testMocoActuators.cpp:

/home/compiler/opensim-workspace/opensim-core-source/OpenSim/Moco/Test/testMocoActuators.cpp:95:26: error: expected type-specifier before ‘MeyerFregly2016Muscle’
   95 |         auto* actu = new MeyerFregly2016Muscle();
      |                          ^~~~~~~~~~~~~~~~~~~~~

I had the same error if I replaced the auto with MeyerFregly2016Muscle, so I think there may be an issue with including the new muscle as a type, even though the code files for the muscle compile correctly. Are there any issues with how I included the muscle in linking or registering the type in 2063d6f?

CHANGELOG.md (choose one)

• updated, but needs a merge conflict resolved I don't have permission for.

---

This change is 

[nickbianco]

Thanks for creating the PR @SpencerTWilliams! I've left a handful of comments to get started. My biggest question right now is whether or not we should create a separate class the muscle or if it is possible to use DeGrooteFregly2016Muscle, since the implementations look very similar. Could you list the primary differences between the two classes? I'd like to get a better understanding before I review.

I ran into an error while trying to compile a test for the MeyerFregly2016Muscle by modifying testMocoActuators.cpp:

I think you might have forgotten an include statement in the test file? Feel free to push the test changes if you want me to take a look.

Finally, we use Reviewable for code reviews (click the purple button in the PR description). Let me know if you have any questions on how it works.

Reviewed 6 of 8 files at r1, all commit messages.
Reviewable status: 6 of 8 files reviewed, 10 unresolved discussions (waiting on @SpencerTWilliams)

---

CHANGELOG.md line 63 at r1 (raw file):

- Fixed bug in `report.py` preventing plotting multiple MocoParameter values. (#3808)
- Added SynergyController, a controller that computes controls for a model based on a linear combination of a set of Input control signals and a set of synergy vectors. (#3796)
- Added `MeyerFregly2016Muscle` and completed the implementation of the `MeyerFregly2016Force` included in the `StationPlaneContactForce` class to support NMSM Pipeline-equivalent muscle and contant models in Moco. 

Typically we append the PR number in changelog entries:

Suggestion:

- Added `MeyerFregly2016Muscle` and completed the implementation of the `MeyerFregly2016Force` included in the `StationPlaneContactForce` class to support NMSM Pipeline-equivalent muscle and contant models in Moco. (#3877)

---

OpenSim/Actuators/MeyerFregly2016Muscle.h line 41 at r1 (raw file):

This muscle implementation is based on the previously implemented 
DeGrooteFregly2016Muscle.

At first glance, this implementation looks very similar to DeGrooteFregly2016Muscle, which makes me wonder again if it's worth making these muscle separate classes. I think it would be worth enumerating the differences between the curves in DeGrooteFregly2016Muscle and the ones here to get a bette sense if an entirely separate class is necessary.

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 6 at r1 (raw file):

 * OpenSim: StationPlaneContactForce.h                                        *
 * -------------------------------------------------------------------------- *
 * Copyright (c) 2017 Stanford University and the Authors                     *

Might as well update the year too:

Suggestion:

 * Copyright (c) 2024 Stanford University and the Authors                     *

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 146 at r1 (raw file):

B. J. (2016). Muscle Synergies Facilitate Computational Prediction of
Subject-Specific Walking Motions. Frontiers in Bioengineering and
Biotechnology, 4, 1055–27. http://doi.org/10.3389/fbioe.2016.00077 

Now that this contact element is being fleshed out, it would be good to expand the class description here (and possibly update the reference, if relevant). Consider explaining how the contact forces are compute and include equations. Refer to MocoControlGoal on how to write equations using Doxygen comments.

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 178 at r1 (raw file):

        const SimTK::Real y = pos[1];
        const SimTK::Real velNormal = vel[1];
        SimTK::Real velSliding = sqrt(vel[0] * vel[0] + vel[2] * vel[2]);

The OpenSim convention is that y-direction is vertical, but it would be good to note this somewhere in the class documentation above.

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 185 at r1 (raw file):

        const SimTK::Real h = 1e-3;
        const SimTK::Real c = 5e-4;
        const SimTK::Real ymax = 1e-2;

It would be good to explain these hardcoded parameters in the class docs.

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 198 at r1 (raw file):

        if (SimTK::isNaN(Fspring) || SimTK::isInf(Fspring)) {
            Fspring = 0;
        }

I'm wondering if this if-statment should be removed. I'm not sure if we added these originally or borrowed them from the Meyer et al. implementation.

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 211 at r1 (raw file):

        if (velSliding < 1e-10) {
                velSliding = 0;
        }

Why not just let this be near zero? We try to avoid if-statements like this, as they can introduce discontinuties in the model.

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 217 at r1 (raw file):

        if (SimTK::isNaN(horizontalForce) || SimTK::isInf(horizontalForce)) {
            horizontalForce = 0;
        }

Similar to my comment above: it would be ideal to not require this if-statement. I'm also wary about hiding NaN or Inf values that might indicate deeper issues with the model.

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 219 at r1 (raw file):

        }
        
        const SimTK::Real slipOffset = 1e-4;

What is this? Should it be a property or grouped with the other hardcoded parameters above?

[nickbianco]

@SpencerTWilliams, just checking in on this PR. Let me know if you have any questions about the review.

[SpencerTWilliams]

Thanks for the suggestions! We found the documentation for our contact force equation and the meanings of the constants, so I'll push my changes with the updated class description.

The MeyerFregly2016Muscle class is based on DeGrooteFregly2016Muscle and is very similar, but it uses different active and passive force-length and force-velocity curves with different numbers of coefficients. I made this class since we want a muscle that is completely consistent with the muscle model in our code, though it will look similar to the already existing DeGrooteFregly2016Muscle.

For your comments on the if statements, I included those to guarantee that the contact force elements work exactly the same as they do in our code. We haven't had issues with discontinuities, but if these statements will cause problems for Moco, we can try removing them. The slipOffset constant likely accounts for any of these potential issues already as it prevents dividing by zero when calculating horizontal forces. This constant is not a property, so I'll move it up to be with the other hardcoded values for clarity.

[nickbianco]

Thanks for the update @SpencerTWilliams!

The MeyerFregly2016Muscle class is based on DeGrooteFregly2016Muscle and is very similar, but it uses different active and passive force-length and force-velocity curves with different numbers of coefficients.

Could you be more specific about these differences? My understanding is that the curves in the MeyerFregly2016Muscle entirely different, or do they use the same general formulations (e.g., sum of Gaussians for the force-length curve) with slightly different coefficients? If it's the latter, then I'm somewhat inclined to not introduce a new muscle class.

For your comments on the if statements, I included those to guarantee that the contact force elements work exactly the same as they do in our code. We haven't had issues with discontinuities, but if these statements will cause problems for Moco, we can try removing them.

I'm generally concerned about hiding Inf or NaN values in a component accessible to the wider userbase. If your force is produce Inf or NaN values something should probably be changed about your model, initial guess, etc. Zeroing out velSliding in that if-statement will introduce a discontinuity, but that doesn't mean your probably won't converge, especially with finite differences. I would want to know if these if-statements are merely safeguards or if they provide some functional benefit before agreeing to leave them in.

[nickbianco]

Also, just a heads up that the Force API will be changing based on #3891 soon. It might make sense to break up these two changes: we can use this PR for the contact force addition and then (if needed) open a second PR for the new muscle class.

[SpencerTWilliams]

I looked at the curve differences again, and the active force-length curves do have a similar formulation. However, the passive force-length and force-velocity curves are very different between the two muscles, with different numbers of coefficients. I think these differences would require a new muscle class.

I'll test our model more to see whether it's possible to remove the Inf and NaN values. I think this check would be safe to remove at this point, so I'll remove it if it isn't possible to get Inf and NaN values from actual double inputs.

Let me know if it ends up making more sense to make a different PR for the muscle. Would I need to change the contact force class to account for the Force API changes?

[adamkewley]

The force API changes are non-breaking so, in general, code can still override+use the "raw" computeForce APIs. The ForceProducer + ForceConsumer APIs are extra - they're to facilitate other systems, such as visualizers, force reporters, and so on.

Specifically for this PR, I don't think you need to change what you're doing in order to receive the benefits of the new API. I've already ensured that the new API is rolled out to OpenSim::PathActuator, which OpenSim::Muscle inherits from:

• https://github.com/opensim-org/opensim-core/pull/3891/files#diff-5531d7f57434db355a0f850f1439e04caeab955f43dc90a8b5ee22fb740c4ef5R201

Which ultimately means that any forces/actuations produced by MeyerFregly2016Muscle muscle class will automatically be report-able, visualize-able, etc.

Your changes to the StationPlaneContactForce look like they're separate to the API changes I made in order to ensure StationPlaneContactForce emits its forces as point forces to a ForceConsumer (so that they could potentially be visualized). The change I made is here:

• https://github.com/opensim-org/opensim-core/pull/3891/files#diff-70c5e57686d093bb16e1765aeb5cb014cb317a50baaef60508ba2c8c83812faaR76

But it's far away enough from your modifications that you may not even see a merge error.

[nickbianco]

@adamkewley, thanks for the info!

@SpencerTWilliams, yes, let's focus on the contact force additions here and have a separate PR for the muscle.

[SpencerTWilliams]

I've pushed my documentation and changes based on your suggestions. The NaN and Inf checks have been deleted from the force element. Should I go ahead and make the second PR for the muscle?

[nickbianco]

The updates look good overall! I just have one comment to address since I'm still unsure about the if-statement for velSliding.

Should I go ahead and make the second PR for the muscle?

Yes, go for it!

Reviewed 2 of 2 files at r4, all commit messages.
Reviewable status: 6 of 8 files reviewed, 5 unresolved discussions (waiting on @SpencerTWilliams)

---

CHANGELOG.md line 26 at r4 (raw file):

  means that API users can now now ask many `Force` components what forces they produce (see #3891 for a
  comprehensive overview).
- Added `MeyerFregly2016Muscle` and completed the implementation of the `MeyerFregly2016Force` included in the `StationPlaneContactForce` class to support NMSM Pipeline-equivalent muscle and contant models in Moco. (#3877)

Update the changelog when removing MeyerFregly2016Muscle from this PR.

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 211 at r1 (raw file):

Previously, nickbianco (Nick Bianco) wrote…

Why not just let this be near zero? We try to avoid if-statements like this, as they can introduce discontinuties in the model.

I'm still unsure about keeping this if-statement. What would be the downside of removing it?

[SpencerTWilliams]

It looks like it's okay to remove that if-statement, so I took it out. I'll make the new muscle PR now.

[nickbianco]

@SpencerTWilliams, I've left a few minor formatting suggestions. Now that the component looks complete and the muscle is no longer in the PR I have a few additional thoughts:

1. It looks like there are no testing changes currently in the PR. It looks like the contact force is not tested at all in testMocoContact. At the very least, the contact model should be added to the test there. It might make sense to beef up the tests for the StationPlaneContactForce elements, we have not looked at them in a while.
2. Now that it is being formalized, we should probably move StationPlaneContactForce.h/.cpp next to SmoothSphereHalfSpaceForce under OpenSim/Simulation/Model. The OpenSim/Moco/Components folder is mostly reserved for internal classes. (Note that OSIMMOCO_API should be changed to OSIMSIMULATION_API when moved.) I would also move AckermannVanDenBogert2010Force to the bottom of the file so that MeyerFregly2016Force is the first concrete class under StationPlaneContactForce.
3. As a follow on from 2., now that MeyerFregly2016Force is becoming a proper OpenSim force element, we should probably also spruce up the base class StationPlaneContactForce. This would involve adding some class documentation and addressing the TODO comments regarding caching. I can address these TODOs in a later PR, but if you want to take a stab at the class documentation that would be great.

Reviewed 8 of 8 files at r5, all commit messages.
Reviewable status: all files reviewed, 3 unresolved discussions (waiting on @SpencerTWilliams)

---

CHANGELOG.md line 26 at r5 (raw file):

  means that API users can now now ask many `Force` components what forces they produce (see #3891 for a
  comprehensive overview).
- Completed the implementation of the `MeyerFregly2016Force` included in the `StationPlaneContactForce` class to support NMSM Pipeline-equivalent contant models in Moco. (#3877)

Suggestion:

- Completed the implementation of the `MeyerFregly2016Force` included in the `StationPlaneContactForce` class to support NMSM Pipeline-equivalent contact models in Moco. (#3877)

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 153 at r5 (raw file):

Biotechnology, 4, 1055–27. http://doi.org/10.3389/fbioe.2016.00077 

Following OpenSim convention, this contact element assumes that the y direction 

Suggestion:

Following OpenSim convention, this contact element assumes that the y-direction

---

OpenSim/Moco/Components/StationPlaneContactForce.h line 197 at r5 (raw file):

with a tanh function) and viscous (modeled with a linear function) friction 
models may be used. 
 */

Suggestion:

models may be used. */

[SpencerTWilliams]

I made your suggested change to the documentation and added to the contact element test. The MeyerFregly2016Force is now uncommented, and I commented out the other concrete types for now since they aren't fully implemented.

I wrote a test for known kinematics based on how the others were formatted. I'm haven't been able to run these tests with the code changes, but the values I gave are what I would expect the model to produce. Does the test I added make sense?

[nickbianco]

The added test looks good! I added a few comments.

Reviewed 2 of 2 files at r7, all commit messages.
Reviewable status: all files reviewed, 3 unresolved discussions (waiting on @SpencerTWilliams)

---

OpenSim/Moco/Test/testMocoContact.cpp line 261 at r7 (raw file):

// set of input kinematics and default parameters. 
template<typename T>
void testKnownKinematics() {

This test looks good, but you'll need to wrap it in a TEST_CASE to enable it. Would you expect it to pass with the other models? If so, then you could add it testStationPlaneContactForce<T>. Otherwise, you could create a standalone test just for MeyerFregly2016Force.

---

OpenSim/Moco/Test/testMocoContact.cpp line 541 at r7 (raw file):

TEMPLATE_TEST_CASE("testStationPlaneContactForce", "[tropter]", 
        /*AckermannVanDenBogert2010Force, EspositoMiller2018Force,*/

The tests are still passing with the other contact models, then I would leave those models uncommented here.

---

OpenSim/Simulation/Model/StationPlaneContactForce.h line 32 at r7 (raw file):

Vertical force is calculated using stiffness and dissipation coefficients, 
while horizontal force is calculated using a friction model. Multiple concrete 
implementation of this class are available, but currently MeyerFregly2016Force 

Suggestion:

implementations of this class are available, but currently MeyerFregly2016Force

[SpencerTWilliams]

I added a TEST_CASE based on how the others were set up. These calculations are specific to the type of model, so I made the test specific to the MeyerFregly2016Force.

[nickbianco]

I just enabled the CI given the recent test changes, let's see if they pass! You should also be able to run the test suite locally.

Reviewed 2 of 2 files at r8, all commit messages.
Reviewable status: complete! all files reviewed, all discussions resolved (waiting on @SpencerTWilliams)

[nickbianco]

@SpencerTWilliams, the configuration step in these steps are failing because the CMakeLists.txt file in OpenSim/Moco is looking for StationPlaneContactForce, but it is no longer there. The lines for those components need to be removed.

The CI failures on Windows and Mac are unrelated, a fix is in the works for those (#3921).

[nickbianco]

Reviewed 2 of 3 files at r9, all commit messages.
Reviewable status: 10 of 11 files reviewed, 1 unresolved discussion (waiting on @SpencerTWilliams)

---

OpenSim/Moco/Test/testMocoContact.cpp line 646 at r9 (raw file):

TEST_CASE("testMeyerFregly2016ForceValues", "[casadi]") {
    testKnownKinematics();

Since testKnownKinematics() only tests the MeyerFregly2016Force component, I think you could just move that code directly in to this TEST_CASE.

[nickbianco]

@SpencerTWilliams I've resolved the merge conflicts and left another minor comment. Let's see if the tests pass!

Reviewed 1 of 3 files at r9.
Reviewable status: all files reviewed, 1 unresolved discussion (waiting on @SpencerTWilliams)

[nickbianco]

Reviewed 1 of 1 files at r10, all commit messages.
Reviewable status: all files reviewed, 3 unresolved discussions (waiting on @SpencerTWilliams)

---

OpenSim/Simulation/Model/StationPlaneContactForce.h line 180 at r10 (raw file):

        /// Normal force.
        y = y - resting_length;

y was marked as const above, so you cannot update it like this. You either need to remove the const qualifier above or assign this value to a new variable. (This is causing the CI to fail)

---

OpenSim/Simulation/Model/StationPlaneContactForce.cpp line 39 at r10 (raw file):

        const auto& pt = getConnectee<Station>("station");
        const auto pt1 = pt.getLocationInGround(s);
        const SimTK::Vec3 force = calcContactForceOnStation(s);

CI is failing here because you need to update to the new method name:

Suggestion:

        const SimTK::Vec3 force = computeContactForceOnStation(s);

[nickbianco]

@SpencerTWilliams, I left a few comments on why the CI builds are failing. These are compile time errors, so I would recommend recompiling locally to catch these and any other compile errors still lingering. I would also recommend running the tests locally so you don't have to wait on CI to see if they pass or fail.

[SpencerTWilliams]

Thanks for the advice on handling the tests. For running tests locally, I've been having trouble modifying the build script to compile my fork with tests on Linux, which I had thought would be more convenient at first.

I remember you had mentioned that you used CMake when you ran tests locally, so do you have directions I can follow for using CMake to compile and run OpenSim tests on Windows?

[nickbianco]

The C++ tests will be built by default, I don't think we have any CMake options to turn them off. I usually run tests locally through VS Code, but you can always run the compiled testMocoContact from the terminal.

[SpencerTWilliams]

Are there any changes I might need to make to compile with CMake locally? When I try to build through VS Code, I get this error:

[cmake] CMake Error at CMakeLists.txt:713 (find_package):
[cmake]   Could not find a package configuration file provided by "spdlog" with any
[cmake]   of the following names:
[cmake] 
[cmake]     spdlogConfig.cmake
[cmake]     spdlog-config.cmake
[cmake] 
[cmake]   Add the installation prefix of "spdlog" to CMAKE_PREFIX_PATH or set
[cmake]   "spdlog_DIR" to a directory containing one of the above files.  If "spdlog"
[cmake]   provides a separate development package or SDK, be sure it has been
[cmake]   installed.

I haven't changed the directory structure from the main branch of OpenSim, so I was wondering if it's a machine-specific issue.

[nickbianco]

spdlog is a dependency. You will need to build the CMake project located in the dependencies folder first.

But you have a working build in your Linux virtual machine? If so, the dependencies should be built there already and you should be able to run the tests.

[SpencerTWilliams]

Thanks, that error was fixed by building dependencies first! I did build OpenSim on Linux previously, but that was from the .sh build script, which clones OpenSim from the main branch before building. I also don't think it includes tests, so I wasn't sure how I needed to modify that script.

With CMake, I"m still getting a similar error:

[cmake]   Could not find a package configuration file provided by "casadi" (requested
[cmake]   version 3.4.4) with any of the following names:
[cmake] 
[cmake]     casadiConfig.cmake
[cmake]     casadi-config.cmake
[cmake] 
[cmake]   Add the installation prefix of "casadi" to CMAKE_PREFIX_PATH or set
[cmake]   "casadi_DIR" to a directory containing one of the above files.  If "casadi"
[cmake]   provides a separate development package or SDK, be sure it has been
[cmake]   installed.

It looks like it would be the same sort of issue of needing to find a directory to build first, but I can't tell where this is. Do you have a guide for which directories I should build from before the project will work?

[nickbianco]

To build CasADi, you'll need to enable the OPENSIM_WITH_CASADI CMake variable. I'd try to follow the CMake settings in the build scripts for other settings.

[SpencerTWilliams]

From the build script, I can see that CMake is only used to build dependencies and then at the opensim-core, level, so this should be the last use of CMake, then. For the OPENSIM_WITH_CASADI variable, I have this on lines 192-193 in CMakeLists.txt:

option(OPENSIM_WITH_CASADI
        "Build CasADi support for Moco (MocoCasADiSolver)." ON)

If OPENSIM_WITH_CASADI is already enabled, is there something else CMake could be missing at this point, or something with CasADi I need to install before this will work?

[nickbianco]

Sorry, you also need to set that flag when building the dependencies.

[SpencerTWilliams]

Thanks, that was turned off in the dependencies. I built the dependencies again with that flag enabled, but I'm getting the same error at the core level. In the dependencies, am I supposed to see a directory for CasADi? I'm only seeing catch2, docopt, simbody, and spdlog.

[nickbianco]

Hmm, you might need to delete your CMake cache and reload for the changes to take effect.

[SpencerTWilliams]

I'm not getting a CMakeCache.txt file, so I don't think it's a caching issue. I met with Dr. Fregly, and he'll email you about our next steps for finishing this in a minute.

[nickbianco]

Closing. Work continuing on #3949.