2013-12-21

数学わかんねーからプログラムでズルしてみた　微分編（Python：Sympy）

微分ってなんぞ？

知るわきゃねーだろ！
コンビニの会計で使うくらい一般的になって出なおせや！

Wikipediaを読む

数学、とくに解析学における微分法（びぶんほう、differentiation, derivation）は、
空間やその上に定義される関数・写像を各点の近傍で考え、
その局所的な振舞いを調べることによって、それらの特徴を記述する方法である。

微分法 - Wikipedia

ふむ、局所的な振る舞いが空間で写像するわけだな。
振る舞い関数の近傍を記述する方法であると。

はい、お手上げ！

とりあえず微分する

In [4]: diff(3*x, x)
Out[4]: 3

変数が一つなので、何も指定しなくても良い。

In [3]: diff(x**2)
Out[3]: 2⋅x

定数は消える。

In [3]: diff(x**2 + 10000)
Out[3]: 2⋅x

こっからは偏微分（変数が2個以上あると、一部だけ微分するので偏微分というらしい）

In [7]: a = x**2 * sin(x + y)

In [8]: a
Out[8]: 
 2           
x ⋅sin(x + y)

In [9]: diff(a, x)
Out[9]: 
 2                            
x ⋅cos(x + y) + 2⋅x⋅sin(x + y)

In [10]: diff(a, y)
Out[10]: 
 2           
x ⋅cos(x + y)

変数が2つあるので、明示的に何で微分するか渡す必要がある。
第2引数を渡さないとエラー。

なにやら正解っぽいのが出たが、これっぽっちも微分がなんたるかを理解していない。
ググったところ、ぐにゃぐにゃしてるグラフのある一点での傾きを求めるものっぽい。

さっき定数が消えたのは、傾きと関係ないからか？

newton_methodなるものが、微分使ってがんばるらしいのでやってみた。

newton_methodとは

りんごが落ちそうな雰囲気が漂ってる。
Newton-Raphson method（ニュートン・ラフソン法）とも呼ばれるらしい。
名称はアイザック・ニュートンとジョセフ・ラフソンに由来するらしいが、ラフソンさん省略されることがあるとかかわいそう。

Wikipediaを読む

方程式系を数値計算によって解くための反復法による求根アルゴリズムの1つである。
http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95

方程式系の系ってなんだよ系って。
細かいことは気にせず先に進もう。

やってみる

Wikipediaに $sqrt(2)$ の近似値の求め方が書いてあったので、こいつを題材にしよう。

1.41421356（ヒトヨヒトヨニヒトミゴロ）を求めるのが目標。
9桁までそれっぽくでりゃいいっしょ。

In [1]: sqrt(2)
Out[1]: 
   ___
╲╱ 2 

In [2]: sqrt(2).evalf(9)
Out[2]: 1.41421356

よっしゃぁ！！！
クリアぁああああ！！

じゃなくて、newton_method使わないとね、せっかくだし。

まず、平方根をグラフで表現してみる。
方程式にするとこんな感じ。
$y=x^2$

In [1]: plot(x**2)
Out[1]: <sympy.plotting.plot.Plot at 0x108989dd0>

f:id:kaerouka:20131221172216p:plain
x=1の時、y=1
x=1.4142... y=2
x=2の時、y=4
x=3の時、y=9 ...
となる。

で、こんなグラフ書いてどうするわけ？ってとこなんだが、どうにもならない・・・
この方程式をちょっといじれば、newton_methodで解ける形になる。

$y=x^2+C$ がその方程式。

Cは切片。C=-2としてプロットしてみる。

In [1]: plot(x**2 - 2, (x, -5, 5))
Out[1]: <sympy.plotting.plot.Plot at 0x106f5af10>

f:id:kaerouka:20131221175750p:plain
この関数のyに0を代入すると
$0=x^2+2$
$x^2=-2$
$x=\pm \sqrt 2$
となり、yが0の時に、求めたい値が出るって仕組み。

こいつをsympyで見てみると

In [1]: C = symbols('C')

In [2]: solve(x**2 - C, x)
Out[2]: 
⎡   ___    ___⎤
⎣-╲╱ C , ╲╱ C ⎦

になる。
ずれてるってレベルじゃねーぞ。

