NaNs when using OpenSim::GeometryPath::getPointForceDirections()

[Richardo Khonasty]

Hi All,

I'm trying to calculate a muscle Jacobian for an upper limb model and encountered an issue. In some states, certain muscles in the model would result in NaNs when calling getPointForceDirections() which in turn causes NaNs in my muscle Jacobian. Would someone be able to point me towards where the issue could potentially lie? If it makes a difference, the model I'm using is the Upper Extremity Model.

Thanks,

Richardo

[Dimitar Stanev]

Hi,

If you are trying to calculate the muscle moment arm matrix you can use the following code (OpenSim v4.0)

Code: Select all

Matrix calcMomentArm(const State& s, const Model& model) {
    auto& coordinateSet = model.getCoordinateSet();
    auto& actuatorSet = model.getActuators();
    vector<vector<double> > momentArm;
    for (int i = 0; i < actuatorSet.getSize(); i++) {
        auto pathActuator = dynamic_cast<const PathActuator*>(&actuatorSet[i]);
        if (pathActuator) {
            vector<double> column;
            for (int j = 0; j < coordinateSet.getSize(); j++) {
                column.push_back(
                    pathActuator->computeMomentArm(s, coordinateSet[j]));
            }
            momentArm.push_back(column);
        }
    }
    // convert to Matrix
    Matrix R(momentArm[0].size(), momentArm.size());
    for (unsigned int m = 0; m < momentArm.size(); m++) {
        for (unsigned int n = 0; n < momentArm[0].size(); n++) {
            R[n][m] = momentArm[m][n];
        }
    }
    return R;
}

If not can you be more specific in what exactly are trying to calculate (muscle Jacobian??). It is possible for the specific model that you are calculating this for coordinates that are constrained. Can you post which coordinates cause this problem?

[Richardo Khonasty]

Hi Dimitar,

Thanks for the reply. What I am trying to calculate is essentially a partial derivaties of MTU lengths over the unconstrained coordinates of the model. I've attached the code below where the getPointForceDirection() function is called. In the function iGenCoords is an array of indexes to the unconstrained coordinates of the model. The matrix C is a coupling matrix relating generalized speeds to generalized coordinate derivative.

To investigate whether the issue lie in the code or whether I am doing something else wrong, I've calculated the moment arm using the code that you attached where I was able to obtain the muscle moment arms. In a state where my code resulted in some NaNs, the moment arm calculated using your code also resulted in NaNs across the 20th column. I am using OpenSim 3.3 so could that make a difference? Have you encountered NaNs when using the code you provided? This may help me determine the root of the problem.

Regards,

Richardo

Code: Select all

SimTK::Matrix strength_server::calcL(OpenSim::Model& model, SimTK::State& si, SimTK::Array_<int> iGenCoords, SimTK::Matrix C)
{
    int nGenCoords = iGenCoords.size();
    int nMuscles = model.getMuscles().getSize();

    // Get the body forces equivalent of the point forces of the path
    OpenSim::Array<OpenSim::PointForceDirection *> pfds;
    SimTK::Vector generalizedForces;
    SimTK::Matrix F(si.getNU(),nMuscles);

    // Cycle through muscle, finding the generalized force each produces
    for (int m = 0; m < nMuscles; m++) {

        SimTK::State s = si;

        // Reset speeds to zero
        s.updU() = 0;
        model.getMultibodySystem().realize(s, SimTK::Stage::Velocity);
        SimTK::Vec3 v0;
        v0.setToZero();
        // Calculate the body forces due to a unit force in the muscle
        pfds.setSize(0);
        //Get the the path as PointForceDirections directions, which can be used to apply tension 
        //to bodies the points are connected to.
        model.getMuscles().get(m).getGeometryPath().getPointForceDirections(s,&pfds);	

        // Apply body forces along the geometry described by pfds due to a tension of 1N
        SimTK::Vector_<SimTK::SpatialVec> bodyForces(model.getNumBodies(), SimTK::SpatialVec(0));
        for(int i=0; i < pfds.getSize(); i++)
        {
            model.getMatterSubsystem().addInStationForce(s, SimTK::MobilizedBodyIndex(pfds[i]->body().getIndex()), ...
             pfds[i]->point(), pfds[i]->direction(), bodyForces);
        }

        // Convert body spatial forces F to equivalent mobility forces f based on
        // geometry (no dynamics required): f = ~J(q) * F.
        model.getMultibodySystem().getMatterSubsystem().multiplyBySystemJacobianTranspose(s, ...
        bodyForces, generalizedForces);

        // Store generalized forces into F[nCoords,nMuscles]
        F(0,m,s.getNU(),1) = generalizedForces;

    }

    // Now compute L[m][q]
    SimTK::Matrix L = (~F*C) * -1;
    return L;
}