MyoSuite Tutorial No.3

# はじめに
本記事は以下の人を対象としています．

- 筋骨格シミュレータに興味がある人
- MuJoCoに興味がある人

:::note warn
本記事では[3_Analyse_movements.ipynb](https://github.com/MyoHub/myosuite/blob/main/docs/source/tutorials/3_Analyse_movements.ipynb)をGoogle Colabで実行した際の手順と結果に関してまとめる．
:::

# MyoSuiteとは

https://sites.google.com/view/myosuite

MuJoCoを学習モジュールとして組み込んだ筋骨格シミュレータ

:::note info
2024年にICRAにてworkshopを展開
https://sites.google.com/view/myosuite/myosymposium/icra24
:::

# Tutorial Notebook

https://github.com/MyoHub/myosuite/tree/main/docs/source/tutorials

# 実施事項
Tutorialに記載のコードを1つずつ試して, Errorになったところは都度修正版を記載する.

## ライブラリのinstallと環境設定

```bash
!pip install -U myosuite
!pip install tabulate matplotlib torch gym==0.13 git+https://github.com/aravindr93/mjrl.git@pvr_beta_1vk
!pip install scikit-learn

%env MUJOCO_GL=egl
```

## ライブラリのimportと結果表示の関数

```python
from myosuite.utils import gym # importが上手くいくと"MyoSuite:> Registering Myo Envs"と表示される
import skvideo.io
import numpy as np
import os

from IPython.display import HTML
from base64 import b64encode

def show_video(video_path, video_width = 400):

  video_file = open(video_path, "r+b").read()

  video_url = f"data:video/mp4;base64,{b64encode(video_file).decode()}"
  return HTML(f"""<video autoplay width={video_width} controls><source src="{video_url}"></video>""")
```

## modelのloadとpolicyの準備
```bash
# stable-baselines3のinstall
!pip install stable-baselines3
```

```python
from stable_baselines3 import PPO
env = gym.make('myoElbowPose1D6MRandom-v0')

env.reset()
pi = PPO("MlpPolicy", env, verbose=0)

pi.learn(total_timesteps=1000)
```

stable_baselines3についてはnpakaさんの資料を参照

https://note.com/npaka/n/nb615fd590274

## actionによる筋負荷データの出力

```python
frames = []
data_store = []
for _ in range(10): # 10 episodes
    for _ in range(100): # 100 samples for each episode
        frame = env.sim.renderer.render_offscreen(
            width=400,
            height=400,
            camera_id=0
        )
        frames.append(frame)
        o = env.get_obs()
        a = pi.predict(o)[0]
        next_o, r, done, *_, ifo = env.step(a) # take a random action

        data_store.append({"action":a.copy(),
                            "jpos":env.sim.data.qpos.copy(),
                            "mlen":env.sim.data.actuator_length.copy(),
                            "act":env.sim.data.act.copy()})

env.close()

os.makedirs('videos', exist_ok=True)
# make a local copy
skvideo.io.vwrite('videos/temp.mp4', np.asarray(frames),outputdict={"-pix_fmt": "yuv420p"})
show_video('videos/temp.mp4')
```

:::note alert
残タスク
- action, jpos, mlen, actがそれぞれ何を示すのか記載
:::

出力結果

![3_analise_movements-default.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2043769/a4125850-b1ed-ebe6-4135-770b10eaa885.gif)


## 解析結果の表示

```python
import matplotlib.pyplot as plt
from sklearn.decomposition import NMF

def VAF(W, H, A):
    """
    Args:
        W: ndarray, m x rank matrix, m-muscles x activation coefficients obtained from (# rank) nmf
        H: ndarray, rank x L matrix, basis vectors obtained from nmf where L is the length of the signal
        A: ndarray, m x L matrix, original time-invariant sEMG signal
    Returns:
        global_VAF: float, VAF calculated for the entire A based on the W&H
        local_VAF: 1D array, VAF calculated for each muscle (column) in A based on W&H
    """
    SSE_matrix = (A - np.dot(W, H))**2
    SST_matrix = (A)**2

    global_SSE = np.sum(SSE_matrix)
    global_SST = np.sum(SST_matrix)
    global_VAF = 100 * (1 - global_SSE / global_SST)

    local_SSE = np.sum(SSE_matrix, axis = 0)
    local_SST = np.sum(SST_matrix, axis = 0)
    local_VAF = 100 * (1 - np.divide(local_SSE, local_SST))

    return global_VAF, local_VAF

act = np.array([dd['act'] for dd in data_store])

VAFstore=[]
SSE, SST = [], []

sample_points = [1,2,3,4,5,10,20,30]
for isyn in sample_points:
    nmf_model = NMF(n_components=isyn, init='random', random_state=0);
    W = nmf_model.fit_transform(act)
    H = nmf_model.components_

    global_VAF, local_VAF = VAF(W, H, act)

    VAFstore.append(global_VAF)

plt.plot(sample_points,VAFstore,'-o')
plt.xlabel('Number of Muscle Synergies')
plt.ylabel('Explained Variance R^2')
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
```

:::note info
- ```sample_points```はNon-negative Matrix Factorization(NMF)実行時の次元数
    - 一応>1000でも実行はできる
- Muscle Synergyとは筋グループの数
    - ```sample_points```の値が増えるということは筋グループの数が増えるということ
    - 今回だとMuscle Synergyは10未満が適度に分解されていい値？
:::

出力結果
![global-VAF.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2043769/a3b9ad91-daf2-c07f-2230-cbbd08f613ea.png)

![local-VAF.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2043769/dcb3c549-5597-67b5-7051-ec31812957b1.png)


# 参考資料

https://qiita.com/takechanman1228/items/6d1f65f94f7aaa016377

https://movement-design.com/%E7%9F%A5%E3%81%A3%E3%81%A6%E3%81%8A%E3%81%8D%E3%81%9F%E3%81%84muscle-synergy