pythonの区間演算ライブラリ進捗2016-01-05
実は大体できた(唐突)。
今回はvmathクラスの紹介
import pint as pn
from pint import vmath
a = 2.
itv_a = pn.interval(a)
print(itv_a)
[2.0,2.0]
vmath.sqrt(a)
1.4142135623730951
floatを代入するとfloatが返ってくる
vmath.sqrt(itv_a)
[1.414213562373095,1.4142135623730951]
intervalを代入するとintervalが返ってくる
じゃあ、少し変なことをしてみる
itvmath = pn.vmath()
itvmath.arg_type = pn.interval
このitvmathのmath関数を使う時は代入する前にinterval型に変換するように設定
itvmath.sqrt(a)
[1.414213562373095,1.4142135623730951]
返ってくるのはinterval型になる。
itvmath.sqrt(pn.interval(a))と同じ
itvmath.sqrt(itv_a)
[1.414213562373095,1.4142135623730951]
itvmath2 = pn.vmath()
itvmath2.return_type = pn.interval
itvmath2.sqrt(a)
[1.4142135623730951,1.4142135623730951]
return_typeはreturnする直前に変換を行う
return return_type(answer)となる
そのため、変換するタイミングが違うitvmath.sqrt(2.)とitvmath2.sqrt(2.)は結果が少し違う
このようにmath関数の振る舞いを自由に変更することができるのがvmathクラスの特徴です。
今はinterval型とビルトインされている型にしか対応してませんが、そのうち増やす予定
メモ:as_interger_ratio()では正確に元の分数には戻せない
as_interger_ratio()では正確に元の分数には戻せない
a = 0.1 を区間型に代入する際にうまく包含できないので、色々実験してみた。
何が問題なのか?
import pint as pn
x=pn.interval(0.1)
と、2進数小数で表現不可能な数字なので下限と上限は
format(x.inf,'.17g')
'0.10000000000000001'
format(x.sup,'.17g')
'0.10000000000000001'
となり、0.1はうまく包含できない
解決策になりそうだったもの
pythonのfloat型にはas_integer_ratio()という関数があり、
a = 0.1
x,y = a.as_integer_ratio()
print(x)
3602879701896397
print(y)
36028797018963968
x/y
0.1
となり、0.1の分母と分数を取得できる
これを用いれば、さきほどのinterval型への代入問題を解決できそうだが
sup = pn.roundfloat.rddiv(x, y, pn.roundmode.up)
inf = pn.roundfloat.rddiv(x, y, pn.roundmode.down)
format(sup,'.17g')
'0.10000000000000001'
format(inf,'.17g')
'0.10000000000000001'
となり、0.1を包含するように計算できない。
なぜならば a=0.1 とした時点で 0.1 から、0.1よりずれた小数に近似されるため、as_interger_ratio で取得できる分母と分子は 0.1よりずれた小数に合うように取得されるためである。
誤解回避
0.1...となる分数を取得することは可能だが
厳密に0.1になるような分数はas_integer_ratioでは取得できない
(このように厳密に計算するとき以外は大した問題にならないと思われる)
メモ:twosumアルゴリズムの研究
twosumアルゴリズムの研究
from pint import roundfloat as rf
from pint import roundmode as rdm
a=0.1
b=0.2
誤差が発生するパターン
rf.rdadd(a, b, rdm.nearest)
0.30000000000000004
最近点丸め
rf.rdadd(a, b, rdm.up)
0.30000000000000004
上方向丸め。最近点丸めと一緒なので、最近点丸めでは上方向に丸め込まれた
rf.rdadd(a, b, rdm.down)
0.3
下方向丸め。最近点丸めと別方向に曲がったことがわかる
rf.twosum(a, b)
(0.30000000000000004, -2.7755575615628914e-17)
twosumを実行してみると、最近点丸めの際の、真値(2進数表現できない)との誤差の方向が検出できているのがわかる。
詳しく見ていこう
x = a + b
print(x)
0.30000000000000004
普通の足し算だ
tmp = x - b
print(tmp)
0.10000000000000003
tmp = x - b = (a + b) - b = a
であるのでこれはaとなるはずである。
だが、ならない
xは最近点丸めで計算されているため、bを引き算したときに、aに誤差sが加算されている。
tmp = a + 誤差s という形になっている
y = a - tmp
print(y)
-2.7755575615628914e-17
y = a - tmp = a - (a + 誤差s) = 誤差s という式になっている
pythonの区間演算ライブラリ進捗2015-11-16
成果
github.com
ソースコードが長くなりそうだったので分離した。
本当は一つのファイルで行きたかったけど、管理や読み込みなどを考えたら分離せざるを得なくなった。
- list型を継承したinterval型も入る便利な配列を作る。
に該当するmcmatrixクラス - 標準のmathを継承し、interval型対応と精度保証された値を出力するmathクラスを作る。 の該当する予定のvmathクラスを作成
反省
独自の配列は結局、list型を継承して作ることにしたので、パフォーマンスは期待できない。
今後改良していく予定はもちろんあるが、numpyのndarrayと連携しようとすると、死にそう。
今後の方針
とりあえず、mcmatrixクラスに内積(dot)や転置、zerosやonesやらの実装をしないとまともに使えないので、さっさと実装する。
pythonの区間演算ライブラリ進捗2015-10-25
成果
反省
PEP8対応のため、手動で途中まで直したが、autopep8コマンドで十分だった。(勉強にはなったけど)
今後の方針
- list型を継承したinterval型も入る便利な配列を作る。
- 標準のmathを継承し、interval型対応と精度保証された値を出力するmathクラスを作る。
pythonの区間演算ライブラリ進捗2015-10-13
成果
numpyのndarray型に対して、numpyのndarrayの加減乗っぽく計算する、仮想丸め方向付き計算クラスroundarrayを実装
反省
- ndarrayの除算っぽい計算が実装できなかった(数値計算にとって必要なのかもよくわからない……)
- double型の行列に対してtwosum,twoproductで誤差を拾おうとしているが、たぶん無理そう(そのうち別記事でまとめる予定)。
今後の方針
最初に掲げていた「interval型で簡単にplotしたい」は怪しくなってきた。
今後はnumpy依存部分と、標準ライブラリのみの部分を分けて開発をしていく。
numpy依存部分
引き続き「interval型で簡単にplotしたい」を実現するためにmatplotlibやらndarrayと戦っていく。
(でも、効率の良いinterval型の配列は中心値形式を利用しなければできそうにないので、時間がかかりそう)
しかし、使えるものは何でも使っていく予定
標準ライブラリのみ部分
「標準ライブラリでどこまで精度保証できるか?」を追求していく。
しかし、アルゴリズム以外の部分での高速化(BLASを使うetc...)はあまり求めないでいく。
ここまで来たら、「pythonのみで実装する」も方針に入れていこうと思う。
少し、コード周りを整理したら、githubで公開する予定。
ブラウザで動くIDE codenvyで遊ぶ
なんか他に書くこと溜まってるような気がするけど、便利そうだったので書いておく。
数年前にニュースサイトで目にしたのでお気に入りに放り込んでおいてたが、全く使っていなかった。
お気に入りの整理をしている最中に目についたので、アクセスしてみたら、超進化していたので、使ってみることにした。
いくつかの言語(Java,C++,ruby,python...)に対応しており、自分のPCに一切影響を与えることなく、開発が可能。
python3.4を使ってみたが、標準ライブラリーは一通り揃っており、パッと思いついたことをブラウザ一つで試せるのは楽しい。
しかし、numpyなどの追加ライブラリーは無かったので、本格的な運用には向かないようだ。
とりあえず、今開発しているpythonのライブラリも一部は動作した。
標準ライブラリー郡のみで実装すると、こういう場所でも動作するので、いつかやってはみたいなぁ。
コードはここに置いておきます。
github.com