ここからが本題。
Wikipediaの導入の項目を見てみる。

この方法の考え方は以下のようである：まず初めに、予想される真の解に近いと思われる値をひとつとる。
次に、そこでグラフの接線を考え、その x 切片を計算する。
http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95#.E5.B0.8E.E5.85.A5

予想される真の解・・・1.4142...ってわかってんだけど、今回は2にしてみる。

$y=x^2 - 2$ のグラフと、左記方程式のx座標が2の時の接線を同時にプロットしてみる。

まず、x座標が2の時の方程式を求める。

In [5]: Eq(y, x ** 2 - 2)
Out[5]: 
     2    
y = x  - 2

を微分すると

In [6]: diff(x ** 2 - 2)
Out[6]: 2⋅x

x座標が2の時の傾きは

In [8]: (2 * x).subs(x, 2)
Out[8]: 4

これぞ微分パワー！

点(2, 2)を通る傾き4の1次方程式のグラフは
bを切片とすると

In [10]: b = symbols('b')

In [11]: Eq(2, 4 * 2 + b)
Out[11]: 2 = b + 8

これをbについて解くと

In [12]: solve(Eq(2, 4 * 2 + b))
Out[12]: [-6]

※戻り値がlistなのは、解が複数ある方程式があるからね。 $y=x^2$ とか

これで傾きと切片が出せたので、1次方程式が書ける。
$y=4*x-6$

$y=x^2 - 2$ と一緒にプロットすんぜ！

In [13]: plot(x**2-2, 4*x-6, (x, -4, 4))
Out[13]: <sympy.plotting.plot.Plot at 0x1078a2890>

f:id:kaerouka:20131221202840p:plain
ちなみに赤文字はmac付属のソフトでちまちま加筆してます・・・

√2に近い何かが何かというと・・・
$y=4*x-6$ のyが0の時なので

In [31]: Eq(y, 4 * x -6)
Out[31]: y = 4⋅x - 6

In [32]: Eq(y, 4 * x -6).subs(y, 0)
Out[32]: 0 = 4⋅x - 6

In [33]: solve(Eq(y, 4 * x -6).subs(y, 0), x)
Out[33]: [3/2]

In [34]: solve(Eq(y, 4 * x -6).subs(y, 0), x)[0].evalf()
Out[34]: 1.50000000000000

ようやく出たよ $\sqrt 2$ に近い何か！
1.5だった！！確かに1.41421356...に近い！！

Wikipediaの導入の項目の続きをを読むと

このx切片の値は、予想される真の解により近いものとなるのが一般である。
以後、この値に対してそこでグラフの接線を考え、同じ操作を繰り返していく。
http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95#.E5.B0.8E.E5.85.A5

って書いてあるけど、最初に使ったx切片が2だね。
同じ操作を繰り返していくってのは、次はx切片を1.5にしてやってみようって事。
これを繰り返すと $\sqrt 2$ に近づきまっせってことだ。

定式化しないと使い物にならんので、定式化する。

だが、ここでまさかの気力切れ。

あるじゃない、同じ事解説してる動画。
めっちゃわかりやすい。
しかもイケメンボイス。

さて、手を抜いたところで最後にプログラム貼り付けて終わり。

def sqrt_newton_method(value, cnt):
	 # x切片
	preval = value

	for i in range(cnt):
		preval = 1/2 * (preval + value/preval)

	return preval

if "__main__" == __name__:
	# √2を1回計算
	print(sqrt_newton_method(2, 1)) #=>1.5
	# √2を2回計算
	print(sqrt_newton_method(2, 2)) #=>1.4166666666666665
	# √2を4回計算
	print(sqrt_newton_method(2, 4)) #=>1.4142135623746899

	print()
	# √3を4回計算
	print(sqrt_newton_method(3, 4)) #=>1.7320508100147274

2013-12-17

setup.py でインストールしたモジュールをアンインストールする（Python）

Python プログラミング

やり方

インストールし直す

sudo python setup.py install --record test.txt

--recordを付けると生成したファイルの一覧をtest.txtに記載してくれる。

消す

cat files.txt | sudo xargs rm -rvf

sudoを付けとかないと、うざいエラーメッセージNo.1の名をほしいままにしている
Permission deniedが顔を出すかも知れないので注意。

参考

