科学と家事とプログラミング (python を中心に)

python 温度計測 湿度計測 DS18B20 USB9097

角度を扱うプログラムは、要注意なのです

天文計算や宇宙機の姿勢軌道制御では、角度の計算が欠かせませんが、 角度を扱う時は、色々と気が抜けません。 代表的には以下の3点。

  • 単位
  • 極性
  • 定義域 (周回の不確定性の扱いを含む)

単位

初等教育や日常生活では、deg を使う場面が多いけれど、プログラミング言語(C 言語や python)の 数値演算ライブラリの引数は rad に統一されています。 角度の単位の間違いを減らす一つの方法は、「deg 単位の変数は使わない」 という方針を禁欲的に守ることかな...と思っています。 例えば、引数で与えられた角度 [deg] と sin の値を表示するなら

import sys
import math
DEG = math.pi/180

x = float(sys.argv[1])*DEG # 変数は rad で保持
s = math.sin(x)
print("x %8.3f [deg]  sin(x) %6.3f [-]" % (x/DEG, s)) # 出力時に deg 単位に変換

x を deg 単位で保持した方が無駄が少なそう... といった、ケチな事を考えてはいけません。 ここでブレると、いずれ間違いが生じます。

もちろん、例外はあって、10deg 単位のメッシュを切りたいとか、 そういった場合には deg 単位の変数をその目的に限定して使うことはありますが、 角度一般の変数としては、必ず rad を使うことにしています。 プロジェクトや開発チームで統一されているとなお吉。

参考までに、print 出力の x/DEG という表記は、 x を deg 単位で表した値という意味を直接表現していて分かり良い表記と思っています。 入力側の x = 3*DEG も、1 deg が 3 つ分の角度ということで、こちらも 時々みかける deg2rad や rad2deg といった定数を乗ずる変換より、コードの 意味が明瞭と思います。このあたりは慣れや好みの問題かもしれませんけど。

定義域

定義域や周回の不確定性の扱いが問題なる例として、位相の区間の表現を取り上げます。 与えられた位相 x が位相区間 a, b に含まれるかどうかを判定することを考える。 単位円上の2点 a, b の位相角で位相の区間が与えられるものとする。 (問題を簡単にするため、ここでは a, b で定まる狭い側の円弧を採用する) 数直線であれば、a < b となるように選んで、a < x < b を判定すれば良いけれど、 相手が位相の場合には、そう簡単ではありません。 問題は -720, -360, 0, 360, -360, 720 [deg] といった位相角が、位相としては同じだということにつきます。 a, b の符号や大小関係を駆使して、場合分けによる判定式を構築できなくはないけど、 なかなか複雑になるので、出来れば避けたいところです。 ここでは、位相の定義域や a, b の数値的な大小関係に依存しない解法を考えてみます。 理由は、設計も検証も楽だから。

成す角で勝負

簡単に思いつくのは、a と b の中点 c をとって、c と x の成す角を計算する方法でしょうか。 c と a の成す角を基準値として判定できそうです。位相ではなくベクトルの成す角を 使って定式化しているので、位相の定義域や周回の不確定性には影響されません。 中点 c を計算するときは、かならずしも円周上で考える必用はなくて、ベクトル a, b の 単純平均(を規格化)で十分な場合が多いですね。

判定用に局所座標を導入する

他には、a または b を位相の原点に採った座標を導入して、その座標で判定することでしょうか。 よく設計されたベクトル行列演算ライブラリがあれば、こんな感じ

u = vector(az=a, el=0) # a 方向単位ベクトル
v = vector(az=b, el=0) # b 方向単位ベクトル

x = u.unit()       # a 方向を位相原点にとる
z = (u * v).unit() # * は外積
y = (z * x).unit() # y 方向 (位相=+π/2)
loc_fm_ref = mat(x, y, z)  # 基準座標 ==> 判定局所座標

v_loc = loc_fm_ref*v_loc
x_loc = loc_fm_ref*vector(az=x, el=0) # x 方向

# v_loc.az() と x_loc.az() の大小関係を調べる。具体的な判定規則は、
# .az() の仕様に依存する

計算量は多いけれど、a や b の定義域や極性をプログラム中で定義しなおす感じですね。 結果の良し悪しや限界も分かりやすいので、利用しやすい方法と思います。

定義域については、0~360deg と -180~180deg のどちらかという場合が多いけれど、 周回の不確定性を(部分的に)許容する場合もあるので、必ず確認が必用です。

極性

