数値解析日記

数値解析関係のネタ帳。C++,pythonなど

Unite 2017に行ってきたよ

Unite 2017に行ってきたよ

見たセッションとかメモとか

Day1

「魔法使いプリキュア!」EDでのUnity映像表現の詳細解説

非公開なのでおおまか。 導入前後での作業フローの比較など。グラフィック寄りの話で素敵なシェーダーの話など。途中で抜けてしまったが、最後の人がおもしろかったらしい(友人談)。
展示ブースでもEDをVR体験。友人が下から覗こうとしたが、当然健全なので、視点は不必要に下がらないようになっていた。治安が良い。

Unity UI最適化ガイド 〜ベストプラクティスと新機能

講演動画、資料リンクから。
プロファイラーなどの使い方を詳しく実践してくれる。学習段階で言うと、Unityでゲームを作り始める前に見たい。

バグを殲滅!Unityにおける実践テスト手法

講演動画、資料はリンクから。
テストツール全般に対してのお話。やっぱUIのデバッグはきついか……。

シェーダープログラミング入門!カスタムシェーダー、作るで!

講演動画、資料はリンクから。
シェーダーも正直魔境だと思って食わず嫌いをしていたが、これ見てるとなんとかできそうな気もする。もう一回見直さねば。

Day2

実演!Timeline+Cinemachineによるカットシーン作り

講演、資料はリンクから。
そのうち搭載されるであろうTimeLineの紹介。今ライブラリ作ってるけど、これできたら用済みになる気がしてきた。南無。

最新 PS VR コンテンツ制作事例紹介 with Unity

非公開なので、おおまかに。 「なぜDWがマシュVRにUnityを使ったのか?」「カヤック傷物語VRで使った表現技法」など説明。これ見ると、両方やりたくなっちゃう。

『Shadowverse開発事例』 ~美麗カードが動く!制作テクニックのすべて~

講演、資料はリンクから。
みんな大好きCygames。
プレミアムカードのエフェクトに使っているシェーダーのお話など。今回のUniteだけでも、シェーダー自作のお話が3, 4つあった気がする。やっぱ軽くて良いんだね。

展示ブース

一覧はここ

日本マイクロソフト株式会社

サイクロップス先輩みたいな写真が撮れる。職業に「VR」とかあったけど、それって職業なんでしょうか?

株式会社Cygames

AGE AGE VR。正直楽しい。
少しルールに慣れがいるが、ゲーセンにあったらやっちゃいそう。

イリュージョン

身体検査をした。やましい気持ちはない。
ブースのお兄さん「誰も見てないから安心して心置きなく検査してね」

魔法使いプリキュア!(東映アニメーション株式会社)

目の前で踊ってくれるのはなんだかんだドキっとしてしまう。

こんな感じ。来年はもういけないかなー。

n=5000兆でヤコブ・ベルヌーイのeの定義は使えるか?

n=5000兆でヤコブ・ベルヌーイのeの定義は使えるか?

https://twitter.com/nrmti/status/868980739047530496 の検証
見づらかったら、ここ見てね 素のコード jupyter ipynb

言語はPythonです。

import decimal
decimal.Context()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])

5000兆円欲しい

want = decimal.Decimal(5000 * 10 ** 12)
want
Decimal('5000000000000000')

不要かもだけど

one = decimal.Decimal('1')
one
Decimal('1')

wikipediaの収束数列による定義を参照

e = (one + one / want) ** want
e
Decimal('2.718281828459044963532104625')
import math
math.e
2.718281828459045

うーん、まあまあか

おまけ: doubleでやると

double_want = 5000 * 10 ** 12
e = (1 + 1 / double_want) ** double_want
e
3.0350352065492636

ダメダメですね

武蔵野アブラ学会の美味しいトッピング

武蔵野アブラ学会の美味しいトッピング

 一番槍再び。 この記事はみす裏アドベントカレンダーの一日目です。
テーマなんぞ与えられてないので好きに書きます。

 武蔵野アブラ學会というお店をご存じでしょうか? 公式サイト
 油そばの専門店でここ数年虜になっているので紹介したいと思います。
 油そばの説明とかは省くので各自ぐぐってください。 (学会のQ&Aにも載ってます)

美味しいところ

濃いです。 やたらと濃い味が好きな人はすぐに虜になります。

基本メニュー

これにトッピングや調味料をプラスしていく。

魅惑のトッピングたち

海苔

 私が一番オススメするトッピングは海苔。

  • パリパリの状態で麺に巻いて食べる
  • ひたひたにした状態でご飯に巻き込む
  • 最後までとっておいて、残ったタレを全てすくうように食べる

