三次元空間における2つのベクトルのなす角度の計算

とある Python プログラムのお披露目のサイトで，「三次元空間における2つのベクトルのなす角度の計算」があった。

パッケージの関数を使いまくりで，わかりやすいのかわかりにくいのか。

import numpy as np

def angle(a, b, c):
    # ベクトルを定義
    vec_a = a - b
    vec_c = c - b

    # コサインの計算
    length_vec_a = np.linalg.norm(vec_a)
    length_vec_c = np.linalg.norm(vec_c)
    inner_product = np.inner(vec_a, vec_c)
    cos = inner_product / (length_vec_a * length_vec_c)

    # 角度（ラジアン）の計算
    rad = np.arccos(cos)

    # 弧度法から度数法（rad ➔ 度）への変換
    degree = np.rad2deg(rad)
    return rad, degree

http://w3e.kanazawa-it.ac.jp/math/category/vector/henkan-tex.cgi?target=/math/category/vector/naiseki-wo-fukumu-kihonsiki.html

にあるように，式自体は簡単なのと，numpy.ndarray は要素ごとの演算ができるのが特徴なので，以下のように簡単になる。

内積：　sum(x * y)
ノルム：np.sqrt(sum(x**2))
ノルム：np.sqrt(sum(y**2))

長すぎる変数名というのも，読みにくいなあ。

def angle2(a, b, c):
    x = a - b
    y = c - b
    cos_theta = sum(x * y) / np.sqrt(sum(x**2)) / np.sqrt(sum(y**2))
    rad = np.arccos(cos_theta)
    degree = rad * 180 / np.pi
    return rad, degree

>>> a = np.array([0, 1, 2])
>>> b = np.array([10, 20, 30])
>>> c = np.array([5, 7, 9])

>>> print(angle(a, b, c))   # (0.09657114452339467, 5.533119003938428)
>>> print(angle2(a, b, c))  # (0.09657114452339467, 5.533119003938429)

>>> a = np.array([2.123, 3.215, 4.325])
>>> b = np.array([1.547, 2.487, 6.543])
>>> c = np.array([3.234, 1.124, 4.328])

>>> print(angle(a, b, c))   # (0.8548133509459923, 48.97719728064064)
>>> print(angle2(a, b, c))  # (0.8548133509459923, 48.97719728064063)

range() と numpy.range() -- Python

range と numpy.arange は速度が大違い。numpy.range を使うべし。

import numpy as np

def Eratosthenes(n):
    """ n 以下の素数をエラトステネスの篩で生成 """
    tbl = np.arange(n+1, dtype=int)
    tbl[1] = 0
    for i in range(2, int( np.sqrt(n))+1):
        if tbl[i]:
            tbl[range(2*i, n+1, i)] = 0
    return tbl[tbl > 0]

def Eratosthenes2(n):
    """ n 以下の素数をエラトステネスの篩で生成 """
    tbl = np.arange(n+1, dtype=int)
    tbl[1] = 0
    for i in range(2, int(np.sqrt(n))+1):
        if tbl[i]:
            tbl[np.arange(2*i, n+1, i)] = 0
    return tbl[tbl > 0]

from time import time as t
s = t();a = Eratosthenes(10000000);print(t()-s)　# 4.958sec
s = t();b = Eratosthenes2(10000000);print(t()-s) # 0.268sec

な，な，なんと！18 倍も違う。

演算の意味をよく考えて，簡潔な記述を

Python のプログラム例を披瀝するページで以下のようなプログラムを見かけた。

余計なカッコがあったりして，演算がごちゃごちゃしているのと，「整数除算演算子」と「整数剰余演算子」をご存じないようなのが仇をなしているのだ。

# 秒を入力するとHH:MM:SSの形式の文字列で出力する関数
def sec_to_time(sec):
    h = int(sec / 3600)
    m = int((sec - (3600 * h))/60)
    s = int(sec - (3600 * h) - (60 * m))
    rtx = '{:02}h {:02}m {:02}s'.format(h, m ,s)
    return rtx

これを以下のように書き換えるとわかりやすい。

def sec_to_time2(sec):
    h = sec // 3600
    m = sec % 3600 // 60
    s = sec % 3600 % 60
    rtx = '{:02}h {:02}m {:02}s'.format(h, m ,s)
    return rtx