CreatePeriodicTrajectory in 3D - Why negate and shift?

[Pasha van Bijlert]

Hi all,

I don't understand the default behaviour of CreatePeriodicTrajectory, and I'm wondering if I'm overlooking something now that I've made the switch to full 3D predictive simulations.

I'm simulating half strides with a 6DOF pelvis. To get bilaterally symmetric periodic locomotion, I think I want to enforce a constraint so that pelvis list (R_x), pelvis rotation R_y) and pelvis tz values and speeds at the end of the half stride are negative wrt to the start of the stride. In other words, I set them all using negatedStatePairs in MocoPeriodicityGoal.

According to the documentation, if I do:

fullStride = opensimMoco.createPeriodicTrajectory(gaitPredictionSolution)

pelvis list & pelvis Tz values get fed to "negateAndShift" patterns. This leads to the graphs below for pelvis values and speeds, where both pelvis_list_value and pelvis_tz_value jump at 50% of the gait cycle. I don't understand in what situation "negateAndShift" is the desired behaviour. Am I setting my constraints incorrectly, thereby leading to solutions where "negateAndShift" doesn't work, or should I manually set all the patterns to "negate", because "negateAndShift" is useful in some other situation I'm overlooking?

Cheers,
Pasha

graph.PNG (206 KiB) Viewed 575 times

[Ross Miller]

Hi Pasha,

Your setup on the half-stride symmetry constraint sounds correct to me.

I've had trouble getting half-stride / one-step simulations of walking to converge with this same setup. I think the issue is the pelvis_tz coordinate. For example, suppose the initial value of pelvis_tz at the start of the step is 100 km. Obviously we don't want the final value at the end of the step to be -100 km.

Negation seems to be assuming the coordinate will necessarily oscillate around zero. I would guess the shifting is accounting for this but it seems like it isn't working as intended in your problem.

Ross

[Nicholas Bianco]

Hi Pasha,

The default arguments for createPeriodicTrajectory() may not necessarily apply to your solutions depending on how you set up the periodicity constraints. Seems like here you can just move those coordinates to the "negate" argument to get the trajectory you're looking for.

Best,
Nick

[Pasha van Bijlert]

Hi Ross & Nick,

Thanks for debugging this with me. As I understand it, my initial hunch was correct and I need to feed all the pelvis coordinates into "negatePatterns", including the ones that by default go into "negateAndShift". Does that mean that it is impossible to get a periodic 3D trajectory in Matlab? Last year when dealing with something similar, it turned out that you cannot define a standard vector string pair in Matlab, so you cannot manually set the SymmetryPatterns.

About negateAndShift: am I correct in saying that the shifting is for a periodic coordinate that oscillates about a non-zero value? So if your pelvis list is 0 in neutral position, that implies that shifting is only appropriate for asymmetric gaits, correct? I still don't understand when "negateAndShift" is used. Is this for when you're doing a tracking problem, and your subject was slightly asymmetric during a trial?

Cheers,
Pasha

[Nicholas Bianco]

Hi Pasha,

If it's helpful, here's the code for creating a periodic half gait cycle from the tracking simulation in the Moco paper: code here.

If you specify your symmetry constraints in the same way as this example, then the default arguments for createPeriodicTrajectory work as expected. Lower down in that example you can see the arguments specified manually in Python, but these arguments are the same as the C++ defaults.

As for the issue related to custom arguments in Matlab, we'll create an issue in OpenSim internally and hopefully get a fix in near future.

Best,
Nick

[Pasha van Bijlert]

Hi Nick,

Thanks for the code link. I'm not familiar with C++ or python, but the code is understandable anyway: I think the only difference is that in that example, pelvis_list & pelvis_tz are added as regular StatePairs, while the speeds are added as a negatedStatePairs. I had set them all as negatedStatePairs. I will experiment with setting the values as regular StatePairs, but I think this could result in asymmetric gaits (pelvis list does not necessarily have to oscillate around 0, even though it will be periodic).

