The model maintains its own state

[mitkof6]

Hi,

I am having this issues that was not present in an older version of OpenSim that I use. I am porting my code for the new version v4.0.0_beta.

I maintain a model and a working model which is an identical copy of the original model. I noticed that when I try to compute quantities such as the system Jacobian the two models returned different quantities given the same state:

_model->getMatterSubsystem().calcSystemJacobian(s, J1); workingModel.getMatterSubsystem().calcSystemJacobian(s, J2);

Then I looked at the coordinates of the model and realized that they are different.

    for (int i = 0; i < _model->getCoordinateSet().getSize(); i++) {
        cout << _model->getCoordinateSet()[i].getName() <<
            _model->getCoordinateSet()[i].getValue(s) - 
            workingModel.getCoordinateSet()[i].getValue(s) << endl;
    }

I had the impression that the model do not have a state and the state is used to compute things. But in the new version for some reason this does not hold true. What can be the cause of the problem?

Best

[chrisdembia]

I think that you cannot use the state from one instance of a model with a different instance of the model. This is because the state's cache contains pointers back into the model. I'm not sure if this is new to 4.0 (https://github.com/simbody/simbody/pull/414/files#diff-f37d0f0d5ba0cad81170fb9f2cee1150) or was also true for 3.3. If it's new to 4.0, we should fix it. Otherwise, it would be great if we could document this better and give good error messages.

What happens if you use a different state with workingModel that you get from workingModel.initSystem()?

[mitkof6]

Basically, I have a generalized force controller which computes all the forces that act on the model (that is why I use getRigidBodyForces). For me to get them I have to realize to dynamics stage and to disable the ForceController so that I will not get an infinite loop. But since the state is const I can't disable any force so I have a working copy of the state. Then I figured out that I need realize the workingModel and not the actual model (_model) to a dynamics stage because I was getting an error in the finish method of ParallelExecutor::Task::finish when the integrator realizes to dynamics twice (ones in the calcTotalGeneralizedForcesand once during numerical integration), as you correctly pointed out. The following code worked well in the previous version but for the current version it is not.

Vector OpenSim::DynamicCompensator::f(const SimTK::State& s) const {
    PRE_PROCESS_CACHE(s, CACHE_F, value, Vector);

    calcTotalGeneralizedForces(s, workingModel, value);

    // these must be the same when no external forces are applied 
    // and muscles are disabled
    cout << c(s) + value << endl;
    cout << c(s) - g(s) << endl;// actual

    value = c(s) + value;// add Coriolis because they are not computed

    POST_PROCESS_CACHE(s, CACHE_F, value, Vector);
}

void calcTotalGeneralizedForces(const State& s, OpenSim::Model& model, Vector& f) {
    State workingState = model.getWorkingState(); // the following may be redundant but I tried everything
    workingState.updTime() = s.getTime();
    workingState.updQ() = s.getQ();
    workingState.updQDot() = s.getQDot();
    workingState.updQDotDot() = s.getQDotDot();
    workingState.updU() = s.getU();
    workingState.updUDot() = s.getUDot();
    workingState.updQDotDot() = s.getQDotDot();
    workingState.updZ() = s.getZ();
    workingState.updZDot() = s.getZDot();
    workingState.updY() = s.getY();

    // disable any actuators when computing the total force
    const OpenSim::Set<OpenSim::Actuator>& as = model.getActuators();
    for (int i = 0; i < as.getSize(); i++) {
        as[i].setAppliesForce(workingState, false);
    }

    // generalized forces and torques should be accessed at a Dynamics stage
    model.realizeDynamics(workingState);

    const Vector_<SpatialVec>& forces =
        model.getMultibodySystem().getRigidBodyForces(workingState, Stage::Dynamics);
    const Vector& torques = model.getMultibodySystem().getMobilityForces(workingState,
        Stage::Dynamics);

    model.getMatterSubsystem().multiplyBySystemJacobianTranspose(workingState, forces, f);

    // in our conviction we subtract their contribution
    f = -1.0 * torques - f;

    // enable any actuators when computing the total force
    /*for (int i = 0; i < as.getSize(); i++) {
        as[i].setAppliesForce(workingState, true);
    }*/
    //model.getForceSet()[0].setAppliesForce(workingState, true);
}

For example if I get the gravity forces of the mode for state s and workingState they are different. Then again if I pass the actual model (_model) to calcTotalGeneralizedForces I am getting the correct results but I get the error in the Parallel::Task. For some reason either the states are not identical or the model has an internal state. Then I tried to compare the coordinate values (given the same input state) on the two models in f and they were different as mentioned in my previous comment (similarly for the system Jacobean).

Can you think of a workaround because I have been straggling with this for a couple of days.

I have re-updated the comment.

[chrisdembia]

I can read your code in more depth later but can you try invalidating the cache in your workingState before you start using it?

[mitkof6]

Can you please tell me which method should I use to invalidate cache?

[chrisdembia]

Perhaps state.invalidateAllCacheAtOrAbove(SimTK::Stage::Instance);.

[chrisdembia]

See related issue https://github.com/opensim-org/opensim-core/issues/1455.

Also, I just opened an issue in Simbody (https://github.com/simbody/simbody/issues/567).

[chrisdembia]

As you can see, from Model::getWorkingState(), the model does have a member variable that is a State, but I don't know if it is used for any calculations.

[mitkof6]

Invalidating the cache did not produce the correct result. I tried different combinations, both on the actual state and the working state. My intuition tells me that the problem originates in the fact that the two models produce different results given the same state. Please note that I have calculated variables that do not depend on the parallel executors (e.g. system Jacobian). This was not a problem in the previous version. I will try to come up with a simple example to demonstrait the problem. If you come up with some idea please let me know.

[chrisdembia]

I have calculated variables that do not depend on the parallel executors (e.g. system Jacobian).

Ah good point.

@aseth1 and I have been chatting about this...he suggests trying to replace

State workingState = model.getWorkingState(); 

with

State workingState = s;

Also, it should not be necessary to update the cache variables UDot, UDotDot, etc.; just Y should need to be copied over.

How different are the results between the two situations? Could you show numbers?

[mitkof6]

I can't change to this because I am getting the errors in the Parrallel::Task. I have removed the unnecessary UDot. I also tried to link OpenSim v4.0.0_beta with a different version of Simbody which I use (version). But it seems that the error is present, thus I think that some new changes to OpenSim and a combination of misuse cause this strange behavior.

I also tried comparing only the gravity forces and they were different.

Total forces vs gravity + Coriolis (this should be the actual result):

~[-100.855 -43.9548 10.0187 0.383468 719.386 -0.0467892 -13.2727 0.0627176 0.0174332 -10.6339 17.2063 -0.469689 -11.6627 -34.1427 -3.75768 -2.67722 -3.03028 0.848464 0.639069 -0.0349104 -0.0783089 0.0332138 -0.129238]
~[-43.0502 9.70286 3.45835 0.383468 719.386 -0.0467892 -15.3565 -0.322575 -0.0970621 -4.71244 0.973876 -0.239469 0.0459451 21.4472 -1.32968 0.0794113 4.69688 1.06089 -0.209704 0.0385918 12.5435 4.48172 0.235493]

System Jacobian diff:

Matrix J1, J2;
    _model->getMatterSubsystem().calcSystemJacobian(s, J1);
    workingModel.getMatterSubsystem().calcSystemJacobian(s, J2);
    cout << J1 - J2 << endl;

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 -0.0623424 0.962096 0.258198 0 0.165566 -0.890993 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0.0846026 -0.259955 0.965579 0 0.144764 -0.443539 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0.994463 0.0824288 0.0314957 0 -0.975516 0.096977 0 0 0 0 0 0 0 0 0 0 0]
[0.063787 0.00391049 0.0807834 0 0 0 0.434867 0.0398596 1.47266e-05 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[-0.0745171 -0.0866956 0.0104139 0 0 0 -0.121031 0.00378832 -0.000115834 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 -0.0603645 0.0770193 0 0 0 0.0375582 -0.453288 0.00343046 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 -0.0623424 0.962096 0.258198 0.00429041 -0.165566 0 0.600344 0.0623424 0.907665 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0.0846026 -0.259955 0.965579 -0.0337469 -0.144764 0 -0.12566 -0.0846026 0.419158 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0.994463 0.0824288 0.0314957 0.999421 0.975516 0 -0.789808 -0.994463 0.0212418 0 0 0 0 0 0 0 0]
[0.502291 0.00326667 0.0282063 0 0 0 0.815073 0.0746199 0.000349271 0.372825 0 0 0 0 0 0 0 0 0 0 0 0 0]
[-0.195327 -0.072422 0.0203404 0 0 0 -0.303994 0.000383088 -0.00274725 -0.174231 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 -0.492979 0.216927 0 0 0 0.0769584 -0.869744 0.0813603 -0.00748367 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0.00429041 -0.165566 0.890993 0 0 0 -0.419456 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 -0.0337469 -0.144764 0.443539 0 0 0 0.907694 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0.999421 0.975516 -0.096977 0 0 0 0.012145 0 0 0 0 0 0 0]
[0.821862 -0.00119167 -0.0982108 0 0 0 0.815073 0.0746199 0.000349271 0.408074 0.033604 0 0 0 0 0 0 0 0 0 0 0 0]
[-0.304028 0.0264193 0.0254079 0 0 0 -0.303994 0.000383088 -0.00274725 -0.236682 -0.05974 0 0 0 0 0 0 0 0 0 0 0 0]
[0 -0.807327 0.339447 0 0 0 0.0769584 -0.869744 0.0813603 -0.00974374 -0.00316199 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 -0.0623424 0.962096 0.258198 0 0 0 -0.600344 0 0 0 -0.0216488 0 0 0 0 0 0]
[0 0 0 0 0 0 0.0846026 -0.259955 0.965579 0 0 0 0.12566 0 0 0 0.00337265 0 0 0 0 0 0]
[0 0 0 0 0 0 0.994463 0.0824288 0.0314957 0 0 0 0.789808 0 0 0 -0.99976 0 0 0 0 0 0]
[0.482656 0.00368574 0.039179 0 0 0 0.888003 0.076473 0.0177783 0.418795 0.033604 0 0 0 0 0 0.0105541 0 0 0 0 0 0]
[-0.320192 -0.0817128 0.0316806 0 0 0 -0.186327 -0.00505155 -0.00310578 -0.245512 -0.05974 0 0 0 0 0 -0.000792973 0 0 0 0 0 0]
[0 -0.467738 0.340306 0 0 0 0.0715201 -0.908512 -0.050529 -0.0100879 -0.00316199 0 0 0 0 0 -0.000231213 0 0 0 0 0 0]
[0 0 0 0 0 0 0.0623424 -0.962096 -0.258198 0 0 0 0 -0.0623424 -0.907665 0 0 -0.0618229 0 0 0 0 0]
[0 0 0 0 0 0 -0.0846026 0.259955 -0.965579 0 0 0 0 0.0846026 -0.419158 0 0 -0.195849 0 0 0 0 0]
[0 0 0 0 0 0 -0.994463 -0.0824288 -0.0314957 0 0 0 0 0.994463 -0.0212418 0 0 -0.978683 0 0 0 0 0]
[-0.0992984 -0.00423359 -0.0848126 0 0 0 -0.106697 -0.0291876 0.0699965 0 0 0 0 0 0 -0.0747269 0 0 0 0 0 0 0]
[0.136973 0.0938586 -0.0161992 0 0 0 0.130358 -0.0629033 -0.0156973 0 0 0 0 0 0 -0.0357469 0 0 0 0 0 0 0]
[0 0.0930255 -0.140767 0 0 0 -0.0177788 0.142296 -0.0925825 0 0 0 0 0 0 0.0907841 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 -0.00429041 0.165566 -0.890993 0 -0.0623424 -0.907665 0.419456 0 0 0.572382 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0.0337469 0.144764 -0.443539 0 0.0846026 -0.419158 -0.907694 0 0 0.811974 0 0 0 0]
[0 0 0 0 0 0 0 0 0 -0.999421 -0.975516 0.096977 0 0.994463 -0.0212418 -0.012145 0 0 0.114356 0 0 0 0]
[-0.39673 -0.00349742 -0.0425365 0 0 0 0 0 0 -0.439856 -0.439454 0.0423595 0 0.411698 -0.00652382 -9.63119e-05 -0.414373 0 0 0 0 0 0]
[-0.122835 0.0775378 0.00726099 0 0 0 0 0 0 -0.375022 -0.357959 0.0382893 0 0.185251 -0.00890127 1.50044e-05 -0.186697 0 0 0 0 0 0]
[0 0.401862 0.104431 0 0 0 0 0 0 -0.0107749 -0.127705 0.564308 0 0.0100492 0.45441 -0.00444777 0.00834299 0 0 0 0 0 0]
[0 0 0 0 0 0 0.0623424 -0.962096 -0.258198 0 0 0 0 -0.0623424 -0.907665 0 0.0216488 0 0 -0.587451 0 0 0]
[0 0 0 0 0 0 -0.0846026 0.259955 -0.965579 0 0 0 0 0.0846026 -0.419158 0 -0.00337265 0 0 -0.158582 0 0 0]
[0 0 0 0 0 0 -0.994463 -0.0824288 -0.0314957 0 0 0 0 0.994463 -0.0212418 0 0.99976 0 0 -0.79357 0 0 0]
[0.315155 -0.00447868 -0.12648 0 0 0 -0.504587 -0.0633516 0.0616647 0 0 0 0 0.809588 -0.0131721 -0.0748232 0.389847 -0.391077 0 0 0 0 0]
[0.323597 0.0992923 -0.0327683 0 0 0 -0.0109661 -0.0789064 -0.0213269 0 0 0 0 0.326575 -0.0159507 -0.0357319 0.143243 -0.139628 0 0 0 0 0]
[0 -0.329416 -0.307883 0 0 0 -0.0306994 0.490583 0.148309 0 0 0 0 0.0229698 0.877597 0.0863364 -0.0079585 0.0526456 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 -0.00429041 0.165566 -0.890993 0 -0.0623424 -0.907665 0.419456 0 0.0618229 0 0 0.0623424 0 0]
[0 0 0 0 0 0 0 0 0 0.0337469 0.144764 -0.443539 0 0.0846026 -0.419158 -0.907694 0 0.195849 0 0 -0.0846026 0 0]
[0 0 0 0 0 0 0 0 0 -0.999421 -0.975516 0.096977 0 0.994463 -0.0212418 -0.012145 0 0.978683 0 0 -0.994463 0 0]
[0.00375238 -0.00369723 -0.0819816 0 0 0 0 0 0 -0.500909 -0.499986 0.0519802 0 0.869346 -0.0109445 0.00594734 -0.0139714 0.0578645 0.000136121 0 0 0 0]
[0.0195541 0.0819676 -0.00539356 0 0 0 0 0 0 -0.33706 -0.319479 0.0271215 0 0.288242 -0.0227638 0.00412919 -0.0442462 -0.0366906 0.000457872 0 0 0 0]
[0 -0.00462968 -0.0192909 0 0 0 0 0 0 -0.00923096 -0.132268 0.601621 0 0.0299771 0.916851 -0.103203 0.000153273 0.00368705 -0.00393238 0 0 0 0]
[0 0 0 0 0 0 0.0623424 -0.962096 -0.258198 0 0 0 0 -0.0623424 -0.907665 0 0.0216488 0 -0.572382 0 0 -0.975385 0]
[0 0 0 0 0 0 -0.0846026 0.259955 -0.965579 0 0 0 0 0.0846026 -0.419158 0 -0.00337265 0 -0.811974 0 0 0.205997 0]
[0 0 0 0 0 0 -0.994463 -0.0824288 -0.0314957 0 0 0 0 0.994463 -0.0212418 0 0.99976 0 -0.114356 0 0 -0.0786714 0]
[0.37596 -0.00485701 -0.140154 0 0 0 -0.450599 -0.0608334 0.070432 0 0 0 0 0.815358 -0.00673234 -0.0679407 0.397038 -0.333212 0.0120883 0.0437495 0 0 0]
[0.285577 0.10768 -0.0298194 0 0 0 -0.185027 -0.100423 -0.0287465 0 0 0 0 0.462303 -0.0331442 -0.0317484 0.28102 -0.176318 -0.0242635 -0.143569 0 0 0]
[0 -0.388447 -0.267318 0 0 0 -0.0125069 0.393333 0.303901 0 0 0 0 0.0117847 0.941699 0.0263207 -0.00764943 0.0563327 0.111776 -0.00369626 0 0 0]
[0 0 0 0 0 0 0 0 0 -0.00429041 0.165566 -0.890993 0 0 0 0 -0.0216488 0 0.572382 0 0 0.975385 0]
[0 0 0 0 0 0 0 0 0 0.0337469 0.144764 -0.443539 0 0 0 0 0.00337265 0 0.811974 0 0 -0.205997 0]
[0 0 0 0 0 0 0 0 0 -0.999421 -0.975516 0.096977 0 0 0 0 -0.99976 0 0.114356 0 0 0.0786714 0]
[-0.950251 0.00088281 0.102621 0 0 0 0 0 0 -0.423474 -0.425891 0.0503279 0 0 0 0 -0.386771 0 -0.019586 0 0.078625 0 0]
[-0.207588 -0.0195719 0.0190383 0 0 0 0 0 0 -0.23056 -0.213267 0.00491957 0 0 0 0 -0.043356 0 -0.0041141 0 0.106858 0 0]
[0 0.95864 0.163928 0 0 0 0 0 0 -0.00596728 -0.103931 0.484896 0 0 0 0 0.00822885 0 0.127245 0 -0.00416183 0 0]

[chrisdembia]

Thanks. I wonder if the structure in the error can help us pinpoint the bug. Why is there no error in the upper block?

[mitkof6]

Is it possible that the state layout is accessed differently between the two models? The actual model (_model) contains a couple of ModelComponents that have a set of cache variables that are used to store some extra computations for efficiency.

for (int i = 0; i < _model->getNumCoordinates(); i++)
    {
        double a = _model->getCoordinateSet()[i].getValue(s);
        double b = workingModel.getCoordinateSet()[i].getValue(s);
        cout << _model->getCoordinateSet()[i].getName() << " " << a - b << endl;
        //cout << a << endl;
    }

Differences (the same 9 first components have zero difference for the whole simulation):

pelvis_tilt 0
pelvis_list 0
pelvis_rotation 0
pelvis_tx 0
pelvis_ty 0
pelvis_tz 0
hip_flexion_r 0
hip_adduction_r 0
hip_rotation_r 0
knee_angle_r -0.145004
ankle_angle_r 0.299732
subtalar_angle_r -0.0158161
mtp_angle_r -0.0513976
hip_flexion_l 0.568373
hip_adduction_l -0.375621
hip_rotation_l -0.023667
knee_angle_l 0
ankle_angle_l -0.00486526
subtalar_angle_l 0.224013
mtp_angle_l -0.0165329
lumbar_extension -0.258534
lumbar_bending -0.332488
lumbar_rotation 0.131809

Thanks for your time!

[chrisdembia]

Is it possible that the state layout is accessed differently between the two models?

I don't know of any way that this would be the case.

How exactly do you create workingModel?

[mitkof6]

At some point I load the model form a file. I remove a set of excluded forces which in this case is the generalized force controller (to solve the infinity loop problem in case of a force controller). I have another type of controller which is not a Force but a Controller. I have two modes one for performing torque driven simulation and one for muscle driven (by the Controller). I also copy the force set since the user may have appended external forces to the model (e.g. when I simulate a gait movement).

    auto self = const_cast<DynamicCompensator*>(this);
    // a model may be programmatically constructed
    if (_model->getInputFileName() == "Unassigned") {
        _model->print("temp.osim");
        _model->setInputFileName("temp.osim");
    }
    self->workingModel = Model(_model->getInputFileName());
    self->workingModel.updForceSet() = _model->updForceSet();

    for (auto force : excludedForces) {
        try {
            int index = self->workingModel.updForceSet().getIndex(
                force->getName(), 0);
            self->workingModel.updForceSet().remove(index);
        }
        catch (...) {
            cout << "Can't remove force: " << force->getName() << endl;
        }
    }

    self->workingModel.initSystem();

[chrisdembia]

@mitkof6 the corresponding bug in Simbody is now fixed; could you see if you still have this issue?

[mitkof6]

Yes thank you. I will let you know asap.

[mitkof6]

@chrisdembia it seems that I am getting the same error. Another thing that may be related to this bug is that in the previous version I made use of addModelComponent. When I moved to v4 I had to change to addComponentbecause I was getting an Exception on a random component (depending on the model):

Exception: Component PinJoint::ground_body1 cannot extendAddToSystem until it is up-to-date with its properties.

[chrisdembia]

@chrisdembia it seems that I am getting the same error.

Darnit.

[mitkof6]

Hi @chrisdembia.

I managed to solve this problem when I was re-implemented my task space work. I can't figure out why this was happening, it may be that I didn't use the ModelComponents correctly.

I believe that the new API is modular and permits different combinations and extensions (see use of lambdas as closures in the examples). I have managed to implement a version of Task Space Computed Muscle Control for non-linear muscles. There is a documentation (doxygen) and examples are provided. Please feel free to comment or test the code if you have time.

New API