など可能性に満ちている。 しかも、刻み海苔変更が可能だ。刻み海苔は香りが高く、海苔とタレがより高度に絡むようになる。

f:id:o108minmin:20150730172326j:plain

 一部だけ刻みになったのりたまも両方楽しめるので良い。 ちなみに海苔を追加して、すべて刻み海苔に変更するとこうなる。

f:id:o108minmin:20161028110142j:plain

最高。
ただし、刻み海苔はひたひたになると風味が失われるので、さっさと食べるといい感じ。

味玉

f:id:o108minmin:20160905114059j:plain

 濃いめの味玉は黄身が濃縮されていて、生の状態よりもアクセントを強くつけてくれる。 ご飯のお供にも最適。

生卵

f:id:o108minmin:20161109111414j:plain

 白身と合わさって、タレが薄くなるのでは? と懸念した私を待っていたのは卵黄だった。 卵を溶かして和えると、タレが少し水っぽくなって食べやすくなるので食欲が無いときに向いてるかもしれない。

メンマ

f:id:o108minmin:20150509104758j:plain

 メンマを増すと、タレを残さず食べるためのワイパー代わりになるので、残さず食べられる。メンマの甘さで塩辛さとの緩急がつくので飽きやすい人にも向いている。

肉マシ

f:id:o108minmin:20140520110855j:plain

 學会のチャーシューは固くもなく柔らかすぎないため、いくらでも食べれる。
だが、そんなにチャーシューが好きなら、肉そばにすればもっと幸せになれるぞ。

ねぎゴマ攻め

f:id:o108minmin:20150904132959j:plain

 白ねぎが大量に追加される。大盛り程度だと、しなしなになりきらないので、ダブル以上での使用をオススメする。
ゴマ? あいつは良いやつだったよ。ゴマの風味が欲しいならラー油をかけよう。

ニンニク増し

f:id:o108minmin:20150120111220j:plain

 店内で刻まれたとてもフレッシュな生ニンニクがトッピングできる。 とても香りが強く、大盛り程度だと味が大幅に崩れることもあるので注意。 しかし、特濃油そばはそれに負けない 強さ を持っているので、必ずトッピングすること。 ちなみに、昨日特濃ニンニク増しで食べたら、不調気味だったのが回復したので、これからの風邪予防にぜひ。

f:id:o108minmin:20161129114354j:plain

 ほかにもキムチやマヨネーズなどのトッピングがあるが、研究できてないので今後の課題とする。

調味料

 食べるときには、食卓の脇においてある調味料も欠かせない。
これらも紹介していこう。

 油そばの定番。かけすぎるとタレが薄くなってしまうが、食べやすくなる。 食べるのがつらくなってきたら、酢を増してサッパリと食べていくのが良い。

胡麻ラー油

 ゴマ油の香りが強く押し出されたラー油。これと酢の相性は抜群なので是非やってほしい。

山椒ラー油

 山椒の香りがつけられたラー油。胡麻ラー油よりも控えめで爽やかさが増す。胡麻ラー油をかけずに山椒ラー油だけ使うのもまた一興。

以下のものは食べる分だけかけながら食べると、味の変化をより効果的に楽しめる。

コショウ

 あらびきのものが銀色の容器に入っている。酢とラー油をかけたあとの最初の変化をつけたいときにかけてみて。

山椒

 ミルでひいてかけるので、舌にビリビリとくる感触が濃い味と絡んで心地よい爽やかさを演出してくれる。 山椒ラー油と合わせるのが最高。

たまり

 ラー油を作ったときのスパイスらしい。 辛さは控えめだが、麺とよく絡めると辛さが増す。 香りが独特なので最初は加減しながら入れよう。

最後に

 いろいろ言ったが、最初に食べるなら「海苔味玉油そば大盛り + ご飯」がオススメ。
次からは特濃なり、自分だけのトッピングの組み合わせを模索していこう。

f:id:o108minmin:20160706110828j:plain

良い油そばライフを

nextafter関数の使い方

nextafter関数の使い方

import numpy as np

numpyのnextafter関数を使っていたらこんな事態に遭遇した。

nextafter(x1, x2[, out])

Return the next floating-point value after x1 towards x2, element-wise.

np.nextafter(2, 1)
1.9999999999999998

え? x2に+の値入れたのに、下方向に行った!?   慌てて、リファレンスを見なおしてみると、  

after x1 towards x2,

となっている。このtoward x2がポイントで、「x2の方向に向かって」とある。
これは、「x1からx2への方向」なので、上の例だと「2から1の方向」なので、マイナス方向になる。
よって、正しく2の隣の浮動小数点数を見つけたかったら

np.nextafter(2, 2+1)
2.0000000000000004

としなければならない。これを一般化して、succ関数とすると以下のようになる。