位相の測り方は、CCW が多いけれど、CW も無くはないので、明示されていないときは 必ず確認を。

桁落ち誤差は、時に牙をむくのです

数値計算には誤差がつきものですが、扱いをあやまると酷い目に会うことになります。 特に桁落ち誤差には注意が必要です。 天文計算の分野で有名な例は、universal kepler 方程式を解く時に出てくる stumpff 関数の計算式でしょうか。例えば、こんな感じ。

 \frac{x-sin(x)}{x^{3}}

見るからに危なそうな形をしていますね。そのまま計算した時の誤差の様子を調べてみましょう。

f:id:sken20k:20180311191455j:plain

横軸は対数にとってあります。 ゼロに近いところの正解は 1/6 のはずですが、x が 1e-7 より小さいところで値が暴れて、 まともに計算できなくなるのが分かります。 「誤差がある」という言葉を聞いたときに想像する誤差の大きさは、分野によっても違うでしょうけど、 大きくても 3割ぐらいを想像する人が多いのではないでしょうか。 桁落ちに代表される計算機の数値誤差は、ときどきデタラメな返として現れるので、注意が必要です。

桁落ちの対処方法はケースバイケースです。上記の場合、x が小さい領域では級数展開に切り替える ような工夫がなされます。例えば、奥村晴彦さんの「C言語によるアルゴリズム」あたりに、 双曲線関数を対象としたチューニング例が出ています。

実際の場面では、元の問題(universal kepler eq.) を解くときに、stumpff 関数の使用そのものを避けるよう な研究(定式化)もあるようです。

ケプラー方程式をニュートン法で解くときの開始点の選び方は、悩みが尽きないのです

ケプラー方程式

 M = U - e*sin(U)

を解くとき、よくニュートン法が使われます。 離心率が大きい長楕円軌道の場合、開始点の選択を誤ると繰り返し計算がバタついて、 収束しないことがあります。 ケプラー方程式は、下に凸な関数なので、 大きい側から近づくと動作が安定します(作図すれば明らかですね?)。

とにかく安定させたいなら

範囲を 0~π に限定したうえで、開始点としてπを選べば、安定した動作が期待できます。 繰り返し回数が多くなりそうですが、安定した動作には代えられません。 ニュートン法は二次収束するので、開始点が離れていても、すみやかに正解に到達します。

もう少し正解に近い初期値が欲しいなら

平均近点角と離心近点角の差に外接する3つの直線を使って、 正解より少し大きい初期値を求めることができます。

 U0 = min(\frac{M}{1-e}, M+e, \frac{M+πe}{1+e})

これは、Nijenhuis の初期値と呼ばれます。 式の意味が幾何学的に明確なところがお勧めポイントです。

近点角が小さいところの相対誤差が気になるなら

Nijenhuis の初期値が万能かというとそうでもなくて、e が 1 に近くて(例えば e=0.9 以上) M がゼロに近いところだと、初期値の相対誤差が大きくなります。 この誤差を抑えるのに、以下の式が役に立ちます。

 M + (6M)^{1/3}

離心率 e = 1 を仮定した近似解をベースにしています。実際の離心率を反映 すれば、さらに改善が見込めます。ただし、近似解は正解より大きい値を返す 必要があります(イタレーションの安定な収束を保証するため)。初期値の計算に三乗根を 使うのはやりすぎの感もありますが、全区間に渡って相対誤差を一定以下に抑えることができる点は魅力的です。

f:id:sken20k:20180313215559j:plain  U3:nijenhuis U4:三乗根

M ゼロに近い領域の挙動を確実におさえておきたいならば、 M を対数でサンプリングして、上のグラフを相対誤差でプロットすると良いと思います。

収束判定の一工夫

上で紹介した初期値を使えば、大きい側から正解に近づきます。 そうはいっても数値誤差の関係で修正方向が反転することもあります。 そういう場合に備えて、「解の振動を検出したらそこで停止する」のも一案です。

python の デコレーターを使ってみよう !

定義域(値域)を 0~π に限定した関数を実装したあと、 定義域(値域)を-π~πに拡張することを考えます。 python なら、デコレーターを使って、ちょっとおしゃれに実装できます。

def deco_expand_sign(f):
    # 符号拡張
    def _f(x):
        if x < 0:
            y = -f(-x)
        else:
            y = f(x)
    return _f

@deco_expand_sign
def keq_eli(mu, m):
    # m
    ...
    return u
    

