Welcome to the [fourth] NeuroMechFly video tutorial. In this tutorial, we'll cover how to run
a simulation where the simulated fly replays the behavior of a real fly that's being recorded
in a real experiment.
Again, the current version of FlyGym is 2.1.0. Future versions might deviate from this
tutorial, therefore it's recommended that you always follow the tutorials on the website instead of repeating
what I will do in this video recording verbatim.
For the purpose of this tutorial, we'll follow code tutorial notebook number two, because we are
going to run actual simulations and performance kind of matters. We're not going to do it
on Google Colab, which can be up to 50 times slower than running the notebook locally.
So I will use a local installation of FlyGym and I'm going to do it in
VS Code. So this is a cloned version of FlyGym. You can go to tutorials and
open this tutorial notebook.
Select a kernel in this case.
UV environment that we have created.
The first sub block is for it only runs if you're if you're on Google Colab,
because we are not, it simply finishes in zero second and does nothing.
And then we go to the main content of the tutorial.
One way to use or rather one advantage of running neuromechanical simulations is that we can
replicate the observable behavior from real experiments and infer latent quantities that we cannot directly measure.
For example, we can record a movement sequence in a real fly.
But that recording itself doesn't tell us, for example, ground reaction forces that the fly is experiencing.
But if we replay that behavior snippet, then we can measure ground reaction forces from the
physics simulation. But because the simulation being in silico, we have free access to basically
everything we want.
Of course, validation is crucial.
It's enough. It's not enough. But one validation that we might use is to see if
the observed behavior is faithfully replicated in the simulation.
So in this particular example, we will use a, we'll use FlyGym to replay a
recording that's been recorded in our, in our fly behavior quantification system called Spotlight.
You can read more about it here.
In essence, this experimental system allows us to record behavior of a freely
walking fly at very high resolution, but also recording muscle activities at the same time.
In this case, we don't use muscle data. We only use the behavior data.
So the first thing we want to do is to import from flygym_demo.spotlight_data,
this pre-programmed utility class called Motion Snippet.
And this is going to load a default data snippet with metadata that you can see
from this object.
So in this case, we have behavior data on all six legs. We have these degrees
of freedom per leg.
The data is recorded at 330 frames per second. This is the name, et cetera, etc.
And in this data set, we have joint angles. So this many frames over two seconds,
six legs and seven degrees of freedom per leg.
We have egocentric X, Y, Z position of each leg joint or rather. Is it?
Yes. So we tracked five key points per leg and we have six legs. So we
have 30 key points and for each of them, we have X, Y and Z.
The difference between forward kinematics and raw prediction is not as critical. We're just going to
skip it in this video tutorial.
We can visualize the data that we have.
For example, here we're visualizing the X, Y, Z positions of the tip of the leg
or the claw of the left front leg.
So over time, this is how the fly is, how this claw is moving in
space in egocentric coordinates relative to the thorax of the fly.
We can also plot this in 3D, but we only plot a single frame. But in
this frame, this is the pose of the fly.
So if we don't look at the colored lines for now and only look at the
dashed line, this is the results from our pose estimation model from the experiment directly.
The problem with it is that you have X, Y and Z position for each joint
key point, but you don't have joint angles.
And in fact, it's difficult to have it because there could be errors in the
3D pose prediction. The distance between two key points,
which reflects the physical length of a leg segment, which should remain constant over time, might
change over time because of these inaccuracies.
So we apply a step called inverse kinematics, which basically involves fitting or optimizing for joint
angles from raw X, Y, Z predictions.
And if we reapply the set of joint angles, we can reconstruct the pose of each
leg showing in colored lines and see that even though it's not perfect, it agrees mostly
with the raw predictions.
And similar to the previous plot, we can generate a figure showing the joint angles over
time. In this case, it's the seven degrees of freedom of the left front leg.
Now we basically repeat what we did in tutorial one. I'm
just going to skip the details of it and let the viewer refer to tutorial one
if it's not unclear.
And different from tutorial one, we're also going to add leg adhesion for this simulation.
So in the biological reality of fly locomotion is that flies have adhesive structures at the
tips of their legs, which allow them to maintain a stronger grip on surfaces and work
on walls and ceilings.
Again, we define a spawn position and orientation. We create a flat ground world and attach
the fly to the flat ground world.
Once that is done, we can create a simulation object imported from FlyGym. So from now
on, this is new compared to tutorial one.
We import a simulation, we create a simulation on the world, and we can also set a renderer.
Recall that when we created the fly, we attached a camera to it. Where is it?
Here. We attached a tracking camera to the fly.
And now once we have created the simulation, we can set a renderer. Let's just run that.
Now once we have the simulation created, the next thing we need to do is to
extrapolate the experimental data to the time grid of the simulation.
So for physics simulation, we typically use very, very small time steps. For example, 0.1 millisecond.
This is because if you integrate the physics at a coarser level, then the physics becomes
very unstable.
However, the behavior data is not recorded at 1000 FPS or rather 10,000 FPS.
And therefore what we need to do is to extrapolate, excuse me, interpolate the experimental data
to this timestamp.
Again, this function has been, the interpolation function has been provided as part of
this utility class called MotionSnippet from flygym_demo.spotlight_data.
Once we created the snippet class or rather the snippet object, we can simply call
get_joint_angles by giving it a desired output time step of 0.1 millisecond.
And by specifying the order of actuators from NeuroMechFly from the fly simulation, we can
get an array with 20,000 time steps to simulate. Recall that our behavior snippet is
2 seconds long and we are simulating at 10,000 hertz. So this is what we
would expect.
And the shape of target joint angles is 20,000 by 42 because we have 42
actuators, 7 per leg and 6 legs.
Now we have everything we need to run the simulation.
What we can do is that we can reset the physics simulation.
We can create some empty buffers to record the state of the simulation as we go.
And to start, we can already enable the leg adhesion actuators.
We set them to a constant of 1. So in other words, the leg adhesion forces
are always on as we simulate.
So when the fly lifts and like it, it actually needs to apply a force that
is stronger than the adhesion force in order to lift the legs.
We can warm up the simulation. Basically, what this does is that it drops the fly
from the spawn point and the fly lands on the ground.
And I think for a 0.05 second time period, the latent forces of spawning the
fly, the fly landing on the ground, mostly settle.
And after this period, we can start running our simulation from a cleaner state.
Now, what we do next is to implement our main physics simulation loop.
So we want to write a for loop that iterates over this many time steps, in
our case 20,000.
If you don't know this already, using the tqdm package, you can replace range with trange,
which allows you to basically get a progress bar as you iterate over this loop,
which is handy. And within each iteration of the loop, we fetch the target joint angles.
Recall that we have this this array called joint angles, interpolated to NeuroMechFly time that
is 20,000 by 42.
And we are iterating over the simulation time steps. So for each simulation time step, the
target joint angles is simply this, with the step index being applied to the
0th dimension of this array.
Once we have these target angles, we can call the set_actuator_inputs method to the
simulation object.
We apply it to the fly identified by its name. We apply it to the position actuators.
So actuator type here, I think it's defined a little bit earlier.
It's simply position, yes, here.
And we set the target angles, target positions.
Once we have set the target joint actuator target states, we can step the physics simulation.
You can run this sim.step with no argument. Or alternatively, you can run step_with_profile,
which simply runs some time counter to report the model performance to you.
But in reality, you can completely, it's recommended that you run sim.step
because it's a little bit faster, because in order to time the simulation,
there's a little bit of overhead.
Now we have told the simulator what we want the actuators to do.
We have stepped the physics forward by a time step. The next step,
what we can do now is to get the information that we want to have, such
as joint angles, actuator forces, and positions of different key points or sites on the fly.
And we can write these values to these three pre-existing buffers that we created
before this loop.
And finally, we have provided a render_as_needed function or method on the simulation class
that you can call.
It captures the frame if necessary.
The reason why a frame is not always captured for every time step is simply that
we're running 20,000 time steps and we don't need an output video of 20,000 FPS.
It's only going to make the simulation very, very slow.
So here, this method encapsulates the calculation of at which simulation time steps
a render is needed.
Again, we can use render_as_needed with profile to get some timing information.
And now let's run the simulation. At
the progress bar, we're simulating 4,000 iterations per second.
In other words, we're simulating roughly half a second of physics per second.
And because we used the with profile versions of step and render_as_needed, we can
call this little method to get a table indicating how fast our simulation is.
You can see a breakdown that 39% of the simulated time is used to actually
advance the physics simulation.
It's running at just a bit more than 10,000 iterations per second.
In other words, it's about 107% of real time.
But the simulation is spending 61% of the time rendering.
So 61% of the time is just spent on video generation. And overall, with rendering,
the overall throughput is 42% of real time.
There's a note saying that out of the 20,000 steps, only 250 frames were rendered.
You can actually adjust this by going back to the set_renderer call.
Here you can actually see that by default we're generating the video that's played at 20% speed.
So five times slower than real time and we're rendering the video at 25 frames per second.
In other words, we are rendering 125 frames —
25 times 5 — per second simulated, which would
make sense because over two seconds we simulated 250 frames.
By adjusting the playback speed and output FPS, you can potentially reduce the number of frames
that need to be rendered.
Now going back to the renderer class that we created when we called sim.set_renderer.
We can use the show_in_notebook method to display the video describing the simulation.
So two seconds of simulation played at 20% speed, we have a 10-second video.
We see a somewhat natural-looking behavior of flywalking and executing a turn.
Alternatively, you can also call this save_video method on the renderer.
So after we save the video, instead of in addition to showing it in the notebook.
Now it's an mp4 file on your hard drive.
Now, recall that we didn't just set the actuator input and simulate and render.
We also recorded joint angles, forces applied by position actuators and key point positions as we simulated.
So now we can create a plot that compares the trajectory of joint angles to the target.
So in this case, the blue dotted lines indicate the target positions.
If we go back to this interactive viewer, it's basically where the sliders are, where we
want the joints to be.
And the orange line indicates the actual position in the simulation.
In other words, where the green triangles would be in our interactive demo.
We see that it's pretty close, which is good.
And we can also plot the pose of the fly over time by looking at the
key point positions in X, Y, Z.
So in this case, we selected just a couple of snapshots.
And for each of them, we plotted all six legs.
And there you go.
We have taken some data from a real experiment and replayed it in FlyGym.
For more information, you can actually see the Spotlight —
where is it —
yes, this website describing our Spotlight system.
Let's see a brief preprint on bioRxiv.
There's also a little bit of discussion on replay of behavior in NeuroMechFly.
And that's it in the next tutorial.
We'll basically do the same thing, but this time with GPU acceleration. Thanks for watching.