馬の脳: python setup.py installしたモジュールのアンインストール
installation - python setup.py uninstall - Stack Overflow

2013-12-17

ファイルを読み込んで重複を削除する（Python）

Python Python3.x プログラミング

ソース

In [1]: import codecs
In [2]: set([x for x in codecs.open("yakuman.txt", "r", "utf-8")])
Out[2]: 
{'九蓮宝燈\n',
 '四暗刻\n',
 '四暗刻単騎\n',
 '国士無双\n',
 '国士無双１３面\n',
 '地和\n',
 '大三元\n',
 '大三元 四暗刻\n',
 '大三元 字一色\n',
 '大四喜\n',
 '大四喜 字一色\n',
 '天和\n',
 '字一色\n',
 '小四喜\n',
 '小四喜 字一色\n',
 '数え役満\n',
 '清老頭\n',
 '緑一色\n'}

In [3]: codecs.open("test.txt", "w").writelines(set([x for x in codecs.open("yakuman.txt", "r", "utf-8")]))

In [4]: cat test.txt
小四喜 字一色
大三元 字一色
大四喜
地和
国士無双
四暗刻単騎
大三元
九蓮宝燈
清老頭
国士無双１３面
大三元 四暗刻
四暗刻
緑一色
字一色
天和
大四喜 字一色
数え役満
小四喜

途中でさりげなく登場してるcatは、catコマンドのalias的存在。

In [6]: %alias
Total number of aliases: 16
Out[6]: 
[('cat', 'cat'),
 ('clear', 'clear'),
 ('cp', 'cp -i'),
 ('ldir', 'ls -F -G -l %l | grep /$'),
 ('less', 'less'),
 ('lf', 'ls -F -l -G %l | grep ^-'),
 ('lk', 'ls -F -l -G %l | grep ^l'),
 ('ll', 'ls -F -l -G'),
 ('ls', 'ls -F -G'),
 ('lx', 'ls -F -l -G %l | grep ^-..x'),
 ('man', 'man'),
 ('mkdir', 'mkdir'),
 ('more', 'more'),
 ('mv', 'mv -i'),
 ('rm', 'rm -i'),
 ('rmdir', 'rmdir')]

モノホンのコマンドとして実行したい時は!catとする。
最初に!をつければ、headとかpingとか、osのコマンドが使える。

データ

四暗刻単騎
国士無双
大三元
国士無双
国士無双
小四喜
小四喜
大三元
....以下略

新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）〜Part3〜

Python Python2.x やってみたプログラミング

f:id:kaerouka:20131213081901p:plain
画像の出典：paiza
https://paiza.jp/poh/ec-campaign

前回の続き

100点になった。

二分探索(bisect)使った。

ソース

# -*- coding:utf-8 -*-
import sys, bisect


#prices = ["4 3","8000","4000","9000","6000","3000","14000","10000"] #=>0, 14000, 10000
#prices = ["5 2", "4000", "3000", "1000", "2000", "5000", "10000", "3000"] #=>9000, 3000
#prices = ["4 1","1000","3000","7000","8000","10000"] #=>10000
#prices = ["10 2", "4311", "89", "2836", "9904760", "12435", "8437", "345", "4356", "43287", "9648", "8999", "43222"] #=>8782, 22083
prices = [x.rstrip() for x in sys.stdin]

N, D = [int(x) for x in prices[0].split(' ')]
m_s = [int(x) for x in prices[N+1:]]
prices = [int(x) for x in prices[1: D * -1]]
prices.sort()

for n in xrange(D):
    result = 0
    m = m_s[n]
    current_prices = [x for x in prices if x + prices[0] <= m]
    current_len = len(current_prices)
    for i, x in enumerate(current_prices[:-1]):
        # 二分探索
        j = bisect.bisect_left(current_prices, m-x, i+1)

        # j < current_len → m-xがmax(current_prices)を超えると、範囲外のindexが返されるため、その対処
        # current_pricesはソート済みなので、current_prices[i]はcurrent_prices[i-1]より必ず大きいため
        # 前回までの値（result値）との比較は必要ない。
        if j < current_len and current_prices[i] + current_prices[j] == m:
            result = m
            break
        # 一つ上のif文の条件を満たさず、j == i+1を満たす場合、current_prices[i] + current_prices[j] <= mを満たさなかったと判断
        elif j == i+1:
            continue
        else:
            # 1つ前。
            j -= 1
            total = current_prices[i] + current_prices[j]
            if result < total <= m:
                result = total
                if result == m:
                    break

    print result