deco_expand_sign は、双曲線にも放物線にも適用できるあたりががミソです。 楕円軌道の場合は、+/-n 周回の不確定性を解決するデコレーターの採用も見通しが良いですね。 私的には、戻り値に +/-n 周回の情報を残した設計(定義域も値域も -∞~∞)が好みです。

ケプラー方程式の解法に関しては、天文台の福島登志夫さんによる1980年代の仕事が完成度高いです。 そのぶん、実装も検証も片手間で...とはいきませんので、採用するにはそれなりの覚悟がいります。

時刻を秒の少数まで表示するのは、ちょっと難しいんです

時刻表示の落とし穴

天文や宇宙関連のプログラムを作っていると、時刻表示で落とし穴にはまることがありますね。 例えば、秒の少数1桁まで表示したいのだけど、四捨五入されて秒が 60.0 になってしまったり...

>>> "%02d:%02d:%04.1f" % (23, 59, 59.95)
'23:59:60.0'

事前に四捨五入をチェックすれば大丈夫か ?

ありがちな対応策は、事前に四捨五入をチェックすることかと思います。 こんなコードでしょうか

if s >= 59.95:
    s = int(s) + 1 # 少数以下を切り捨てて 1 秒追加
    # 必用に応じて、繰り上がりを上位に伝搬してから文字列表現を編集...

これで、回避できそうですが、計算機誤差(計算機ε)を考慮すると 59.95 のまるめ誤差も気になります。 たとえば(?) 、少数 3 桁の時は、繰り上がりがなかったりして、ちょっと不安になります。

>>> "%02d:%02d:%06.3f" %(23,59,59.9995)
'23:59:59.999'

まるめ誤差はやっかいですね。

やってみて、だめならやり直す!

こういう時は、実際に文字列を作ってみて、繰り上がったかどうかを調べるのも一つの方法です。

if float('%4.1f' % s) >= 60:
    s = int(s) + 1 # 少数以下を切り捨てて 1 秒追加
    # 必用に応じて、繰り上がりを上位に伝搬してから文字列表現を編集...

試しに文字列を作るのは無駄が多い感じもしますが ^^; やりたいこと(文字列変換の四捨五入が発生するかどうかで場合分けしたい)を直接表現したコードに見えてきませんか?処理系が有する書式変換の細かな仕様を気にする必要もない。 少なくとも、検証の観点から有利な設計と思います。

落穂ひろい

繰り上げが発生した場合は、必用に応じて 分,時,日,月,年 まで上位に繰り上げ を伝搬しないといけません(うるう年の判定も必要)と、時刻の印字は、なかなか、奥が深いですね。 秒の少数以下を表現できるライブラリがあれば、それを使った方が良いけれど、 うるう秒の補正や表示(60.0~60.999...)まで対応したものはなさそう。 自前で実装する場合は、秒の整数部と少数部を完全に分離して管理することをお勧めします。 以上、時刻の表示は、秒の少数の表示上の繰り上がりが年にまでおよぶ可能性がある。 ということで、ちょっと難しいんです。

超構造化(09) プログラムは手順の集まりか?

手続き vs 宣言

「プログラムは手順の集まりである」という表現は、賞味期限つきの標語の良い例だ。 初めてプログラムを書くときに、この助言が役立つ場面は確かにあると思う。 しかしその一方で、さまざまなプログラミング言語が目指すものは 「人間から見た分かりやすさ/扱いやすさ」であり、そこに通底する主題は、 「いかに計算機のように考えないで済ませるか」ということに他ならない。 通常の構造化プログラミングの範囲は超えるが、 python を用いた関数型プログラミングっぽい手法を紹介する。

無限ストリーム

例として、減衰振動のシミュレーションを取り上げる。 力学系を表す微分方程式を関数としてを動的に生成して、 状態量を無限ストリームとして与える。 こういう技法を適用すると、いたって宣言的なプログラムが出来上がる。 python の場合、無限ストリームはジェネレータと呼ばれる。 文法的には、return の代わりに yield という構文を使って値を返す。 例題中の関数 states() は yield t や x を返すが、これらの t, x は、 参照される時に(必要になって初めて)生成される。 一方、使う側では for t, x in stream: という形で、stream が既に出来上がっている (完成している存在である)かのように記述できる。 つまり、未来を先取りできるというわけだ。

コードは、scheme の文法書(r4rs)の末尾にある例題を真似たものである。

"""減衰振動(LCR 回路)の時間シミュレーション"""