def succ(x):
    return np.nextafter(x, x+1)

追記: 上記の例は大きなxに対しては機能しない場合がある。後にnextafterを使った完璧なsuccを公開する予定です。

succ(2)
2.0000000000000004

となる。
ちなみにStackOverFlowでも同様のミスをしている人を発見。

JuliaのFloatとRoundingの奇妙な挙動[追記]

JuliaのFloatとRoundingの奇妙な挙動

with_rounding(Float64, RoundNearest) do
    exp(-0.999)
end
0.36824750461366296
with_rounding(Float64, RoundUp) do
    exp(-0.999)
end
0.3682475046136629
with_rounding(Float64, RoundDown) do
    exp(-0.999)
end
0.3682475046136629

本来、 Down <= Nearest <= Upとならなければならないところが
見た目上は Up < Nearest となってしまっている。


追記: Up <= Nearest を確認してみると成立しているため、計算上の問題はない 。表記がおかしいだけで、内部的には影響は小さいようだ


もちろん、Roundingは四則演算と平方根のみなのだが、結果が逆になったのは興味深い。
JuliaのRoundingのドキュメントを参照すると

Waring

Rounding is generally only correct for basic arithmetic functions (+(), -(), *(), /() and sqrt()) and type conversion operations. Many other functions assume the default RoundNearest mode is set, and can give erroneous results when operating under other rounding modes.

とあるので、expはおかしくなってしまったんだろう。issueをさらに調べると
exp(1) != e #3370  

Since exp(1) just calls openlibm's exp function, this appears to be an openlibm issue.

とあるので、openlibm側の問題のようだ。
当たり前だがexpにはRoundingは効かず、場合によっては妙な挙動を起こすことがわかる。
なお、これは利用者側が気をつけなければならない仕様であるので、Julia側に落ち度はない。
with_roundingやset_rounding内では、四則演算と平方根以外は利用しないほうが良いだろう。
(できれば、四則演算と平方根以外は最近点丸めで計算してくれるとありがたいんだけど…… )

なお、mpfrを利用するBigFloatではこのような問題は起こらない

with_rounding(BigFloat, RoundNearest) do
    exp(BigFloat(-0.999))
end
3.682475046136629215368118018604961977815208270355409127685022306741116356975955e-01
with_rounding(BigFloat, RoundUp) do
    exp(BigFloat(-0.999))
end
3.682475046136629215368118018604961977815208270355409127685022306741116356975998e-01
with_rounding(BigFloat, RoundDown) do
    exp(BigFloat(-0.999))
end
3.682475046136629215368118018604961977815208270355409127685022306741116356975955e-01

見づらいんで、大小比較すると

with_rounding(BigFloat, RoundNearest) do
    global near
    near = exp(BigFloat(-0.999))
end

with_rounding(BigFloat, RoundUp) do
    global up
    up = exp(BigFloat(-0.999))
end

with_rounding(BigFloat, RoundDown) do
    global down
    down = exp(BigFloat(-0.999))
end
3.682475046136629215368118018604961977815208270355409127685022306741116356975955e-01
down <= near <= up
true

となるのでmpfr上では大丈夫
なお、この現象はJulia 0.4.5で確認し、windows 10(64bit)、OS X 10.11.5で確認した

numpyを使って、行列に対するtwosum, twoproductを高速化する

 柏木先生の考案された改良型twosum, twoproductを行列に適用し、それをnumpy流にベクトル化し、高速化した。
コードと結果は以下
github.com
github.com 環境は

  1. Mac Book Pro 2012 mid (ただし、メモリを16GBに増やしている)
  2. anaconda 4.0.8 (以下anaconda内のパッケージ)
     * python 3.5.1  
     * scipy 0.17.0 np110py35_4  
     * numpy 1.10.4 py35_2  
     * mkl 11.3.3  
    

結果は

%timeit twosum(a, b)
10 loops, best of 3: 80.7 ms per loop
%timeit twosum_old(a, b)
1 loop, best of 3: 2 s per loop
%timeit twoproduct(a, b)
10 loops, best of 3: 76.8 ms per loop
%timeit twoproduct_old(a, b)
1 loop, best of 3: 4.22 s per loop

 特に特殊な手法は用いていないが、ベクトル化するだけで大幅に高速化した。

分数区間を用いた、浮動小数点数の区間包含についてのアイデア

分数区間を用いた、浮動小数点数の区間包含について

python 3.5.1で動作を確認
区間演算の実装について(2)python上解決する場合のアイデア
いろいろと検証が足りていないが公開はする。
jupyter notebook用のコードもアップロードした。

