機械学習勉強ログ

機械学習に関連した勉強を記録していこうと思います。

線形代数(回転行列)とPythonを使ってアナログ時計を表現してみる。

今回のテーマ

前回は、拡大行列に続いて、回転行列についてアウトプットしてきました。 今回は、応用としてアナログ時計の針を回転させるというのをやってみたいと思います。

実は、これ以前から、やってみたかったんですよね。それも随分昔から、、。

昔話

大学生のころ(1997年頃)、当時Windows上でゲーム作りにハマっていまして、 ちょうど当時流行っていた3D技術(Direct 3D)を使って何か作りたいと思っていました。(世の中的にも3Dというだけで持て囃されるような時代でした) それで、買ったのがこの本です。

Direct3Dプログラミングガイドブック

Direct3Dプログラミングガイドブック

この本には「地球を回してみよう」 という、エンジニアだったら誰しもあこがれるような魅力的なお題が設定されていました。 自分も回せるだろうと思って読み始めたのですが、当時は3Dの原理となる数学の部分が理解できなくて 結局回せなかったのでした orz。(もちろん、サンプルは動かしましたがw)

しかし、20年たった今、「線形代数」を学んだことで「あ、たぶん自分地球回せるわw」と思いまして。 そして、いきなり問題を3次元に拡張するのも、自分にとっては無謀なので、まずは2次元で時計の針を回そうと思い、 このテーマに行き当たった訳です。

時計の針を回してみる

まずは、問題を式として、表現することから手を付けていきます。

まず、時計の針は 秒針長針短針 とありますが、これは ベクトル としてあらわせそうです。

それぞれの長さは、3, 2, 1としておきましょうか。 そうすると、0時0分0秒の位置を初期位置とすると、それぞれのベクトルは。

こんな感じで表せそうです。

で、これらの針達が 躍動 回転する。

回転の速度はそれぞれ違うので、一気に考えるのではなく、 まずは、秒針の動きから考えていくことにします。

秒針の回転を考えてみる。

秒針の座標は、秒針の初期位置からの回転運動なので、次の式になると思います。

ただ、これだと角度を与えて秒針の位置を求める関数なので、

時間(s)の関数に直したほうがよさそうです。

そう考えると、秒針60秒時計回りに1週( -2\pi) するので

 \theta(s) = -2\pi \frac{s}{60}

時間と角度の関係は上式になり、時間sの時の秒針の位置は、

 n_s(s) = R(-2\pi \frac{s}{60}) \cdot n_{s0}

あっけないほど、簡単に表すことができてしまいました

短針長針に関しては、秒針との違いだけ考えればよさそうです。

違いは次の2点なので、

  • 針の長さ(針の初期位置のベクトル)
  • 針の回転速度

針の式は以下のようになります。

  • 短針の式

    •  n_m(m) = R(-2\pi \frac{m}{60}) \cdot n_{m0}
  • 長針の式

    •  n_h(h) = R(-2\pi \frac{h}{12}) \cdot n_{h0}

Pythonでのシミュレーション

さて、各針の動きの式ができたので、これを Pythonを使って シミュレーションしていきましょう。

シミュレーションで使っている環境の説明を簡単にしていきます。

まず、シミュレーションは、Google Colaboratory上で作成しています。

Google Coraboratoryは、 Jupyter notebookのマネージド環境です。 簡単にPythonの環境を用意できる上、Jupyterなのですぐに結果が確認できるので非常に便利です。 無償で利用できます。

NumPyの使い方

Pythonで行列演算を扱うためのPythonのライブラリです。 軽く触ってみましたが、今まで勉強した行列やベクトルの操作がライブラリ関数 として提供されているので、非常にシンプルに式をプログラムに変換できて便利です。

パッケージのインポート

import numpy as np

NumPy パッケージは、npと省略して使うのが慣習のようなので、 それに倣っていきます。以下、簡単な演算の数式とPythonの簡単な対応を書きます。

横ベクトル

 A=\begin{bmatrix}
1 & 2 & 3
\end{bmatrix}

A = np.array([[1, 2, 3]])

縦ベクトル

 A=\begin{bmatrix}
1 \\\
2 \\\
3
\end{bmatrix}

A = np.array([[1,2,3]]).T # Tで転置(transpose)できる

ベクトルの内積

 \begin{bmatrix}
1 & 2 & 3
\end{bmatrix}
\cdot 
\begin{bmatrix}
1 \\\
2 \\\
3
\end{bmatrix}

np.dot(np.array([1,2,3]), np.array([1,2,3]))

ここまであれば、とりあえず時計の実装に必要な行列演算はほぼ手に入ったと思います。

Matplotlib

こちらは、データを可視化するためのライブラリです。 一般的な使い方はググってもらうとして、今回、ちょっと特殊な使い方をしているので その点だけ補足します。

FuncAnimation関数