class SVEC:
    """ 状態ベクトル (rk4 が要求する演算のみ実装) """
    def __init__(self, x): self.x = x[:]
    def __getitem__(self, j): return self.x[j]
    def __add__(self, o): return SVEC([a + b for a, b in zip(self.x, o.x)])
    def __rmul__(self, f): return SVEC([f*a for a in self.x])
    def __div__(self, f): return SVEC([a/f for a in self.x])

def rk4(t, x, f, h):
    """ 積分ステップ h 後の状態を返す。f は力学系(微分方程式)。"""
    f1 = h*f(t, x )
    f2 = h*f(t + h/2, x + f1/2)
    f3 = h*f(t + h/2, x + f2/2)
    f4 = h*f(t + h , x + f3 )
    return x + (f1 + 2*f2 + 2*f3 + f4)/6

def damped_oscillator(L, C, R):
    """ 減衰振動する系(微分方程式を表す関数)を生成して返す。"""
    def f(t, s):
        v, i = s # s = [C_電圧, L_電流]
        return SVEC([-(i + v/R)/C, # d(s[])/dt = dv/dt
        v/L ]) # di/dt
    return f # 上で生成した関数を返す

def propagation(intg_method, f, h):
    """ intg_method を使った刻み h の伝播 => next_state() """
    def next_state(t, x):
        x = intg_method(t, x, f, h)
        t = t + h
        return t, x
    return next_state
    
def states(derivative, t, x, h):
    """ states は、状態量の無限ストリーム(微分方程式 derivative, 初期値t, x, 刻み h) """
    next_state = propagation(rk4, derivative, h)
    while 1:
        t, x = next_state(t, x) # 次の状態は、直前の状態で決まる。
        yield t, x

if __name__ == '__main__':
    """ 状態量の無限ストリーム """
    s = states(damped_oscillator(L=1e3, C=1e-3, R=10e3),
                0.0, # 初期時刻
                SVEC([1.0, 0.0]), # 初期状態
                0.1) # 刻み幅

    # ストリームを印刷
    for t, x in s:
        if t >= 100: break
        print t, x[0], x[1]

f:id:sken20k:20180212175350j:plain

python 3.6 (32bit on Windows) に numpy / matplotlib をインストール

numpy の *.whl を入手

scipy 使うなら mkl 版の numpy が吉(ダウンロードにちょっと時間がかかった)。 *.whl を、https://www.lfd.uci.edu/~gohlke/pythonlibs/#numpy から入手した。 (pypi の numpy は mkl ではない)。

インストール

念の為、管理者モードでコマンドプロンプトを立ち上げて pip でインストール。

python -m pip install --upgrade pip
python -m pip install wheel
python -m pip install c:\FTP\numpy‑1.14.0+mkl‑cp36‑cp36m‑win32.whl
python -m pip install matplotlib

64bit 環境に移行してもも良いのかも

python 3.x on windows で numpy や scipy を使う場合、 少し前まで 64bit のバイナリパッケージが入手しづらかった。 それで 32bit 環境を選んでいましたが、今は 64bit 環境のバイナリパッケージが安定して入手できるみたいですね。 そろそろ 64bit に移行した方がよいのかもしれない。

python 3.6 on windows 仮想環境~pypi 登録まで

仮想環境は venv

構築 python -m venv myenv
開始 Scripts\activate
終了 deactivate

標準 lib の venv で構築できるのでそれを使う。

(virtualenv は、python 2.x と 3.x を同居したいときに使う)

配布物の生成

構築 python setup.py sdist
python setup.py install で試験もできるが pip install dist/package.tar.gz が吉か ?

pypi は移行期にあるもよう(2018/01)

pypi は新旧の移行期にあるようで情報が混乱している (2108/01 現在)。新旧の URL は以下の通り。

試験 test.pypi.org testpypi.python.org
本番 pypi.org pypi.python.orgi

pypi の登録/装備(インストール)

pypi (pypi.org) への登録は、twine の利用が推奨されている。 リポジトリを %HOMEPATH%.pypirc に記述するとコマンドラインが簡略化できるが、 新旧の pypi がゴッチャになりやすいので、 pypirc による指定を無効にして(行を消去して)、 コマンドラインで明示的に指定した方が混乱が少ないように思う。

test登録 twine upload --repository-url https://test.pypi.org/legacy/ dist/*
test装備 pip install --index-url https://test.pypi.org/simple/  --extra-index-url https://pypi.org/simple  パッケージ

本番登録 twine upload dist/*
本番装備 pip install パッケージ

おまけ

pypi は初回の upload 前の register が不要。

https://packaging.python.org/guides/migrating-to-pypi-org/
"setup.py register" command prior to the first upload is n longer required.