相変わらず入力値の受け取りが汚いが、目を瞑る。
めんどいので動作確認用の入力値もそのまま。

結果

f:id:kaerouka:20131216015643p:plain
画像の出典：paiza
https://paiza.jp/poh/ec-campaign

オールクリア！

f:id:kaerouka:20131216015754p:plain
画像の出典：paiza
https://paiza.jp/poh/ec-campaign

kaeroukaさん、凄いコードですねっ！！私にはこんな効率のいいコード書けませんっ！

かわいい。
彼女にはちゃんと仕様を読むプログラマーになっていただきたいところです（しみじみ

実行速度だが、Case3で５秒弱かかっている。
最速の人が0.09秒とか鬼畜の所業。

ここまでのチューニングはわりと普通なやり方だと思うんだけど、これ以降はどうすんのかねぇ。
いい結果が出たらまた記事にしよう。

過去の記事
 新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python） - Tips
新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）〜Part2〜 - Tips

2013-12-15

新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）〜Part2〜

Python Python2.x やってみたプログラミング

f:id:kaerouka:20131213081901p:plain
画像の出典：paiza
https://paiza.jp/poh/ec-campaign

前回までの流れ

適当にやってみたところ、遅すぎる上に見事にバグりまくってたという醜態を晒した。
真面目にソース読んでみた人ごめんなさい。
しかもソース汚すぎ。
改善した

今回は少しだけ改善してみた。
[100, 200, 250, 120]
こんなデータを真面目にループで回すと
In [42]: for i in [100, 200, 250, 120]:
   ....:     for j in [100, 200, 250, 120]:
   ....:         
こんな感じになるわけだけど
最初のiとjの組み合わせがi=100, j=100
さらに、i=100, j=200、i=100, j=250、i=100, j=120、i=200, j=100
と続くわけだけど、i=200, j=100の組み合わせって、i=100, j=200と同じだからいらないよねって考え。
ようは、j=i+1から始めようって事。
同じ商品を選んじゃいけないので、iは120の行わない。
# python 2.7.3
import sys

prices = [x.rstrip() for x in sys.stdin]

N, D = [int(x) for x in prices[0].split(' ')]
m_s = [int(x) for x in prices[N+1:]]
prices = [int(x) for x in prices[1: D * -1]]

for n in xrange(D):
    result = 0
    for i, x in enumerate(prices[:-1]):
        for y in prices[i+1:]:
            total = x + y
            if total <= m_s[n] and result < total:
                result = total
    print result
結果

画像の出典：paiza
お、Case1クリア
王大人「Case2で死亡確認」

画像の出典：paiza
https://paiza.jp/poh/ec-campaign

kaeroukaさん、、すごーい。

かわいい。
Case2、Case3とクリアしていったら、服も1枚ずつ減ってったりするんだろうか。

次回は2分探索でチューニング予定。
Case3もクリアできたら嬉しいけど、どうだろう。
前回
 新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python） - Tips
続き
 新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）〜Part3〜 - Tips

2013-12-13

新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）

Python2.x Python やってみたプログラミング

新人女子プログラマの...とは

f:id:kaerouka:20131213081901p:plain
画像の出典：paiza
https://paiza.jp/poh/ec-campaign

新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！
ECサイト内の2つの異なる商品（値段は同じでも構わない）を購入し、その合計価格が指定の価格以内で最大になる組み合せを探してください。

だそうです。
いわゆるハッカソンってやつで、みんなでコード書いて実行速度やら可読性やら比べようぜって企画です。
新人女子プログラマという言葉で引っ掛けようって魂胆が見え見えですね。
はい、釣られました。

お仕事の詳細

は、やたら長いので↑のリンク先見てね。

とりあえず動くもの

新人女子プログラマ「納期に間に合わない…」

とか言ってるわけですよ。
チューニングは後回しにして、とりあえず動くものを作らないといけないよね。
遅くても動いてればクライアントへの進捗報告もやりやすい。
おじさんが面倒見てやろうじゃないの。

ソース（バグってるので読み飛ばし推奨）