import pint as pn
from pint import roundfloat as rf
from pint import roundmode as rdm
import fractions
def frac_interval_proto(a):
    aH, aL = rf.split(a)
    if aL ==0:
        answer = pn.interval(a)
    else:
        aS = rf.succ(a)
        aP = rf.pred(a)
        aS_c, aS_p = aS.as_integer_ratio()
        aP_c, aP_p = aP.as_integer_ratio()
        a_c, a_p = a.as_integer_ratio()
        bS_c = a_c * aS_p + aS_c * a_p
        bS_p = a_p * aS_p * 2
        bP_c = a_c  * aP_p + aP_c * a_p
        bP_p = a_p * aP_p * 2
        answer = pn.interval(0.)
        answer.inf = fractions.Fraction(bP_c, bP_p)
        answer.sup = fractions.Fraction(bS_c, bS_p)
    return answer

今回は、代入演算子を用いて区間型を生成した際に、初期値を包含しない区間が発生するのを防ぐ分数区間を生成する。
すなわち

a = 0.1
itv_a = pn.interval(a)
format(itv_a.inf, '.17g')
'0.10000000000000001'
format(itv_a.sup, '.17g')
'0.10000000000000001'

のような区間が発生することを防ぐ。
なお、今回以外の解決方法としては

itv_a.inf = rf.pred(a)
itv_a.sup = rf.succ(a)
format(itv_a.inf, '.17g')
'0.099999999999999992'
format(itv_a.sup,'.17g')
'0.10000000000000002'

のようにsucc predを用いることで解決可能である。
しかし、今回のアルゴリズムのほうが、より区間幅を縮小できる。

まず、上位ビットと下位ビットを分離する

aH, aL = rf.split(a)

もし、下位ビットが存在する場合、else文以降の通りにする。
(下位ビットが存在しない浮動小数点数で点区間を生成した場合、元の初期値を包含する区間になる可能性が高い
元の初期値の下位ビットが全て1で、近似の際にくりあがりで下位ビットが消去された場合、おかしくなるかも)

if aL ==0:
    print(True)
    answer = pn.interval(a)
else:
    print(False)
False

以下はelseの続き

次にaをsucc predし、as_intger_ratio()で近似分数にする

aS = rf.succ(a)
aP = rf.pred(a)
aS_c, aS_p = aS.as_integer_ratio()
aP_c, aP_p = aP.as_integer_ratio()
a_c, a_p = a.as_integer_ratio()

次に、aとsuccしたaの中点と、aとpredしたaの中点を求める。

bS_c = a_c * aS_p + aS_c * a_p
bS_p = a_p * aS_p * 2
bP_c = a_c  * aP_p + aP_c * a_p
bP_p = a_p * aP_p * 2

こうして得た、分数bSとbPとaは以下の性質を満たす

浮動小数点数として評価した場合、同値である。

2016-02-01 例外がたくさん発生した
正しくは bP_c / bP_p <= a <= bS_c / bS_p

bS_c / bS_p == a == bP_c / bP_p
True

分数として評価するとbP <= a <= bSが成立する
(以下は、分母と通分するためbS_pやら、a_pがかけられている。今回は bP < a < bSとなっていることをわかりやすくするため、等号の場合は考慮していない)

a <= bS

a_c * bS_p < bS_c * a_p 
True

bP <= a

bP_c * a_p < a_c * bP_p
True

よって、分数で評価すると大小関係が存在するが、浮動小数点数で評価すると値が一致する分数区間が得られた

なお、succ predで得られたitv_aと比較すると

rf.pred(a) < bP_c / bP_p
True
bS_c / bS_p < rf.succ(a)
True

となり、succ predよりも狭い区間幅が得られた。

answer = pn.interval(0.)
answer.inf = fractions.Fraction(bP_c, bP_p)
answer.sup = fractions.Fraction(bS_c, bS_p)
print(answer)
[14411518807585587/144115188075855872,14411518807585589/144115188075855872]

今回は

  • 引数が正の場合しか考慮していない
  • float.as_intger_ratio()で得られる分数が、最近点丸めで近似されたaと厳密に同じかどうか考慮していない
  • 分数演算の有効桁数に由来する計算誤差を考慮していない(分母、分子ともに巨大な整数になりがちである)
  • その後の分数の演算方法を考慮していない(作ったは良いが、使い方を考えていない)

などの問題がある
最後の分数の演算に関しては以下のようにするアイデアがあるが、検証はしていない

  • 分数 * 分数: 分数演算をしたあと、丸め制御で適切な浮動小数点数にする
  • 分数 * 浮動小数点数(下位ビットがない): 分子 * 小数を計算後、丸め制御で適切な浮動小数点数にする
  • 分数 * 浮動小数点数(下位ビットがある): 小数を分数に変換したあと、分数 * 分数に帰着させる