通常 matplotlib は静的なグラフを表示するのに使いますが、今回は動的に時計の針のアニメーションさせるのに 使っています。これを実現するには、 FuncAnimation という機能を使います。

これは、下記のようにフレーム(n)を引数として、nフレーム目の静的なグラフを返却する関数を定義することで アニメーションを実現するものです。 下記に使用例を示します。

clock_animation = animation.FuncAnimation(
  fig=fig,
  func=update_needle,   # frameに基づいたレンダリングを行う関数
  interval=1000,        # 1秒ごとに描画を更新する
  frames=60*2,          # フレーム数
)

figで指定されるバッファを、funcというフレームを引数とする関数で更新していきます。

上記の例では、1秒ごとにupdate_needleを呼び出して、0フレーム~119フレームまで、合計120フレーム、 2分のアニメーションを生成しています。(clock_animation)

clock_animationは、JSプレーヤーで再生することもできるし、動画ファイルとして出力することも可能です。

注意点

  • 動的にアニメーションしているわけではなく、アニメーションデータを作ってから、データを再生する構造になっています。

    • 毎フレームごとに状態を取得するような処理は、書けません。(例えば、現在時刻を取得するなど)
    • フレーム数が多くなりすぎると、レンダリングに時間がかかることに注意。
  • レンダリング関数では、plot系の関数で全体を再描画してもよいが、公式のマニュアルにあるように、 オブジェクトのプロパティを書き換えて差分更新した方がよさそうです。

完成したアナログ時計

ソースコード

import numpy as np
from sympy import sin, cos, pi
import matplotlib.pyplot as plt
from matplotlib import rc
import matplotlib.animation as animation
from datetime import datetime
from pytz import timezone

# 針ベクトル
HOUR_NEEDLE     = np.array([[0,0], [0,1]])
MINUTE_NEEDLE  = np.array([[0,0], [0,2]])
SECOND_NEEDLE = np.array([[0,0], [0,3]])

def rotate(theta):
  return np.array([
    [cos(theta), -sin(theta)],
    [sin(theta), cos(theta)]          
  ])

# 各針の動きを表す関数
def seconds_needle(s):
  return rotate(-2 * pi * s / 60)

def minutes_needle(m):
  return rotate(-2 * pi * m / 60)

def hours_needle(h):
  return rotate(-2 * pi * h / 12)

# 針の位置をフレーム(1秒)ごとに更新する
def update_needle(frame):
  time = dt.second + dt.minute * 60 + dt.hour * 3600  + frame
  needles[0].set_data(
    np.dot(seconds_needle(time), SECOND_NEEDLE)
  )
  
  needles[1].set_data(
    np.dot(minutes_needle(time / 60),MINUTE_NEEDLE)
  )
  
  needles[2].set_data(
    np.dot(hours_needle(time / 3600),HOUR_NEEDLE)
  )
  
  return needles

def create_needle2d(ax, data, needle_width=3, color='red'):  
  # 針のLine2Dオブジェクトを生成
  needle2d = ax.plot([], [], lw=2)[0]
  needle2d.set_data(data)
  needle2d.set_linewidth(needle_width)
  needle2d.set_color(color)
  return needle2d

dt = datetime.now(timezone('Asia/Tokyo'))
fig, ax = plt.subplots()
ax.set_xlim(-3.5, 3.5)
ax.set_ylim(-3.5, 3.5)
ax.set_aspect('equal')

plt.close() # 余計なプロットがされるため追加

needles = [
    create_needle2d(ax, SECOND_NEEDLE, needle_width=1, color='red'), 
    create_needle2d(ax, MINUTE_NEEDLE, needle_width=2, color='black'), 
    create_needle2d(ax, HOUR_NEEDLE, needle_width=3, color='black')
]

# 針を1秒ごとに描画するためのイベントハンドラを設定
clock_animation = animation.FuncAnimation(
  fig=fig,
  func=update_needle,
  interval=1000,
  frames=60*2
)
rc('animation', html='jshtml')
clock_animation

colab.research.google.com

まとめ

線形代数を使って、アナログ時計を表現することができました。 とりあえず、線形代数を学んできて、やっと現実との接点が見えるサンプルができたので、 面白いですね。

地球は、いずれ回したいです😉。

現時点で把握している、時計と地球との違いは🤔、

  • 座標が3次元であること
  • 地球の回転運動について、自転と公転の2種類がある。
  • 公転運動について、楕円運動かつ面積速度が一定(ケプラーの法則)

機械学習とは離れていってしまうので、すぐにはやりませんが、いずれ気が向いたら やってみようと思っています。

工学基礎 はじめての線形代数学 (KS理工学専門書)

工学基礎 はじめての線形代数学 (KS理工学専門書)

こちらの本の方は、9章あたりまで解き進めましたので、 次回は、行列を使って連立方程式を解く ことについてアウトプットしていければと思います。