# python 2.7.3
import sys, math

prices = [x.rstrip() for x in sys.stdin]

N, D = [int(x) for x in prices[0].split(' ')]
m_s = [int(x) for x in prices[N+1:]]
prices = [int(x) for x in prices[1: D * -1]]
results = [0] * D

for n in xrange(D):
    for x in prices:
        for y in prices:
            a = math.fabs(m_s[n] - (x + y))
            b = math.fabs(math.fabs(results[n]) - math.fabs(m_s[n]))
            if a < b:
                results[n] = x + y

for x in results:
    print(x)

※2013/12/13追記仕様よく読んでなくて、結構バグってる事が発覚

単純に1個ずつ足してって、今まで1番近かった解答と比較し
今回の比較の方が正解に近ければ、解答を書き換えるって仕組み。
計算量は $O(D*N^2)$ かな。
どんだけmath.fabs呼んでんだよ。

結果は・・・・

f:id:kaerouka:20131213072743p:plain
画像の出典：paiza

うおぉおおおい！タイムオーバー！？
15秒以内なんて要件あったんかい。

f:id:kaerouka:20131213074000p:plain
画像の出典：paiza

kaeroukaさん、ある意味凄いコードですね。
新人の私でもそうそうこんなコードかけませんよ。。

この女殺すわ。

続きは後日。

2013/12/13追記
CASE2のサンプルで動かしてみたら、正解と合わない。
よくよく問題見てみたら、2つの異なる商品と書いてあった罠。
さらに設定金額以下でって書いてあるから、math.fabsとかなんの関係も・・・
さて、そろそろ真面目にやらないと・・・
明日から本気出す！

続き↓
新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）〜Part2〜 - Tips
新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）〜Part3〜 - Tips

2013-12-13

複数端末にpingする（Python3）

Python3.x Python プログラミング

隣の人のPC上で動いてるVMに乗り込んで作業するはずだったんだけど
IP変わったらしく、繋がらない。
だいたい範囲わかってるから乱れ打ち。

ちょっと無理する

ソース

In [1]: import os, subprocess

In [2]: [x + ":" + str(subprocess.call(["ping", "-c", "1", "-t", "1", x], stdout=open(os.devnull,"w"))) for x in ["192.168.111." + str(y) for y in range(50,55)]]
Out[6]: 
['192.168.111.50:0',
 '192.168.111.51:2',
 '192.168.111.52:2',
 '192.168.111.53:2',
 '192.168.111.54:2']

わけわからん。

普通に書いた

ソース

import os
import subprocess

# 192.168.111.50 〜 192.168.111.54のリストを作成
L = ["192.168.111." + str(x) for x in range(50, 55)]
result = []
for ip in L:
    # pingを打って結果（ステータスコード）をipと共にリストに格納
    # 標準出力を/dev/nullへ、pingの回数は1回、タイムアウト1秒
    #引数の都合上、Linux系じゃないと動かない
    r = subprocess.call(["ping", "-c", "1", "-t", "1", ip], stdout=open(os.devnull, "w"))
    result.append("%s:%d" % (ip, r))

#表示
[print(x) for x in result]

結果

192.168.111.50:0
192.168.111.51:2
192.168.111.52:2
192.168.111.53:2
192.168.111.54:2

遅い。最悪ping数*1秒かかる。
fpingコマンド使えばわざわざこんなことしｎ・・・

Tips

IT技術系Tips

数学わかんねーからプログラムでズルしてみた　微分編（Python：Sympy）

微分ってなんぞ？

Wikipediaを読む

とりあえず微分する

newton_methodとは

Wikipediaを読む

やってみる

setup.py でインストールしたモジュールをアンインストールする（Python）

やり方

参考

ファイルを読み込んで重複を削除する（Python）

ソース

データ

関連

新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）〜Part3〜

100点になった。

ソース

結果

新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）〜Part2〜

前回までの流れ

改善した

結果

新人女子プログラマの書いたコードを直すだけの簡単なお仕事です！ ←やってみた（Python）

新人女子プログラマの...とは

とりあえず動くもの

ソース（バグってるので読み飛ばし推奨）

結果は・・・・

複数端末にpingする（Python3）

ちょっと無理する

ソース

普通に書いた

ソース

結果