Thank you, also, for the possible future fix. After thinking about it for a bit, I think I can recreate createPeriodicTrajectory in Matlab by extracting each state, control & derivative as a Mat (eg "traj.getStateMat('state_name') and then performing the necessary modifications to mirror the half stride, and inserting them into a resampled trajectory with the desired time vector.
EDIT: after thinking about it further, I probably only have to do these manual computations for the "incorrect" pelvis values (so list and Tz), and can then set these states manually in the "incorrect" periodic trajectory acquired using the default settings.


Cheers,
Pasha

[Pasha van Bijlert]

This code works if you set your periodicity constraints in the way I described in the first post.

Code: Select all

fullStride = opensimMoco.createPeriodicTrajectory(gaitPredictionSolution);

%correct pelvis list and pelvis tz
pel_list_HS = gaitPredictionSolution.getStateMat('/jointset/groundPelvis/pelvis_list/value');
pel_list_FS = [pel_list_HS; -pel_list_HS(2:end)];
fullStride.setState('/jointset/groundPelvis/pelvis_list/value',pel_list_FS)


pel_tz_HS = gaitPredictionSolution.getStateMat('/jointset/groundPelvis/pelvis_tz/value');
pel_tz_FS = [pel_tz_HS; -pel_tz_HS(2:end)];
fullStride.setState('/jointset/groundPelvis/pelvis_tz/value',pel_tz_FS)

fullStride.write([filename '_periodic.sto']);

You might need to rename the joint coordinates if they're not the same as what I used.

Cheers,
Pasha

[Nicholas Bianco]

Hi Pasha, glad that works!

[Pasha van Bijlert]

Hi Nick & Ross,

Thanks! The solution was quite simple but for some reason I didn't think of it until talking it over with both of you.

Cheers,
Pasha

[Pagnon David]

EDIT: I had initially forgotten to include controls in the output file, this is now done.

Just in case this can be helpful for anyone, I had a somewhat similar issue. I wanted to output several strides from a single full stride, and createPeriodicTrajectory's default behabior is used to create a full stride from a half one. I could not find a way to NOT switch right and left sides in Python, so instead of using this I handled it by hand. There may be a better way but this works.

Code: Select all

def multiple_strides(fullStride, stride_number):
    '''
    Create multiple strides from a full stride 
    
    INPUTS:
    - fullStride: MocoTrajectory
    - stride_number: integer
    
    OUTPUT:
    - fullStride_out: MocoTrajectory
    '''
    
    # Retrieve states
    states_numtimes = fullStride.getNumTimes()
    states_times = fullStride.getTime().to_numpy()
    states_names = fullStride.getStateNames()
    states_array = fullStride.getStatesTrajectoryMat()
    controls_names = fullStride.getControlNames()
    controls_array = fullStride.getControlsTrajectoryMat()

    # Multiplie strides 
    states_numtimes_out = states_numtimes * stride_number

    time_shift = states_times[-1] + (states_times[-1]-states_times[0])/(states_numtimes-1)
    states_times_out = np.ravel([states_times + time_shift * n for n in range(stride_number)])

    pelvis_tx_ind = int(np.where(['pelvis_tx/value' in s for s in states_names])[0])
    pelvis_tx_arr = states_array[:,pelvis_tx_ind]
    pelvis_tx_out = np.ravel([pelvis_tx_arr + (pelvis_tx_arr[-1]-pelvis_tx_arr[0])*n for n in range(stride_number)])
    states_array_out = np.tile(states_array,(stride_number,1))
    states_array_out[:,pelvis_tx_ind] = pelvis_tx_out

    controls_array_out = np.tile(controls_array,(stride_number,1))

    # Create trajectory
    fullStride_out = fullStride.clone()
    fullStride_out.setNumTimes(states_numtimes_out)
    fullStride_out.setTime(states_times_out)
    [fullStride_out.setState(states_names[i], states_array_out[:,i]) for i in range(len(states_names))]
    [fullStride_out.setControl(controls_names[i], controls_array_out[:,i]) for i in range(len(controls_names))]

    return fullStride_out

This can be used this way:

Code: Select all

stride_number = 4
fullStride_out = multiple_strides(fullStride, stride_number) # make it 4 strides

OFF-TOPIC but why not share this too:
Then if you want to save it and give it a name in the GUI (instead of an empty string):

Code: Select all

def add_name_to_sto(sto_path, name):
    '''
    Prepend name to sto_path 
    or replace it if already exists
    '''
    with open(sto_path, 'r+') as file: 
        content_without_name = file.read() 
        firstline = file.readline()
        if '=' not in firstline: # if name already exists
            content_without_name = content_without_name[content_without_name.find('\n')+1:]
        file.seek(0, 0)
        file.write(name+'\n')
        file.write(content_without_name)

Code: Select all

fullStride_out.write('multiple_strides.sto')
add_name_to_sto('multiple_strides.sto', 'multiple_strides')

[/size]