読者です 読者をやめる 読者になる 読者になる

pythonで非線形最小二乗法

scipy.optimize.leastsq()の使い方で苦戦したのでメモ

LM(Levenberg-Marquardt)法を実装するのが面倒だったので既存のleastsq(LM法で非線形最小二乗してくれるやつ)を使うことにした.
というわけで、良くわからないまま適当に以下のようなコードを書いた(研究用のコードを載せるのはアレなので適当に書き換えてる)

"""  
minimize y - a x^2 + b x + c 
params : a,b,c
array_of_Point2d : (x,y) data
"""
def cost_func(param,array_of_Point2d):
    a,b,c = params
    sum_of_err = 0.0
    for point in array_of_Point2d:
        x,y = point.x,point.y
        sum_of_err += (y - (a*x**2 + b*x + c))**2
    return cost_func

params = [0,0,0]
array_of_Point2d = [Point2d(i,i) for i in range(10)]
leastsq(cost_func,x0=params,args=array_of_Point2d)

で、実行してみると案の定エラーを吐かれた

Improper input: N=6 must not exceed M=2

ググる
→ http://stackoverflow.com/questions/21592348/scipy-leastsq-fitting-a-square-grid-to-experimental-points-in-2d

The function passed to leastsq (e.g. err_func) should return an array of values of the same shape as xeand ye

ようするにleastsqに渡す関数の戻り値は、関数に渡す引数と同じ形じゃないとアカンとのこと.
x,yデータの数だけ誤差を計算し、その誤差を全て返すのが正しい形っぽい

修正版

"""
minimize y - a x^2 + b x + c 
params : a,b,c
array_of_Point2d : (x,y) data
"""
def cost_func(param,array_of_Point2d):
    a,b,c = params
    errs = []
    for point in array_of_Point2d:
        x,y = point.x,point.y
        errs.append(y - (a*x**2 + b*x + c))
    return errs

