2011-11-13
Pythonで線形SVM。
こちらの記事を参考にしました。
参考というか、もはや写経です。
非常に勉強になりました。
#coding:shift-jis
import numpy as np
from scipy.linalg import norm
import cvxopt
import cvxopt.solvers
from pylab import *
# 全データの総数。
N = 10
# クラス1のデータ。
cls1 = [(1,6),
(2,10),
(3,8),
(4,11),
(5,9)]
# クラス2のデータ。
cls2 = [(1,3),
(2,2),
(3,4),
(4,1),
(5,3)]
# クラス1には1、クラス2には-1を割り振ってarrayにする。
t = []
for i in range(N/2):
t += [1.0]
for i in range(N/2):
t += [-1.0]
t = array(t)
"""データ行列Xをvsstackで作成。以下のようにスタック。X[1,0]で2とか。X[行,列]。
[[ 1 6]
[ 2 10]
・・・
[ 4 1]
[ 5 3]]
"""
X = vstack((cls1, cls2))
# 線形カーネル。x*yの総和。
def kernel(x, y):
return np.dot(x, y)
#10*10の0行列を用意する。
K = np.zeros((N, N))
# Kを作る。
for i in range(N):
for j in range(N):
K[i,j] += t[i]*t[j]*kernel(X[i],X[j])
#Kはこのmatrixにしないとらめ。
Q = cvxopt.matrix(K)
# マイナスにするのでこのpが必要。-1がN個の列ベクトル。
p = cvxopt.matrix(-np.ones(N))
# 不等式制約。
G = cvxopt.matrix(np.diag([-1.0]*N)) # 対角成分が-1のNxN行列
h = cvxopt.matrix(np.zeros(N)) # 0がN個の列ベクトル
# 等式制約。w*xの総和が0。
A = cvxopt.matrix(t, (1,N)) # N個の教師信号が要素の行ベクトル（1xN）
b = cvxopt.matrix(0.0) # 定数0.0
sol = cvxopt.solvers.qp(Q, p, G, h, A, b) # 二次計画法でラグランジュ乗数aを求める
a = array(sol['x']).reshape(N) # 'x'がaに対応する。ラグランジュ乗数aのリスト。
# サポートベクトルのインデックスのリストを。
S = []
for i in range(len(a)):
if a[i] < 0.00001: # ほぼ0ならサポートベクトルじゃない。
continue
else:
S += [i]
# wを計算。
w = np.zeros(2) #行列を用意。
for n in S:
w += a[n]*t[n]*X[n]
# bを計算。
sum = 0
for n in S: # サポートベクトルの数だけforループ。
temp = 0
for m in S: # サポートベクトルの数だけforループ。
temp += a[m]*t[m]*kernel(X[n],X[m])
sum += (t[n]-temp)
b = sum/len(S)
# データをプロット。
for p in (cls1+cls2):
plot(p[0],p[1],'bo')
# サポートベクトルを描画
for n in S:
plot(X[n,0], X[n,1], 'ro')
# 識別関数を描画。
def f(x1, w, b):
return - (w[0] / w[1]) * x1 - (b / w[1])
# 識別境界を描画
x1 = np.linspace(-6, 6, 1000)
x2 = [f(x, w, b) for x in x1]
plot(x1, x2, 'g-')
xlim(0, 6)
ylim(0, 12)
show()
SVMについてはある程度理解できました。
頭のとんでもなく悪い自分としては結構嬉しい。
数式弄りながら実装するの楽しかったなあ。
まあ、最初はなかなか理解できず死にたくなりましたが笑。
それよりも、Pythonでの行列計算に慣れないとらめだと痛感。
cvxoptにつてもいろいろ知っておかないとだめぽ。
線形代数は院試で割と勉強したけれどちょいちょい忘れかけてるし、復習しないと。
次は非線形SVMだー！
これはさらにおもしろそうだー！！！
Permalink | コメント(0) | トラックバック(0) | 20:42 
2011-11-12
交差確認法をPythonで。
「集合知プログラミング」の続きです。8章なう。
とりあえずでけたのでうｐします。
本のサンプルコードではなく、また自分で勝手に実装しました。
交差確認法を使って、以前実装したk-NN法の誤差評価をしています（たぶん）。
身長と体重にはもちろん相関がありますが、
持っているエロゲの数と体重にはなんら相関はありません。
（変数名てきとーすぐる…）
#coding:shift-jis import numpy as np from pylab import * import random import math # 元データは[(身長,持ってるエロゲの数,体重)]。 data = [] for x in range(130,201): r = random.random() n = random.randint(1,30) data += [(x, n, (1+r/5)*22.0*((x/100.0)**2))] # 元データの各軸。 s =[] e = [] t = [] for (x,y,z) in data: s += [x] #身長リスト e += [y] #エロゲ数リスト t += [z] #体重リスト # 各軸の平均。 def heikin(list): return sum(list)/len(list) # 標準偏差の逆数を出す関数。 def hensa(list): a = 0 m = heikin(list) for x in list: a += (m-x)**2 return 1/sqrt(a/len(list)) # リストを正規化。 def seikika(list): a = [] s = hensa(list) for i in list: a += [i*s] return a # 正規化した各々の軸データ。 shinchyo = seikika(s) eroge = seikika(e) taijyuu = seikika(t) # 正規化バージョンの[(身長,持ってるエロゲの数)]。 data1 = [] for d in zip(shinchyo,eroge): data1 += [d] # 正規化バージョンの[(身長,体重)]。 data2 = [] for d in zip(shinchyo,taijyuu): data2 += [d] # ガウス関数（距離が近い時1〜距離が遠い時0近く）までを重みとして取れる。 def gaussian(x,sigma=10.0): return math.e**(-x**2/(2*sigma**2)) # ガウス関数による重み付けK近傍法。 # 入力x（身長）との距離を全データで出して近いもの上位3つを。 def k_nn(data,x,k): dlist = [] #与えたxと距離が近いデータのリスト。 sum = 0 #体重*距離のガウス関数重みの総和。 w = 0 #重みの総和 for d in data: dlist += [( (x-d[0])**2 , d[1] )] #(距離,体重)のタプルをdlistに入れていく。 dlist.sort() #距離でソートする。 for d in dlist[0:k]: #距離が短い上位k個をひとつずつ。 sum += d[1]*gaussian(d[0]) #体重に距離のガウス関数重みを。 w += gaussian(d[0]) #距離のガウス関数重みの総和。 return (x,sum/w) #(身長,k_nnで割り出した体重） # このk近傍法でどれだけ良く未知データを予測できるかを交差確認法（一つ抜き）で評価する。 def cross(data): sum = 0 for train in data: data.remove(train) a = train[1] #正解の体重。 b = k_nn(data,train[0],3)[1] #k_nnで予測する体重。 sum += abs(a-b)**2 #誤差の二乗の総和。 data.append(train) return sum # k_nn法で身長からエロゲの数を算出した際の、実際のエロゲの数との誤差。 print(cross(data1)) # k_nn法で身長から体重を算出した際の、実際の体重との誤差。 print(cross(data2))
結果はこう。
95.1400375937 3.56766330246
交差確認法での誤差が大きければ大きいほど評価方法としては悪いものになります。
この場合、その人が持っているエロゲの数から体重を予測するのはアホということですはい。
今回、持っているエロゲの数と体重に関係がなさそうなことくらい誰でも解りますが、
扱うデータの変数が多く複雑になってくると、
どの変数が重要でどの変数が不要か、事前に解るとは限りません。
交差確認法で得た値（誤差）が大きすぎる場合には、この変数を消去する等していくわけですね。
これを一種のコスト関数として今後いろいろ使っていくことになりそうです。
（うーん…なんかまちがってそうなのであまり信用しないでくださいorz）
Permalink | コメント(0) | トラックバック(0) | 22:55
アニーリングアルゴリズム（改）。
前回のアニーリングアルゴリズムの記事はサンプルデータがアレでした。
というわけでちょっと書き換えました。
#coding:shift-jis import numpy as np from pylab import * import random import math data = [0] for i in range(1000): a = random.randint(1,100) if (a==1): b = random.randint(100,200) for i in range(b): data += [data[-1]-3] else: data += [data[-1]+3] # データを青でプロット。 for (x,y)in zip(range(len(data)),data): plot(x,y,'bo') # アニーリングアルゴリズム。 def annealing(data,T=10000,p=0.95): r = random.randint(0,len(data)-1) #初期値をランダムに選択。 plot(r,data[r],'yo') #初期値を黄色でプロット。 while True: cost_now = data[r] #現時点でのコスト。 cost_new = data[r+1] #次点でのコスト。 if (T >= 0.1): #高コストの解も受け入れる時。 p = e**(-abs(cost_now-cost_new)/T) #高コストの解を受け入れる確率。 if (cost_now > cost_new) or (p > random.random()): r = r+1 #次点へ移動。。 else : r = r-1 #逆のr-1の点へ移動。 T *= 0.95 #Tを小さくする。 elif (T < 0.1): #Tが十分小さくなったら（高コストの解を受け入れなくなったら）、ヒルクライムで局所最小解を出す。 if((data[r-1]) > cost_now < cost_new): break #最終的にループを抜け出す。 elif (cost_now > cost_new): #次点へ。 r = r+1 else : r = r-1 return ([r,data[r]]) #アニーリングアルゴリズムによる最小解を赤でプロット。 a = annealing(data) plot(a[0],(a[1]),'ro') show()
結果は以下の通りです。
データは青線、初期値は黄点、結果は赤点で示しています。
普通のヒルクライムなら初期値の黄点から下降して初めて行き着いた底が赤点になるはずですが、
アニーリングなのでそうなってないですね。一応成功です。
ふぇぇ〜
Permalink | コメント(0) | トラックバック(1) | 13:40
アニーリングアルゴリズム。
（書き直した記事はこちら。）
「集合知プログラミング」の続きです。5章の最適化あたりを読んでます。
ヒルクライム（山登り法）というものは知っていたのですが、
アニーリング（焼きなまし法）というものは知りませんでした。
坂道を下るだけじゃなく時には登ってやろうぜというのが、アニーリングです（てきとう。
本のサンプルコードが自分の頭の悪さが原因でなんだか理解しにくかったので、
簡単な具体例を使って自分で書いてみました。
あと、アニーリングのみだと坂道の途中で終わるかもしれないので、
最終的にはヒルクライムで近傍の極小に確実に向かわせるようにしました。
#coding:shift-jis import numpy as np from pylab import * import random import math # サンプルデータ[[旅費,乗車時間,乗り継ぎ待ち時間]]の作成。500プロット。 data = [] for i in range(500): a = random.randint(10000,30001) #旅費（10000円から30000円）。 b = random.randint(60,301) #乗車時間（60分から300分）。 c = random.randint(5,61) #乗り継ぎ待ち時間（5分から60分）。 data += [[a,b,c]] # コスト関数。[旅費,乗車時間,乗り継ぎ待ち時間]からコスト計算。 # 乗車時間は+60分で1000円、乗り継ぎ待ち時間は+60分で+1500円のコストとする。 def costf(data): return (data[0]+(1000*data[1]/60)+(1500*data[2]/60)) # アニーリングアルゴリズム。 def annealing(data,T=10000,p=0.95): r = random.randint(0,len(data)-1) #初期値をランダムに選択。 while True: cost_now = costf(data[r]) #現時点でのコスト。 cost_new = costf(data[r+1]) #次点でのコスト。 if (T >= 0.1): #高コストの解も受け入れる時。 p = e**(-abs(cost_now-cost_new)/T) #高コストの解を受け入れる確率。 if (cost_now > cost_new) or (p > random.random()): r = r+1 #次点へ移動。。 else : r = r-1 #逆のr-1の点へ移動。 T *= 0.95 #Tを小さくする。 elif (T < 0.1): #Tが十分小さくなったら（高コストの解を受け入れなくなったら）、ヒルクライムで局所最小解を出す。 if(costf(data[r-1]) > cost_now < cost_new): break #最終的にループを抜け出す。 elif (cost_now > cost_new): #次点へ。 r = r+1 else : r = r-1 return ([r,data[r]]) # 確認のため全コストを青でプロット。 for (n,p) in zip(range(500),data): plot(n,costf(p),'bo') # 軸 xlim(0, 500) ylim(10000, 40000) #アニーリングアルゴリズムによる最小解を赤でプロット。 a = annealing(data) plot(a[0],costf(a[1]),'ro') # 散布図を描く。 show()
ここまで書いていまさら思ったけど、サンプルデータがランダムなのでだめだorz
アニーリングを試すのに良さげなサンプルデータ作らないと。
なんかないかなあ。
アルゴリズムはだいじょぶなはず。
- 作者: Toby Segaran,當山仁健,鴨澤眞夫
- 出版社/メーカー: オライリージャパン
- 発売日: 2008/07/25
- メディア: 大型本
- 購入: 72人 クリック: 1,836回
- この商品を含むブログ (262件) を見る
Permalink | コメント(0) | トラックバック(0) | 09:27
2011-11-11
Pythonでパーセプトロン。
最初Haskellで書こうとして死にました。
こういうプログラムはあまりHaskellには向かないようです（重みを変えてループさせたりとか）。
というわけでPythonでやることに。
こちらの記事が参考になりました。
本当に助かりました。ありがとうございます（ぺこり。
まだイマイチ確率的最急降下法というものが解っていないのですが（ぉぃ、
なにはともあれ書けたのでうｐ。
すごーくカンタンな例でやってみました。
（しかしなぜブログにうｐするとタブがずれるんだ…。）
（あと、np.dotうんちゃらとかarrayうんちゃらとかを使うと、
HaskellのzipWith的な計算できるんだほえーって思った。知らなかった。）
#coding:shift-jis import numpy as np from pylab import * # サンプルデータ。[(x,y,クラス)]。クラスは1か-1の値。 data = [(1,3,1), (2,4,1), (3,2,1), (1,5,-1), (2,7,-1), (3,6,-1)] # データを散布図として描画。 for p in data: if p[2]==1: plot(p[0], p[1],'ro') #クラス1ならプロットを赤に。 else: plot(p[0], p[1],'bo') #クラス-1ならプロットを青に。 # 重みの初期値。 w = array([1.0, 1.0, 1.0]) # 学習係数。 p = 0.01 # 分類が正解したデータ数 correct = 0 # パーセプトロンアルゴリズム。重み調整。 while True: for d in data: # 全データについて検討。 if np.dot(w, [1, d[0], d[1]]) * d[2] > 0: #分類が正しいならcorrect+1。 correct += 1 else: # 誤分類なら重みを調整。 w += p * array([1, d[0], d[1]]) * d[2] if (correct==len(data)): #全てが正しく分類されたらbreak。 break else: #だめならcorrectを0に戻してループ。 correct = 0 continue # プロットする。 x = linspace(0,5) y = -w[1]/w[2] * x - w[0]/w[2] plot(x, y) #グラフの縦軸横軸の範囲。 xlim(0, 5) ylim(0, 10) # 最終的にグラフ描く！ show()
結果はこうなります。
データが6つしかないのでアレですが、とりあえずちゃんと線形分離されています。
Pythonでグラフが描けるなんて初めて知りました（今更。
またいろいろ楽しくなりそうです。
今後は「集合知プログラミング」や「フリーソフトでつくる音声認識システム」の続きを読んだりしながら、
個人的に一番興味のあるサポートベクターマシン（SVM）を実装していきたいと思います。
なんとか日曜までに形にしたいなあ。
これもまた、こちらの方の記事がものすごく参考になります。
その後非線形のSVMもやりたいです。これはさらに楽しそう。
本当に、本当にありがとうございます（ぺこり。
Permalink | コメント(0) | トラックバック(0) | 12:14
2011-11-08
k近傍法（k-NN法）を今度はPythonで。
一つ前のエントリーの例があまりに酷かったので、
今回はとんでもなく単純な例だけれどまともな例で、k近傍法をPythonで書きました。
距離による重みはガウス関数でやりました。
ワインの価格を年代から決めます。
集合知プログラミングの本の例をさらにカンタンにしました。
#coding:shiftjis import math # ガウス関数（距離が近い時1〜距離が遠い時0近く）までを重みとして取れる。 def gaussian(x,sigma=10000): return math.e**(-x**2/(2*sigma**2)) # ワインの[(年代,価格)]。5年おきのデータのリスト。 wine_data = [(1950,5000), (1955,8000), (1960,10000), (1965,20000), (1970,45000), (1975,50000), (1980,42000), (1985,38000), (1990,21000), (1995,10000), (2000,5000) ] # ガウス関数による重み付けK近傍法なう。 def k_nn(data,year,k): dlist = [] #指定したyearに距離が近い順のワインリスト。 sum = 0 #価格*距離のガウス関数重みの総和。 w = 0 #重みの総和 for d in data: dlist += [( (year-d[0])**2 , d[1] )] #(距離,価格)のタプルをdlistに入れていく。 dlist.sort() #距離でソートする。 for d in dlist[0:k]: #距離が短い上位k個をひとつずつ。 sum += d[1]*gaussian(d[0]) #価格に距離のガウス関数重みを。 w += gaussian(d[0]) #距離のガウス関数重みの総和。 return sum/w # 1968年のワイン価格を、近い年代3つのワイン価格の平均として予測。 print(k_nn(wine_data,1968,3))
なんか小学生レベルのコードですがとりあえずうｐ。
- 作者: Toby Segaran,當山仁健,鴨澤眞夫
- 出版社/メーカー: オライリージャパン
- 発売日: 2008/07/25
- メディア: 大型本
- 購入: 72人 クリック: 1,836回
- この商品を含むブログ (262件) を見る
Permalink | コメント(0) | トラックバック(0) | 20:40
k近傍法（k-NN法）をHaskellで書いてみました。
（データ構造的なところでなぜか死にかけましたがなんとかうｐ…orz）
k近傍法と言ってもいろいろあるのですが、
一番単純な多数決だとアレ過ぎるので、多数決と距離の塩梅を考えて計算してみました。
ちなみに、その計算法で大丈夫なのかよく解りません（ぉぃ。
ただk近傍法を実装すると言ってもつまらないので、
今回はロリコン判定関数というものを例にして書いてみました。
your_zokusei関数の二つの引数に、
あなたのお好みの身長とお好みのおっぱいの大きさを入れればロリコン判定ができます。
※ロリが身長とおっぱいの大きさで決まるかどうかは僕の専門外なので解りません。
ロリコンではなく貧乳好きということもあるかもしれません。
つまり、あまり考えてません。とりあえず2つのパラメータがほしかっただけです。
（コードのタブずれるのなんで…まいっか…そして変数名てきとうすぐるｗ）
import Data.List
-- 正規化の準備
-- 平均を求める。
heikin xs = realToFrac(sum xs)/realToFrac(length xs)
-- 分散を求める
bunsan xs = (sum $ map (^2) $ map (+(heikin xs)) $ map (*(-1)) xs)/realToFrac(length xs)
-- 標準偏差を求める。
hensa xs = sqrt $ bunsan xs
-- 軸を正規化する関数。
seikika xs = map (/(hensa xs)) xs
-- データを正規化
-- (名前,身長,おっぱい,属性)のデータ構造で扱う。
girls_Data = [("Yui",143,73,"Loli"),
("Mio",163,90,"not Loli"),
("Azunyan",139,75,"Loli"),
("Kyou",160,82,"not Loli"),
("Tomoyo",164,86,"not Loli"),
("Fuko",135,73,"Loli"),
("Ayu",143,77,"Loli"),
("Nayuki",165,84,"not Loli"),
("Akiko",168,92,"not Loli"),
("Koppora",138,63,"Loli")]
-- 名前を取り出す。
aa = map (\(a,b,c,d) -> a) girls_Data
-- 身長を正規化。
bb = seikika $ map (\(a,b,c,d) -> b) girls_Data
--おっぱいを正規化。
cc = seikika $ map (\(a,b,c,d) -> c) girls_Data
-- 属性を取り出す。
dd = map (\(a,b,c,d) -> d) girls_Data
-- 以上をまとめて正規化したデータを作成。
girls_Data2 = zip4 aa bb cc dd
-- 身長の標準偏差の逆数。
xx = 1/(hensa $ map (\(a,b,c,d) -> b) girls_Data)
-- おっぱいの標準偏差の逆数。
yy = 1/(hensa $ map (\(a,b,c,d) -> c) girls_Data)
-- ここからメインのプログラム
-- 入力データと元データとの距離を計算。(距離,属性)。
f :: (Num t1) => (t1, t1) -> (t, t1, t1, t2) -> (t1, t2)
f (x,y) (name,a,b,c) = ((a-x)^2+(b-y)^2,c)
-- 距離が近い学習データ上位n個の属性のリスト。[(距離,属性)]。
g :: (Double, Double) -> Int -> [(Double, [Char])]
g (x,y) n = take n $ sort $ map (f (x*xx,y*yy)) girls_Data2
-- [(属性,距離)]から、(距離の平均*(1/要素数) , 属性)を出す。
-- 距離は小さいほど良いのでそのまま平均を反映させ、
-- 多数決のため要素数は多いほど良いので逆数をかけて調整する
h :: (Fractional b) => [(a, b)] -> (b, a)
h xs = ( (sum $ map snd xs)/(realToFrac(length xs)^2) , (fst $ head xs) )
-- fstが同じものでグループ化する。ここでは[(属性,値)]について、属性でグループ化したい。
ff :: (Eq a) => a -> (a, t) -> Bool
ff a (b,c) = a==b
groupFst :: (Eq a) => [(a, b)] -> [[(a, b)]]
groupFst [] = []
groupFst xs = (takeWhile (ff a) xs) : groupFst (dropWhile (ff a) xs)
where a = fst $ head xs
-- 例えば身長155、おっぱい80の女の子の場合
-- 単純な多数決ではnot Roriになるが、最終的な判定では距離も考慮するのでRoriになったりするよん。
-- 距離が近い上位3つを取得。
a = g (155,80) 3
-- 属性でソート。
b = sort $ map (\(a,b) -> (b,a)) $ a
-- 属性でグループ化。
c = groupFst b
-- 重み付け多数決距離的なものを算出。
d = sort $ map h c
-- 最も距離が近い属性を出す。
e = snd $ head d
-- a〜eをまとめて関数化すると以下の通り
-- xにお好みの身長、yにお好みのおっぱいの大きさを入力して、ロリコンかどうか判別する関数。
your_zokusei :: Double -> Double -> String
your_zokusei x y = "You are " ++ (snd $ head $ sort $ map h $ groupFst $ sort $ map (\(a,b) -> (b,a)) $ g (x,y) 3) ++ "-con!"
例えば身長が140、おっぱいが70の女の子が好きな場合は以下のようにします。
*Main> your_zokusei 140 70 "You are Loli-con!"
ちゃんと判定されました。
元データが10人と少ないので判定はもちろんてきとーです。
また、WORKING!!のぽぷらちゃんのような身長が極端に小さくおっぱいが極端に大きい、
つまり身長的にはロリだけどおっぱい的にはロリじゃないというような特異な女の子が元データに入ってくると、
判定がちょっとおかしくなってきてしまいます。
パターン認識というのはつまりそういうものです（違。
まあ、ここまで書いておいて言うのもなんですが、これは酷い…。
というわけで、もう筆を置きます。
フリーソフトでつくる音声認識システム - パターン認識・機械学習の初歩から対話システムまで
- 作者: 荒木雅弘
- 出版社/メーカー: 森北出版
- 発売日: 2007/10/17
- メディア: 単行本（ソフトカバー）
- 購入: 30人 クリック: 362回
- この商品を含むブログ (35件) を見る
- 作者: Toby Segaran,當山仁健,鴨澤眞夫
- 出版社/メーカー: オライリージャパン
- 発売日: 2008/07/25
- メディア: 大型本
- 購入: 72人 クリック: 1,836回
- この商品を含むブログ (262件) を見る
Permalink | コメント(0) | トラックバック(1) | 20:18
2011-11-06
集合知プログラミング
集合知プログラミングという本を買いました。
つまりは機械学習なんですが、割とおもしろいです。
今後これで遊んでいこうと思います。
- 作者: Toby Segaran,當山仁健,鴨澤眞夫
- 出版社/メーカー: オライリージャパン
- 発売日: 2008/07/25
- メディア: 大型本
- 購入: 72人 クリック: 1,836回
- この商品を含むブログ (262件) を見る
昨日今日で第2章が終わりました。
実際自分でコードを書いたので、一応うｐしておきます。
サンプルコードに頼る部分も結構あったのですが、
自分なりにちょいちょい書き換えました。
※以下のサンプルコードでのしょくでんさん、たうえさんは架空の人物です。
#coding:shift-jis
from math import sqrt
critics = {
'椿':{'G線':9.5,'車輪の国':9.8,'CLANNAD':10.0,'失われた未来を求めて':9.8,'Kanon':9.9,},
'しょくでんさん':{'車輪の国':7.0,'カタハネ':8.45,'CLANNAD':7.0,'G線':6.2,'おねむりしすたー':9.5,'熟女狩り':2.000,'失われた未来を求めて':5.5,'ろりこん！':10.0},
'たうえさん':{'CLANNAD':10.0,'カタハネ':8.9,'車輪の国':9.9,'おねむりしすたー':3.0,'熟女狩り':9.900,'Kanon':9.3,'ろりこん！':1.2}
}
# ーーーーーーーーユーザーベースのフィルタリングーーーーーーーー
def sim_person(dic,p1,p2):
# p1、p2の二人の類似度を出す。
gameList = []
for item in dic[p1]: # p1の美少女ゲームを1つずつitemに入れる。
if item in dic[p2]: # p1が持っているゲームをp2も持っているか。
gameList += [item] #共通するゲームのリストを作成。
# 二人が共通して持っているゲームの数。
n = len(gameList)
# 共通するゲームがなければ0。
if (n ==0):
return 0
# 共通するゲームの評価点のリスト。
pointList1 = [dic[p1][it] for it in gameList]
pointList2 = [dic[p2][it] for it in gameList]
# 評価点の合計。
sum1 = sum(pointList1)
sum2 = sum(pointList2)
# 評価点の平方を合計。
def sq(x):
return (x**2)
sqSum1 = sum(map(sq,pointList1))
sqSum2 = sum(map(sq,pointList2))
# 二人の評価点それぞれのゲームの積を合計する。
pSum = sum([dic[p1][it]*dic[p2][it] for it in gameList])
# ピアソン相関係数を出す。
num = pSum - ((sum1*sum2)/n)
den = sqrt((sqSum1-(sum1**2)/n)*(sqSum2-(sum2**2)/n))
# 分母0なら0
if den == 0:
return 0
return (num/den)
def recome(dic,p1):
# おすすめのゲームとその評価点を出す。
totals = {} #{p1がやっていないゲーム:(他人の評価点*自分との類似度)の合計}の辞書。
simSum = {} #{p1がやっていないゲーム:類似度の合計}の辞書。
for other in dic: # 一人ずつotherに入れていく。
if (other != p1):
s = sim_person(dic,p1,other) #otherがp1じゃないなら類似度を出す。
for item in dic[other]: #次はotherのゲームについて。
if item not in dic[p1]: #p1がやっていないゲームなら。
totals.setdefault(item,0)
totals[item] += dic[other][item]*s #その人の評価点に類似度をかけてtotalに入れる。
# p1のやっていないゲームの類似度を合計。
simSum.setdefault(item,0)
simSum[item] += s
# [(最終的な評価点,ゲーム)]というリストを作る。類似度の合計が0の場合は除く（0除算）。
rankings = [(total/simSum[item],item) for item,total in totals.items() if simSum[item] != 0]
#ソートしてリバースして最終的なランキングへ。
rankings.sort()
rankings.reverse()
return rankings
# 上位のn位までのオススメを出す。
def best(n,ranking):
count = 1
for r in ranking[0:n]:
print (str(count) + '位 '+ r[1] + ' ' + str(r[0]))
count += 1
print(best(3,recome(critics,'椿')))
# ーーーーーーーーーーーユーザーベースここまでーーーーーーーーーーー
# ーーーーーーーここからアイテムベースのフィルタリングーーーーーーーー
# criticsの{ユーザー:{ゲーム:評価点}}から、p1（人）を指定して[(p1との類似度,他のユーザー)]にできる。
# ということはつまり、データ構造の対称性から考えて以下のようにもできる。
# {ゲーム:{ユーザー:評価点}}から、g1（ゲーム）を指定して[(g1との類似度,ゲーム)]と。
# 次からはユーザー同士の類似度ではなく、ゲーム同士の類似度を考える。
def match3(dic,p1):
# p1と類似度の高い上位3人をぴっくうｐする。
#[(p1との類似度,他のユーザー)]というリスト
scores = [(sim_person(dic,p1,other),other) for other in dic if other != p1]
scores.sort()
scores.reverse()
return scores[0:4]
# [ゲーム:{ユーザー:評価点}]というリストにする。
def trans(dic):
result = {}
for person in dic: #人をひとりひとり。
for item in dic[person]: #その人のゲームをひとつひとつ。
result.setdefault(item,{})
# item（ゲーム）とperson（ユーザー）を入れ換え。
result[item][person] = dic[person][item]
return result
# {ゲーム1:[(ゲーム1との類似度,ゲーム2)]}。類似度の高いゲームを上位3つまで。
def itemSimilar(dic):
result = {}
itemDic = trans(dic) # criticsを[ゲーム:{ユーザー:評価点}]にする。
for item in itemDic:
scores = match3(itemDic,item) #itemと類似度の高い上位3つのゲームをぴっくうｐ。
result[item] = scores
return result
# アイテムベースから推薦を行う（人の類似度じゃなくゲームの類似度から選定ということ）。
def getRecomeItem(dic,p1):
scores = {}
totalSim = {}
# p1に評価されたゲームと評価点をタプルで取る。。
for (item,point) in dic[p1].items():
# このゲームに似ているゲームたち（上位3つ）をループ。
# {ゲーム1:[(類似度,ゲーム2)]}のvalueをループ。
for (similarity,item2) in itemSimilar(dic)[item]:
if item2 in dic[p1]: continue # すでにあるならむし。
# 類似度*評価点
scores.setdefault(item2,0)
scores[item2] += similarity*point
# 類似度の合計
totalSim.setdefault(item2,0)
totalSim[item2] += similarity
# あとは前と同じように。類似度が0の場合は除く（0除算）。
rankings = [(score/totalSim[item],item) for item,score in scores.items() if totalSim[item] != 0]
rankings.sort()
rankings.reverse()
return rankings
print(best(3,getRecomeItem(critics,'椿')))
# ーーーーーーーーアイテムベースのフィルタリング終わりーーーーーーーー
Permalink | コメント(0) | トラックバック(0) | 14:03
2011-11-03
Project Euler 74
問題はこちら。
昨日の92問目と同様の手法を使いました。
まず各桁の重複組み合わせを求めてリスト化します。
次にその各桁の階乗和を計算しそれを続けていき、循環しない項数が60のものだけfilterします。
最後にfilterしたものそれぞれについて順列の数を求めます。
これで求まるはずです。
こちらの方のコードよりほんの少し遅いですが、まあ頭の悪い僕としては上出来です。
自分の超絶低速PCでだいたい4.5秒くらいなのでよしとします。
プログラムが動いた時、まるで自分のことのように喜んでくれる妹がいて幸せです。
import Data.List import Data.Char -- 辞書順重複組み合わせ。 -- comb [0,1,2] 2 で [[0,0],[0,1],[0,2],[1,1],[1,2],[2,2]] comb :: [a] -> Int -> [[a]] comb _ 0 = [[]] comb [] _ = [] comb (x:xs) n = map (x:) (comb (x:xs) (n-1)) ++ comb xs n -- 階乗。かいじょう！ k :: Int -> Int k n = product [1..n] -- リストに一つしか現れない数字を除去し、複数現れる数字をまとめる。 --g [1,1,2,3,3,3,4] で [[1,1],[3,3,3]] g [] = [] g xs | (b == 1) = g (tail xs) | otherwise = a : g (drop b xs) where a = takeWhile (== head xs) xs b = length a -- 普通の順列、あるいは同じものを含む順列の数（ただし先頭に0を含まない）。 -- j [1,2,2,3,3,3] で 60。つまり6!/(2!*3!)を計算してる。 j :: [Int] -> Int j xs = div (k (length xs)) (product $ map (k) $ map (length) $ g xs) -- 0を含む時の順列の数！！！（先頭に0がきてはいけない） -- jj [0,1,2,3] で 18。この場合3*3!を計算してる。 jj :: [Int] -> Int jj xs = sum [ j (delete x xs) | x <- tail (nub xs) ] -- リストの先頭に0がある時は関数h、ないなら関数jで。 h :: [Int] -> Int h xs | (head xs == 0) = jj xs | otherwise = j xs -- 1から999999の数字を桁の重複組み合わせで効率良く出す。これだと19440個だけ調べれば良い。 -- [[1],[2],・・・[0,4,5],・・・[1,3,3,4,5],[9,9,9,9,9]・・・] a :: [[Int]] a = concat [ tail $ comb [0..9] x | x <- [1..6] ] --各桁を階乗して和を出し、それを繰り返していった時の循環しない項数を出す。 -- f' [6,9] [] 1 で 5（循環しない項数） f' :: [Int] -> [Int] -> Int -> Int f' xs ys c | elem a ys = c | otherwise = f' b (ys++[a]) (c+1) where a = sum $ map k xs b = map (digitToInt) (show a) f:: [Int] -> Int f xs = f' xs [] 1 -- 循環しない項数が60ならTrue。 f60 :: [Int] -> Bool f60 xs = f xs == 60 -- できたね！おにいちゃん！やったぁ！ main = print $ sum $ map h $ filter f60 a
Permalink | コメント(2) | トラックバック(0) | 10:12
2011-11-02
Project Euler 92
問題はこちら。
これはかなり前に解いていたのですが、
ちょっと（というかかなりｗ）遅かったので今回高速化してみました。
案外うまくいきました。僕の低速PCでも3秒くらい。
import Data.List import Data.Char -- 辞書順重複組み合わせ。 -- comb [0,1,2] 2 で [[0,0],[0,1],[0,2],[1,1],[1,2],[2,2]] comb :: [a] -> Int -> [[a]] comb _ 0 = [[]] comb [] _ = [] comb (x:xs) n = map (x:) (comb (x:xs) (n-1)) ++ comb xs n -- 各桁の二乗を足し続けて89ならTrue。 -- f [8,5] で True f :: [Int] -> Bool f xs | a == 89 = True | a == 1 = False | otherwise = f b where a = sum $ map (^2) xs b = map (digitToInt) (show a) -- 階乗。かいじょう！ k :: Int -> Int k n = product [1..n] -- リストに一つしか現れない数字を除去し、複数現れる数字をまとめる。 --g [1,1,2,3,3,3,4] で [[1,1],[3,3,3]] g [] = [] g xs | (b == 1) = g (tail xs) | otherwise = a : g (drop b xs) where a = takeWhile (== head xs) xs b = length a -- 普通の順列、あるいは同じものを含む順列の数（ただし先頭に0を含まない）。 -- j [1,2,2,3,3,3] で 60。つまり6!/(2!*3!)を計算してる。 j :: [Int] -> Int j xs = div (k (length xs)) (product $ map (k) $ map (length) $ g xs) -- 0を含む時の順列の数！！！（先頭に0がきてはいけない） -- jj [0,1,2,3] で 18。この場合3*3!を計算してる。 jj :: [Int] -> Int jj xs = sum [ j (delete x xs) | x <- tail (nub xs) ] -- リストの先頭に0がある時は関数jj、ないなら関数jで。 h :: [Int] -> Int h xs | (head xs == 0) = jj xs | otherwise = j xs -- 1から9999999の数字を桁の重複組み合わせで効率良く出す。これだと19440個だけ調べれば良い。 -- [[1],[2],・・・[0,4,5],・・・[1,3,3,4,5],[9,9,9,9,9]・・・] a :: [[Int]] a = concat [ tail $ comb [0..9] x | x <- [1..7] ] --こたえだよ！おにいちゃん！ main = print $ sum $ map h $ filter f a
今回のコードはそこまで圧倒的にくそくないです。たぶん。
Permalink | コメント(0) | トラックバック(0) | 17:18
Project Euler 54
問題はこちら。
バグ取りによって一時死去したけれど、復活して完成。とりあえずうｐ。
変にまわりくどいし、Haskellなのに超長いし（他の人々のと比較して2〜4倍長い）、
しかも型書いてないし（ぉぃ、とんでもなく圧倒的にくそいよお。
disらないでください。すいませんorz
本当にくそいです。
import Data.List import Data.Char import Data.Maybe -- データ構造 ([カードの数値],"スート")。 ehuda = zip ['T','J','Q','K','A'] [10,11,12,13,14] f s | isAlpha a = fromJust (lookup a ehuda) | otherwise = read (take 1 s)::Int where a = head s g ss = (reverse $ sort $ map f ss , map last ss) -- ポーカーの役判定。cardの型は([カードの数字],"スート")。 -- リストの要素が全て同じならTrue。 same xs | (length xs < 2) = True | (head xs == xs !! 1) = same (tail xs) | otherwise = False -- nカードなら(True,[役の数字部分],[数字全体])。以下全てこの型。 nCard (xs,ys) z n | (length xs < n) = (False,[],[]) | same t = (True,t,z) | otherwise = nCard (tail xs,ys) z n where t = take n xs -- ワンペア。 onePair card = nCard card (fst card) 2 -- ツーペア（なんかあれｗ）。 twoPair' card z | same a && same b = (True,a++b,z) | same a && same c = (True,a++c,z) | same d && same e = (True,d++e,z) | otherwise = (False,[],[]) where f = fst card a = take 2 f b = take 2 $ drop 2 f c = drop 3 f d = take 2 $ tail f e = drop 3 f twoPair card = twoPair' card (fst card) -- スリーカード。 threeCard card = nCard card (fst card) 3 -- ストレート。 stright card | f == (reverse $ take 5 $ iterate (+1) (last f)) = (True,f,f) | otherwise = (False,[],[]) where f = fst card -- フラッシュ。 flash card | same (snd card) = (True,f,f) | otherwise = (False,[],[]) where f = fst card --フルハウス。 hullHouse' card z | same a && same b = (True,a++b,z) | same c && same d = (True,c++d,z) | otherwise = (False,[],[]) where f = fst card a = take 3 f b = drop 3 f c = take 2 f d = drop 2 f hullHouse card = hullHouse' card (fst card) -- フォーカード。 fourCard card = nCard card (fst card) 4 -- タプルの第1要素と第2第3要素を返す。 fst3 (a,b,c) = a snth3 (a,b,c) = (b,c) -- ストレートフラッシュ。 strightFlash card | (fst3 (stright card) && fst3 (flash card)) = (True,f,f) | otherwise = (False,[],[]) where f = fst card -- ロイヤルストレートフラッシュ。 loyalStrateFlash card | (fst card == [10..14] && fst3 (flash card)) = (True,[10..14],[10..14]) | otherwise = (False,[],[]) -- カードの役の(得点,([役の数字部分],[数字全体])) poker card | fst3 (loyalStrateFlash card) = (9,snth3 (loyalStrateFlash card)) | fst3 (strightFlash card) = (8,snth3 (strightFlash card)) | fst3 (fourCard card) = (7,snth3 (fourCard card)) | fst3 (hullHouse card) = (6,snth3 (hullHouse card)) | fst3 (flash card) = (5,snth3 (flash card)) | fst3 (stright card) = (4,snth3 (stright card)) | fst3 (threeCard card) = (3,snth3 (threeCard card)) | fst3 (twoPair card) = (2,snth3 (twoPair card)) | fst3 (onePair card) = (1,snth3 (onePair card)) | otherwise = (0,(fst card,fst card)) -- ポーカーのルール。 pokerRule [s1,s2] | (point1 > point2) = 1 | (point1 == point2) && (a > c) = 1 | (point1 == point2) && (a == c) && (b > d) = 1 | otherwise = 0 where [(point1,(a,b)),(point2,(c,d))] = map poker [s1,s2] -- 二人でのポーカープレイ。 pokerPlay [] = 0 pokerPlay ss = pokerRule (take 2 ss) + pokerPlay (drop 2 ss) -- テキストファイルからデータ取得してポーカープレイ。 ff [] = [] ff s = (take 2 s) : ff (drop 3 s) gg [] = [] gg s = (take 5 s) : gg (drop 5 s) main = do a <- readFile "poker.txt" print $ pokerPlay $ map g $ (gg.ff) a
大事なことなのでもう1度言いますが、本当に本当に、くそいです。
Permalink | コメント(0) | トラックバック(0) | 09:30
2011-10-30
Project Euler 35
問題はこちら。
この問題はネットを少々徘徊した限りでは、僕のコードが一番速かったです。
ちょっとうれしい。えへへ。すごいね！おにいちゃん！！！
（とりあえずうｐ。コメント、メモはあとでまた書こ）
import Data.List
import Data.Char
f :: [Int] -> Integer
f xs = read (map intToDigit xs) :: Integer
a :: [Integer]
a = 2:5:(map f $ concat $ map (nub.permutations) $ concat [ comb [1,3,7,9] x | x <- [1..6] ])
g:: Integer -> Integer
g n = read (map (intToDigit) (tail ns ++ [head ns]))::Integer
where ns = map (digitToInt) (show n)
h' :: Integer -> Int -> Bool
h' n c
| (c == length (show n)) = True
| isPrime n = h' (g n) (c+1)
| otherwise = False
h :: Integer -> Bool
h n = h' n 0
main = print $ length $ filter h a
--重複組み合わせ。
comb :: [a] -> Int -> [[a]]
comb _ 0 = [[]]
comb [] _ = []
comb (x:xs) n = map (x:) (comb (x:xs) (n-1)) ++ comb xs n
--素数の無限リスト。
primesList :: [Integer]
primesList = map fromIntegral primesList'
where
primesList' = [2, 3, 5] ++ f 5 7 (drop 2 primesList')
f m s (p : ps) = [n | n <- ns, gcd m n == 1] ++ f (m * p) (p * p) ps
where ns = [x + y | x <- [s, s + 6 .. p * p - 2], y <- [0, 4]]
--素数判定。
isPrime :: Integer -> Bool
isPrime n = (n > 1) && isPrime' primesList
where
isPrime' (x : xs)
| x * x > n = True
| rem n x == 0 = False
| otherwise = isPrime' xs
Permalink | コメント(0) | トラックバック(0) | 15:29
Project Euler 34
問題はこちら。
上限は7桁です。それ以上になるとだめです。
算数できないので最初何故だかよく解らなかったのですが（え、
実際計算してみたらそうでした。
import Data.List import Data.Char -- 階乗を求める。 -- f 5 で 120を返す。 f :: Int -> Int f 0 = 1 f n = product [1..n] -- 重複組み合わせ。リストの数字から重複を許してn個選ぶ。昇順になる。 -- comb [1,2,3] 2 で [[1,1],[1,2],[1,3],[2,2],[2,3],[3,3]] を返す。 comb :: [a] -> Int -> [[a]] comb _ 0 = [[]] comb [] _ = [] comb (x:xs) n = map (x:) (comb (x:xs) (n-1)) ++ comb xs n -- 2ケタ以上の数字を、桁のリストで表す。123は[1,2,3]となる。 -- [[0,0],[0,1],・・・[1,2,3],・・・[1,2,3,4]・・・] a :: [[Int]] a = concat [ comb [0..9] x | x <- [2..7] ] -- 問題の条件に合うならTrue。bはそれに合うものをfilterしたリスト。 -- って言っても結局[[1,4,5],[0,4,5,5,8]]なんだけど。 g :: [Int] -> Bool g xs = (sort $ show $ sum $ map f xs) == (sort $ map intToDigit xs) b :: [[Int]] b = filter g a --ふぇぇ〜。つかれたあ〜。でもあともうちょっとだよ！おにいちゃん！ main = print $ sum $ concat $ map (map f) b
ネットを徘徊した限り、幾人かのスゴイ方のコードより若干遅いしメモリ喰ってますが、
まあ合格の範囲内かと。ってか、自分のPC遅すぎorz
「0.1secで出ました」ってコードを走らせると、自分のPCだと20secくらいかかりますｗ
ちなみにこのコードは0.3secくらいです。うん、合格だね！おにいちゃん！
Permalink | コメント(0) | トラックバック(0) | 15:20
Haskell短期強化週間ひとまず終了！
HaskellでEulerを55題クリアし、なんとか目標達成。つかれたあ。
今後ちょいちょい復習をしたり、時々新しい問題を解いたりするけれど、
とりあえず僕の頭の悪さではここまでが限界なようです。
次はHaskellで何やろうかな。何か作りたいけどなんだか無理そうorz
Permalink | コメント(0) | トラックバック(0) | 12:45
2011-10-29
Project Euler 33
問題はこちら。
Raitoモジュール素晴らしいです。
こんなものがあったんですね（今更…orz
でもこの問題、イマイチ解ってません。
分子/分母のそれぞれの2桁をab/cdと表したとして、
キャンセルされるのって常にbとcだけなんですか？
aとdやbとdがキャンセルされて問題の条件満たすことってないのでしょうか。
そこがよく解らないんですが、なんだかみなさんのコード見てるとそうみたいなので、
もうそんな感じにしちゃいました（ぉぃ。
たぶん自分の頭があまりにも悪くて、
単純な算数の問題に気づいていないだけだと確信していますが…。
うーん…本当によく解らない。
どうしてだろ。
import Data.Ratio --約分した形で与えてくれる。f 49 98 なら 1 % 2 って返ってくる。 f :: Ratio Integer -> Ratio Integer -> Rational f n m = n / m --分母・分子のそれぞれの桁に数字を入れていき、 --問題の条件に合うものだけを抽出。 x :: [Rational] x = [ f a c | a <- [1..9], b <- [a..9], c <- [a+1..9], f (10 * a + b) (10 * b + c) == f a c ] --denominatorは分母を返すんだね！おにいちゃん！ main = print $ denominator $ product x
Permalink | コメント(0) | トラックバック(0) | 15:02
2011-10-26
Project Euler 49
問題はこちら。
これは個人的に素直にできたと思います。
でもネットを徘徊するとやたら複雑なコードなものが多くて、
はてにゃ？と思いました。なんで・・・。
まあ自分のコード、微妙に汎用性みたいなものないですが笑
import Data.List
--今回使う範囲の素数リスト
a :: [Integer]
a = takeWhile (<3340) $ dropWhile (<1000) primesList
--要素数3のリストが全部素数ならTrue。
f :: [Integer] -> Bool
f [x,y,z] = all isPrime [x,y,z]
--要素数3のリストの全要素ををshowしてsortして全て同じならTrue。
g :: [Integer] -> Bool
g [x,y,z] = (\[a,b,c] -> (a==b) && (b==c) && (c==a)) $ map (sort.show) [x,y,z]
--項差3300で項数3の数列のリストを作る。
h :: [Integer] -> [[Integer]]
h [] = []
h (x:xs) = [x,(x+3330),(x+6660)] : h xs
b :: [[Integer]]
b = h a
--要素数3のリストの要素を全てshowして連結させて数値に戻す。
i :: [Integer] -> Integer
i [x,y,z] = read (show x ++ show y ++ show z) :: Integer
--らすとおー。
main = print $ i $ last $ filter g (filter f b)
--素数の無限リスト。
primesList :: [Integer]
primesList = map fromIntegral primesList'
where
primesList' = [2, 3, 5] ++ f 5 7 (drop 2 primesList')
f m s (p : ps) = [n | n <- ns, gcd m n == 1] ++ f (m * p) (p * p) ps
where ns = [x + y | x <- [s, s + 6 .. p * p - 2], y <- [0, 4]]
--素数判定。
isPrime :: Integer -> Bool
isPrime n = (n > 1) && isPrime' primesList
where
isPrime' (x : xs)
| x * x > n = True
| rem n x == 0 = False
| otherwise = isPrime' xs
集中的なHaskellリハビリ疲れたorz
でももうリハビリは終わり。このまま突っ走る！笑
Permalink | コメント(2) | トラックバック(0) | 23:28
Project Euler 46
問題はこちら。
微妙にハマりました。
アプローチを変えたらできたので、まずそちらを。
--奇合成数の無限リスト
a :: [Integer]
a = [ x | x <- [3,5..] , (not.isPrime) x]
--2乗して2かけた数の無限リスト
b :: [Integer]
b = map (\n -> 2*n^2) [1..]
--ちょめちょめする。内包表記の中で内包表記を使ってうんちゃらと。
main = print $ head $ [ x | x <- a , all (not.isPrime) $ takeWhile (>0) [ x-y | y <- b ] ]
--素数の無限リスト。
primesList :: [Integer]
primesList = map fromIntegral primesList'
where
primesList' = [2, 3, 5] ++ f 5 7 (drop 2 primesList')
f m s (p : ps) = [n | n <- ns, gcd m n == 1] ++ f (m * p) (p * p) ps
where ns = [x + y | x <- [s, s + 6 .. p * p - 2], y <- [0, 4]]
--素数判定。
isPrime :: Integer -> Bool
isPrime n = (n > 1) && isPrime' primesList
where
isPrime' (x : xs)
| x * x > n = True
| rem n x == 0 = False
| otherwise = isPrime' xs
で、ハマったのが以下。未だによく解らない。
最後の行がprimesListになった途端だめになる。
他の普通のリストとかだと大丈夫なのになんで…。
--二乗数か判定するよ。 f :: (Floating a) => a -> Bool f n = a !! ((length a)-2) == '.' && last a == '0' where a = show (sqrt n) --なんだかんだでこんな感じにやる。 g n (x:xs) | n < x = False | f ((n-x)/2) = True | otherwise = g n xs h n = g n primesList
混乱してきたのでまた今度。
Permalink | コメント(4) | トラックバック(0) | 23:25
Project Euler 45
問題はこちら。
これは結構簡単でした。
五角数と六角数だけでやれば良いので、
単純に同じものがあるかどうか調べました。
--まずは五角数、六角数の無限リストを作る。 p :: [Integer] p = dropWhile (<= 40755) $ map (\n -> div (n*(3*n-1)) 2) [1..] h :: [Integer] h = dropWhile (<= 40755) $ map (\n -> n*(2*n-1)) [1..] --地道にあるかどうか調べる。 f :: [Integer] -> [Integer] -> Integer f (x:xs) (y:ys) | elem x (takeWhile (<=x) (y:ys)) = x | otherwise = f xs (dropWhile (<=x) (y:ys)) --もうおわり。 main = print $ f p h
でも余裕がないので、関数名とかできとーです。
すいませんorz
Permalink | コメント(0) | トラックバック(0) | 23:18
Project Euler 43
問題はこちら。
結構死んだけどできたわあい！
でも余裕は全然ないので、解説というより完全に自分用メモです。
すいませんorz
import Data.List
--bはaで割り切れるか。
f :: Int -> Int -> Bool
f a b = mod b a == 0
--[[17で割り切れる3桁の数リスト],[11で割り切れる3桁の数リスト]・・・]を作る。
g :: [Int] -> [Int] -> [[Int]]
g _ [] = []
g xs (y:ys) = filter (f y) xs : g xs ys
--ちょっとmapのmapに注意。
a :: [[String]]
a = map (map (show)) $ g [12..987] [17,13,11,7,5,3,2]
--2ケタは先頭に0つけないとだめだあー！
zero :: String -> String
zero c
| length c == 2 = "0" ++ c
| otherwise = c
--というわけでaの完成版b。
b :: [[String]]
b = map (map (zero)) a
--リストの要素が全て違うならTrue。
diff xs
| (length ss == 1) = True
| (head ss == ss !! 1) = False
| otherwise = diff (tail ss)
where ss = sort xs
--h [17で割り切れる3桁の数リスト] [11で割り切れる3桁の数リスト] みたいにしてちょめちょめ。
h :: [String] -> [String] -> [String]
h xs ys = [ [head y] ++ x | x <- xs , y <- ys ,
(tail y == take 2 x) && (diff $ ([head y] ++ x)) ]
--scanl1でさっき作ったbのリストにhを順々に適用（これで出るのかあ）。
c :: [String]
c = last $ scanl1 h b
--最後に足りない数字を見つけて先頭に足す（忘れてた）。
j' :: String -> String -> String
j' (x:xs) (y:ys)
| (x == y) = j' xs ys
| otherwise = [y]
j :: String -> String
j k = j' (sort k) "0123456789" ++ k
--これが最終的に欲しいリスト。
d :: [String]
d = map j c
--さいごなのっ！
main = print $ sum $ map (\n -> read n :: Integer) d
ほんと、これは自分としては結構ハマった問題でした（トホホ。
Permalink | コメント(0) | トラックバック(0) | 23:13