params = [0,0,0]
array_of_Point2d = [Point2d(i,i) for i in range(10)]
leastsq(cost_func,x0=params,args=array_of_Point2d

というか、僕がミスったようなことは全部リファレンスに書いてあったscipy.optimize.leastsq - SciPy v0.18.1 Reference
英語だからって面倒くさがらずにちゃんと読もうな!

追記

どうやら使っているScipyが古かった模様 scipy 0.17.1以降はleast_squares()という関数があって、こいつを使うよう推奨されている
least squares - SciPy: leastsq vs least_squares - Stack Overflow

least_squares()ではキーワード引数methodを指定することで以下の三つのアルゴリズムを使い分けることができる
- ‘trf’ : Trust Region Reflective algorithm
- ‘dogbox’ : dogleg algorithm
- ‘lm’ : Levenberg-Marquardt

ついでに言えばleast_squaresの方は解に上下限を付けることができ、解(x,y,z)について
bounds = ([x_l_lim,y_l_lim,z_l_lim],[x_u_lim,y_u_lim,z_u_lim])
のように解の下限リストと上限リストをタプルにして渡せばよい
(注:LM法を使う場合は設定できない)

ディープではないcopy.deepcopyの話

コピーについて考える

最近、強化学習の勉強の一環としてSuttonBookのLispコードをpythonに書き直してうごかしています。
Lispオブジェクト指向ではないけれど、(Common Lisp Object Systemというオブジェクト指向システムを持っているみたい。勘違いしてました(;´・ω・) )どうせならと思って移植の際にクラス化したりしていたら、グリーディな行動選択をさせる際に オブジェクトをコピーする必要が出てきて、安易にcopy.deepcopyを使ったらクッソ遅くなってしまった。
結局、新しいインスタンス生成して必要なメンバのみコピーするように書き換えたらサクサク動くようになりました。

(そういえば研究で書いてた重いコードも安易にdeepcopyつかっていたような…)

そんなわけで、pythonでオブジェクトや変数をコピーするときどういう風にコピーするのが一番いいのか気になって調べたことをまとめる。

実際どれくらい遅いのか

まず、deepcopyを使う方法と使わない方法でどれくらいcopyの速さが違うのか確認しておく


例1 単純なリスト

長さ1000のリストをlist()で複製した場合とdeepcopyした場合の処理時間を計測する。 処理回数は10000とした。


コード

from timeit import timeit
num = 10000
print timeit("list(test_list)",setup='test_list = [i for i in xrange(1000)]',number=num)
print timeit("copy.deepcopy(test_list)",setup='import copy ;test_list = [i for i in xrange(1000)]',number = num)

結果

0.272120951146
1793.11939573

例2 適当なクラスをコピーする

以下に定義したtest_copyクラスをコンストラクタで複製した場合とdeepcopyで複製した場合の処理時間を計測する。 処理回数は同じく10000回


コード
test_copy.py

class test_copy(object):
    def __init__(self,test_obj=None):
        if test_obj is None:
            self.a = 1
            self.l = [1,2,3,4,5,6,7,8,9]
            self.text ="aiueo"
        elif isinstance(test_obj,test_copy):
            self.a = test_obj.a
            self.l = list(test_obj.l)
            self.text = str(test_obj.text)
num = 10000
print timeit("test_copy(test_obj)",setup='from test_copy import test_copy; test_obj = test_copy()',number=num)
print timeit("copy.deepcopy(test_obj)",setup='from test_copy import test_copy;import copy ;test_obj = test_copy()',number = num)

結果

0.814853623342
40.5248817692

というわけで、deepcopyは確かに処理が重いことが確認できた。

pythonレファレンスによると、deepcopyの場合は単にコピーするのではなく、コピー済みオブジェクトを記録するなど、色々と処理を行っているようだ。

実際copy.pyの中身をみるとメモ化や型チェックで100行ほどのコードになっている。
またネストしたリストのような複雑なオブジェクトをコピーするためにcopy.deepcopyを再帰的に呼び出すような処理もある。(deepcopy内の_reconstruct()は内部でdeepcopyを呼んでいる) deepcopyの処理の遅さはこの再帰呼び出しによるものかと思われる。

def deepcopy(x, memo=None, _nil=[]):
    """Deep copy operation on arbitrary Python objects.

    See the module's __doc__ string for more info.
    """

    if memo is None:
        memo = {}

    d = id(x)
    y = memo.get(d, _nil)
    if y is not _nil:
        return y

    cls = type(x)

    copier = _deepcopy_dispatch.get(cls)
    if copier:
        y = copier(x, memo)
    else:
        try:
            issc = issubclass(cls, type)
        except TypeError: # cls is not a class (old Boost; see SF #502085)
            issc = 0
        if issc:
            y = _deepcopy_atomic(x, memo)
        else:
            copier = getattr(x, "__deepcopy__", None)
            if copier:
                y = copier(memo)
            else:
                reductor = dispatch_table.get(cls)
                if reductor:
                    rv = reductor(x)
                else:
                    reductor = getattr(x, "__reduce_ex__", None)
                    if reductor:
                        rv = reductor(2)
                    else:
                        reductor = getattr(x, "__reduce__", None)
                        if reductor:
                            rv = reductor()
                        else:
                            raise Error(
                                "un(deep)copyable object of type %s" % cls)
                y = _reconstruct(x, rv, 1, memo)

    memo[d] = y
    _keep_alive(x, memo) # Make sure x lives at least as long as d
    return y

またdeepcopyの処理速度について調べていたらこんな記事が
python - copy.deepcopy vs pickle - Stack Overflow

cPickleをつかってpickle化したオブジェクトを再びloadすることでcopy.deepcopy()を使うより高速に ディープコピーすることができるんだが、copy.deepcopy()を使う利点てあるの??

という質問なんだけど、そもそもpickle化したオブジェクトをloadすることでコピーするというテクニックを初めて知った。

で、回答としては

Pickleによるコピーがdeepcopyに比べて速度が速いが、deepcopyより一般性に欠ける手法である

とのこと。

確かにpickle化可能な型には制限があるし、(11.1. pickle — Python オブジェクトの直列化 — Python 2.7.x ドキュメント) copy対象のクラスが改変されたらバグになりうるから、使う場合はassertいれるとかコメント入れるとか、とにかくきちんと対応したほうがよさそう。 pickle、unpickleだけでコピーできる手軽さは魅力的ではあるが。

まとめ

コピー対象が比較的単純な場合はcPickleつかう。複雑なオブジェクトを正確にコピーしたい場合はcopy.deepcopyを使う。速度が必要な場合は専用にdeepcopyを作るのがよさそう。

AtCoderやってみた

今回はAtCoderの問題を練習がてら解いてみたよ、という日記です。

部分点しか取れなくてもやもやしてるのでここにあげてみます。

問題はCODE FESTIVAL 2016 Final (Parallel)のB問題、詳細は以下

Exactry N Points

  • N問の問題が与えられる

  • i番目の問題を解くとi点もらえる

  • 問題番号が大きいほど難易度が高い

  • N問のうちいくつかの問題を解き、合計得点がN点になるようにしたい。この時できるだけ簡単な問題を解いてN点にせよ(解答する問題の問題番号の最大値が小さくなるようにせよ)


最初に考えたのはNを2分割していって解を求めると言う分割統治的なアプローチだったんですが、これだと分割していった結果同じ数が現れたときの処理が大変だと思い断念。

とりあえずNが簡単な場合に、和がNになる番号の組み合わせで重複のないものを考えて、N=kの解からN=k+1の解を構成できないか考えます。

f:id:unthocyanidin:20161210091805g:plain:h450

線で結ばれたノードは上位ノードに以下の操作をすることで下位ノードが構成できることを示しています。
- ノードの要素の一つに1を足す
- 要素1を追加する
またピンク色のノードは問題番号の最大値が最も小さいノードを示しています。(つまりピンクのノードがN=kのときの解)

図から、帰納的に解を構成できそうだとわかりました。次にピンクのノードだけ取り出して、解の構成法に規則性がないか見てみます。

f:id:unthocyanidin:20161210092651g:plain:h450

赤字で示したのはN=kの解からN=k+1の解を構成する際に"1"増えた要素です。 図をみるとツリーの層を下るに従って赤字が右から左にずれているのがわかります。これが解の規則性です。 規則性のおかげでN=kのときの解が既知で、N=kの解を構成するときにどの要素に1を加えたかを記録しておけばN=k+1の解も構成することができます。

以上の方法で解答を提出すると部分点がGetできました。しかしNが大きい場合にタイムアウトしてしまう...

せっかく規則性がわかったのだから、そこから一般項を求めることができれば、あるいは一般項とまではいかなくともボトムアップ的な解法を回避できればNが大きい場合でもうまく行くはずです。

そこで、この図を要素数ごとにグループ分けしてみます。 f:id:unthocyanidin:20161210100201p:plain:h450 すると第n群の1番目の組み合わせは{1,2,3,...,n}であるとわかります。このことから第n群のi番目の組み合わせは{1,2,3,...,n}のうち大きいほうから要素をi個選び、それぞれ1を足したものとなります。

というわけで、解法として以下の3ステップで解を構成することができます。 1. 与えられたNが第何群か計算する 2. Nの属する群の初項を求める 3. 初項にたいして適当な数だけ1を加える

この方法では1から解を構成するのではなくN群の属する初項を求めてから解を作るため最初の方法に比べて計算量が少なくなります。

とりあえずこの方法でもう一回、解答なげてみますかね。うまく行けば満点もらえるかと。

どくおぷと!

おはようございます、こんにちは、こんばんは。 久々の更新です。

さて、皆さんはpythonスクリプトを書くときコマンドライン引数のパースはどのようにしていますか?

  1. sys.argv[1:]を自力でパース
    コマンドが簡単で、特に使いまわす予定のないスクリプトならそれでもいいですが、自分以外の人がスクリプトを使う場合はこのやり方は嫌われますね。

  2. argparseやgetopt,optparseモジュールを使う
    なるほど既存のモジュールを使うのは賢い選択です。このやり方ならスクリプトの引数を間違えた場合にusageを表示させるのも楽でしょう
    *getopt,optparseモジュールは推奨されていないらしいです。
    argparseでコマンドライン引数の解析 | Python Snippets

ちなみに僕はJupyter NotebookやVSCodeの上でスクリプトを実行することが多かったので今まで2の存在をしらなかったです(恥ずかしい)
ですが先日コマンドライン引数をパースする機会に恵まれまして、その辺のモジュールを調べていたらよりも簡単にパースしてくれるモジュールがありました。 Pythonのコマンドライン引数処理の決定版 docopt (あとJuliaに移植したよ) - りんごがでている

その名も docopt

docoptの詳しい説明は公式をどうぞ

docopt—language for description of command-line interfaces

docoptの画期的なところは__doc__スクリプトの使用方法を書いておけば、それをもとにパースしてくれるという点です。 なので使い方も超簡単です!
基本的には以下の2stepでパースできます
1. スクリプト__doc__を書く
2. argv = docopt(__doc__)

これだけです。

1. スクリプト__doc__を書く
2. argv = docopt(__doc__)

だけです

あんまり簡単なのでつい2回書いてしまいました。それにしても簡単だ。

もう少し説明をするとdocopt()の第一引数はスクリプトの使用方法を記述した文字列、第二引数にコマンドライン引数のリストargvを渡すことができます。第二引数を指定しない場合sys.argv[1:]が使用されます。
docopt()には他にもhelp,version,options_firstといった名前付き引数があるみたいですが、この辺の使い方はおいおい調べます。

そしてdocoptが正常にコマンドをパースできた場合、受け取ったコマンドにキーワードが含まれるか否か表す辞書を返します。

具体的な例で説明すると。(以下公式より抜粋)

"""Naval Fate.

Usage:
  naval_fate.py ship new <name>...
  naval_fate.py ship <name> move <x> <y> [--speed=<kn>]
  naval_fate.py ship shoot <x> <y>
  naval_fate.py mine (set|remove) <x> <y> [--moored|--drifting]
  naval_fate.py -h | --help
  naval_fate.py --version

Options:
  -h --help     Show this screen.
  --version     Show version.
  --speed=<kn>  Speed in knots [default: 10].
  --moored      Moored (anchored) mine.
  --drifting    Drifting mine.
"""

このような__doc__がかかれたスクリプトnaval_fase.pyがあった時コマンドプロンプトから
"python naval_fate.py ship new Prinz"と入力するとdocoptはsys.argv[1:](すなわちship new Prinz)を解析して,

{
  "--drifting": false, 
  "--help": false, 
  "--moored": false, 
  "--speed": "10", 
  "--version": false, 
  "<name>": [
    "Prinz"
  ], 
  "<x>": null, 
  "<y>": null, 
  "mine": false, 
  "move": false, 
  "new": true, 
  "remove": false, 
  "set": false, 
  "ship": true, 
  "shoot": false
}

という辞書を返します。(公式サイトのTry docopt in your browserに行くとお試しできるよ!)

一行目にはスクリプトの概要を記述します。usageに具体的な文法、optionsにオプション引数を記述します.
usageの各行は <スクリプト名> <スペース区切りのコマンドキーワード> という形式になっています.
上の例に合わせて<スクリプト名>と表記しましたがこの部分は各コマンドで共通であれば実はなんでもよいみたいです.
たとえば

Usage:
  _ ship new <name>...
  _ ship <name> move <x> <y> [--speed=<kn>]
  _ ship shoot <x> <y>
  _ mine (set|remove) <x> <y> [--moored|--drifting]
  _ -h | --help
  _ --version

このようにスクリプトに関係ない記号(_)でもエラーは出ません。 これはdocoptメソッドがデフォルトではsys.argv[1:]をパースするため,もともとスクリプト名がパースされているからと思われます。 プログラム内部でスクリプトを実行する場合なんかに使えるかもですね。

リサンプリング

学業がひと段落したので、久々に更新。

ちょっとパーティクルフィルタを実装する必要があって、pythonでコードを書いているところです。 今回はその時に調べたことのメモ。 パーティクルフィルタの説明はgoogle大先生にお伺いするといろいろ出てきますが、MISTのチュートリアルが実装も書いてあって わかりやすいかと思います。

Tutorial/Practice0 – MIST Project

パーティクルフィルタでは各パーティクルに尤度に基づいた重みをつけ、イテレートのたびに重みに従ってパーティクルをリサンプリングします。

重みに従ったリサンプリングの簡単な実装方法としては累積重みを使った方法がありますが(MISTのチュートリアルもこのやり方)、numpyに都合のいい関数がありそうなので調べてみたらrandomモジュールのchoicを使えば楽に行けそう。

Scipyのドキュメントを見てみると

numpy.random.choice(a, size=None, replace=True, p=None)
Generates a random sample from a given 1-D array

とある。 ふーん一次元配列から適当に要素を選び取ってくれるわけね。

引数についてみてみると
- aは一次元配列かintで、配列の場合要素から適当にchoiceして値が返ってくる。intの場合np.arange(a)からchoiceしてくる
- sizeはintかintのタプルで戻り値のshapeを指定する。intが与えられた場合は指定された数だけ要素を選んで一次元配列として返す - replaceは非復元抽出か復元抽出かのフラグ(非復元は配列から同じ要素をchoiceすることを認めない)
- pは配列aの各要素が取り出される際の確率分布を表す。指定しない場合は一様にchoiceする

というわけでパーティクルの重みをpにあたえれば重みに基づいた非復元抽出ができるはず。 pythonはこういう痒いところに手が届く関数がいろいろ用意されてて助かりますね。 この方向でパーティクルフィルタを実装してみようと思いますが、それは次回。 また、非復元抽出のアルゴリズムは累積重みを使う方法のほかにもいろいろあるっぽいのでそっちも今後調べたいと思います。

プログラミングの原則

プリンシプル オブ プログラミング を読みました

先日、本屋さんをうろうろしていたらこんな本があったので、買ってみてとてもよかったので紹介します。

著者の上田勲さんはキヤノンITソリューションズで開発を行っている方で、技術書読書ブログ「Strategic Choices」を運営されているそうです。 本書では800冊を超える本を読んだ上田さんがプログラミングや開発の初心者に向けて、プログラミングのプリンシプル、先人たちの知恵をまとめて紹介してくださってます。

全体で300p程ですが、一項目につき2~3ページで読みやすいですし平易な文なので一日で読めました。

この本は初学者に向けて書かれており、僕自身、大学でプログラミングの初歩を学んでネットや書籍で勉強している最中のペーペーなので学ぶところの多い本となりました。 本書では具体的なコードを示すといったことはなですが、参考図書や出典を明らかにしてくださっているので、今後の学習にとても役に立つと思います。 この本はいわば地図のような本で、文法を学んだ初学者が次に何を学べばよいのか、そしてそれはどこにあるのかを示してくれる良書です。

一方で初学者以外は読んでも意味がないかというとそんなことはないです。 原理原則のような基本的な事柄はしばしば暗黙の了解として使用されることがあるため、初心者でなくとも「実はよく知らない」といったことがありうるように思います。 そういった言葉の原典を示してまとめているこの書籍は、初学者以外の人も、今一度原点に立ち返って自分の知識をチェックし、勉強するいいきっかけになると思います。

エンジニアリングに銀の弾丸は存在しませんが、本書で紹介されているプリンシプルをコツコツ自分のものとし、いつでも使えるものとして知識の道具箱に入れておけば、確かな力をつけられると信じて、今後も勉強しようと思います。

君の名は。を二回観ました

君の名は。を観たので感想をつらつらと
かなり売れてるようで、予定の一時間前に映画館行ったらチケット取れなかったので次の回に(それでもいい席は全部埋まってた)

これまで同様、背景美術の凄さは圧倒的、加えてキャラが可愛くてヌルヌル動くので、音なくて絵だけでも楽しめそうな作品でした。

音楽も予想以上によかったです。RADWINPSとのタイアップが発表されたときはどうなんだ?と懐疑的だったけど、新海worldにマッチする美しい音楽でした。RADWINPSというと個人的には「×と〇と罪と」の「五月の蠅」みたいなブラックな歌の印象が強いんですけど、綺麗な曲も作れて振れ幅の広いバンドなんだなぁと再認識しました。

作中では3曲もRADの歌が流れるって、監督本当にRAD好きなんだなー(instrumental含めれば20曲以上)

そういえば映画で流れる前前前世とCD音源の歌詞が微妙に違うというのを二回目で気づいた。

ストーリーに関してはシリアスな部分もあるけど会話の掛け合いとかではテンポよく笑わせて来る場面もあって、こういうコミカルさはいままでの作品にはあんまりなかったなぁと感じました。 正直、新海作品で笑ったのって初めてかもしれない。

それとコミカルさもそうだけど全体的に明るい作品ですね。 雲の向こう~や秒速~は主人公のモノローグが結構重くて、「会いたい人も、話したい人も、僕には誰もいなかった。」(雲の向こう~)とか「日々弾力を失っていく心がひたすら辛かった」(秒速)てな具合に心を刺しにくるんですよね。

基本的に以前の新海作品の主人公って孤独で、もちろん友人はいるけど自分の核となる部分に人を寄せ付けないという点は共通でした。 マンガ版の秒速だと主人公は付き合っている相手を友人に教えないって設定が出てきますし。 一方本作の主人公はあんまりそういう部分はなくて、周囲のキャラが結構深入りしてくるし、主人公もそれを口ではなんだかんだいってるけど嬉しそうだったりします。

その辺のキャラクター設定も万人受けするように作っているなぁと思いました。

映画観た後に新海監督のインタビュー記事をよんだら、「秒速のラストについて、自分(監督自身)ではBADENDとして描いたわけではないけれどお客さんからはBADENDだと受け止められていて、君の名は。ではそういう認識の違いがないようにわかりやすい作品になるようにサービスした」 というようなことを仰られていて、確かにわかりやすいし、キャッチーだし、そういう部分も今回のヒットの一因なのかなと思います。

それと、映画観ていて思ったのが、「雲の向こう~」で主人公とヒロインが再開するシーンや、「秒速」のラストで主人公とヒロインがすれ違うシーンなど 過去の作品を彷彿とさせるシーンが結構あって、「君の名は。」は過去の作品をよりわかりやすい形で再構成しているのかなといった印象を受けました。

なので「君の名は。」で新海監督を知った方は過去作品も見てみると色々楽しめると思います。あと「言の葉の庭」については小説読むとテッシーとサヤチン、ユキちゃん先生に会えます。(そういえばもし「君の名は。」がハッピーエンドじゃなかったら言の葉のラストも変わっていたかもしれないのか)

とりあえず、二回映画を観た感想はこんな感じです。もしかしたら三回目観るのでその時また気づいたことがあったら加筆します。 ごちゃごちゃ色々書きましたが、最後に一言

観てない人は観て損はないよ!観た人はもっかい観よう!ついでに過去作品も観よう!