Project Euler 28, 29

28番
数字を渦巻状に配置して対角に来る数字の和を求める問題

#!/usr/bin/env python
# -*- coding: utf-8 -*-
summary = 0
for i in range(1,1001+1,2):
	summary += 4*i*i - 6*i + 6
print summary-3

29番
a^bをa,bそれぞれ2~100に変えた時何個違う数字ができるか数える問題

#!/usr/bin/env python
# -*- coding: utf-8 -*-
ans = []
for i in range(2,101):
	for j in range(2,101):
		ans.append(i**j)
print len(set(ans))

どっちも簡単だった

Project Euler 26, 27

26番
最長の循環小数を求める問題

#!/usr/bin/env python
# -*- coding: utf-8 -*-
max_ite = 0
max_d = 0
for i in range(2,1000):
	iterate = 0
	rest = 1
	rest_list = []
	while 1:
		rest = rest * 10 % i
		iterate += 1
		if rest == 0 or rest in rest_list:
			break
		rest_list.append(rest)
	if iterate > max_ite:
		max_ite = iterate
		max_d = i

print max_d

27番
二次式で素数が続く際の係数を求める問題
aは奇数,bは素数という条件のもとブルートフォースで解いたけどなんかもっといい方法ないかなあ

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#aは奇数,bは素数
def f(n,a,b):
	return n * n + a * n + b

def issosu(x):
	for i in range(3,x,2):
		if x % i == 0:
			return False
	return True

def make_b(n):
	b_list = range(2,n)
	b = []
	while b_list != []:
		sosu = b_list[0]
		b.append(b_list[0])
		sub_list = []
		for i in b_list:
			if i % sosu != 0:
				sub_list.append(i)
		b_list = sub_list[:]
	return b

max_len = 0
b_list = make_b(1000)
max_a = 0
max_b = 0
for a in range(-999,1001,2):
	for b in b_list:
		n = 1
		while f(n,a,b) > 0 and issosu(f(n,a,b)):
			n += 1
		if n > max_len:
			max_a = a
			max_b = b
			max_len = n
print max_a*max_b

Project Euler 24, 25

24番
0~9の順列の100万番目を求める問題

#!/usr/bin/env python
# -*- coding: utf-8 -*-

num_list = range(0,10)
#0から数え始める
Order = 1000000
Order -= 1
ans = 0
for i in range(0,10):
	tmp = 1
	for j in range(1,10 - i):
		tmp *= j
	ans = ans * 10 + num_list[Order / tmp]
	num_list.pop(Order / tmp)
	Order = Order % tmp
print ans

25番
1000桁になるフィボナッチ数列の番目を求める問題

#!/usr/bin/env python
# -*- coding: utf-8 -*-

Fn = [1]
Fn_1 = [1]

iterate = 2
while 1:
	C = 0
	save = Fn[:]
	for i in range(0,len(Fn_1)):
		Fn[i] = Fn[i] + Fn_1[i] + C
		C = Fn[i] / 10
		Fn[i] = Fn[i] % 10
	if len(save) > len(Fn_1):
		Fn[len(Fn)-1] += C
		C = 0
	if C > 0:
		Fn.append(C)
	Fn_1 = save
	iterate += 1
	if len(Fn) >= 1000:
		break

print iterate

リストをコピーするには
a = b
じゃなくて
a = b[:]
でやらないとbの値変えた時aの中身も変化してしまう

変数に1000桁入らないんじゃないかと思ってリストでやったけどそんなことなかったみたい。だから

def f(n):
	x = 1
	y = 1
	z = 2
	while 1:
		x,y,z = x+y,x,z+1
		if len(str(x)) >= n:
			break
	return z

print f(1000)

の方が簡単だし速い

Project Euler 22, 23

今回もProject Euler

22番
テキストファイルの名前を並び替えてスコア計算

pythonでテキスト分解するのは結構簡単
sort()もかなり便利

#!/usr/bin/env python
# -*- coding: utf-8 -*-

line = open('names.txt').read()
itemList = line[:-1].split(',')
newitemList = map(lambda x: x.strip(" \" "),itemList)
newitemList.sort()

score = 0
iterate = 0
for item in newitemList:
	iterate += 1
	subscore = 0
	for j in item:
		subscore += (ord(j) - 64)
	subscore *= iterate
	score += subscore

print score

23番
2つの過剰数の和で表せない整数の和

過剰数かどうかは前回の関数を使いまわした
時間は22秒とちょっとかかってしまった...なんかうまい方法ないかなあ

#!/usr/bin/env python
# -*- coding: utf-8 -*-

def yakusu(num):
	row_num = num
	yakusu_list = []
	sosu = 2
	while num != 1:
		sosu_count = 0
		while num % sosu == 0:
			sosu_count += 1
			num /= sosu
		if sosu_count > 0:
			yakusu_list.append((sosu,sosu_count))
		sosu += 1
	insu_wa = 1
	for i in map(lambda x: (x[0]**(x[1]+1)-1)/(x[0]-1),yakusu_list):
		insu_wa *= i
	return insu_wa - row_num

kajyo = []
for i in range(1,28124):
	if i < yakusu(i):
		kajyo.append(i)

kajyo_wa_list = []
for i in kajyo:
	for j in kajyo:
		if i + j < 28124:
			kajyo_wa_list.append(i+j)

print sum(range(1,28124)) - sum(set(kajyo_wa_list))

Project Euler 20, 21

前々から少しずつ解いてきたProject Euler(https://projecteuler.net/)の自分なりの解答をここに残しておく。

20番
100の階乗の各桁の和を求める問題

#!/usr/bin/env python
# -*- coding: utf-8 -*-
ans = [1]
for i in range(1,101):
	C = 0
	for j in range(0,len(ans)):
		ans[j] = ans[j] * i + C
		C = ans[j] / 10
		ans[j] = ans[j] % 10
	while C > 0:
		ans.append(C%10)
		C = C / 10
print sum(ans)

21番
10000未満の友愛数の和を求める問題

#!/usr/bin/env python
# -*- coding: utf-8 -*-
def yakusu(num):
	row_num = num
	yakusu_list = []
	sosu = 2
	while num != 1:
		sosu_count = 0
		while num % sosu == 0:
			sosu_count += 1
			num /= sosu
		if sosu_count > 0:
			yakusu_list.append((sosu,sosu_count))
		sosu += 1
	insu_wa = 1
	for i in map(lambda x: (x[0]**(x[1]+1)-1)/(x[0]-1),yakusu_list):
		insu_wa *= i
	return insu_wa - row_num

def isyuai(num):
	yuai = False
	if num == yakusu(yakusu(num)) and num != yakusu(num):
		yuai = True
	return yuai

ans = 0
MAX = 10000
for i in range(2,MAX):
	if isyuai(i):
		ans += i

print ans

19番までもいつか残しとこう

独立成分分析

春休み「これならわかる最適化数学」

これなら分かる最適化数学―基礎原理から計算手法まで

を読んで主成分分析(PCA)に興味をもって、それの応用(?)の独立成分分析(ICA)も勉強してみることにした。

参考書は「[詳解] 独立成分分析」

詳解 独立成分分析―信号解析の新しい世界

結構分厚い…
とりあえず14章まで頑張ろう。

とりあえず今日は1章
レビューには基本から書かれているみたいなこと書いてあったけどしょっぱなからムズイ！
もう心折れそう

1章でまず何よりもわかってなかったのは「無相関」と「独立」の違いについて。
これについてはここを読んで分かった。
「相関」は一方の値が決まった時にもう一方の値を推定するに役立つ規則性やら何なりがあるかないかを示している。
「独立」は一方の値が決まった時にもう一方の値に影響を与えるかどうかを示している。

相関⊃独立
という関係みたい

またICAを実現する方法としては2通りあって

1. 非線形無相関化
   • s1とs2が独立ならば任意の非線形変換g(s1)とh(s2)は無相関らしくその性質を使うらしいことしかわからなかった…
   • 非線形関数h,gを選ぶには最尤推定を使うみたい。とりあえず9章まで我慢
2. 極大非ガウス性
   • 線形結合の分散が一定という制約条件のもとで、その極大を求める。こっちはなんとか理解できそうな気がする。

らしい、これからゆっくり読み進めていきたいと思う。

K-means

とある事情でK-meansを実装してみた。

大まかな流れとしては
1. 適当にK個のクラスに割り振る
2. それぞれのクラスの重心を代表点とする
3. それぞれの点に最も近い代表点を割り出しクラス分けする
4. 2と3を代表点が動かなくなるまでやる

以下ソース

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import sys
import numpy as np
import cv2
import matplotlib.pyplot as plt

#代表点の個数
K = 3

#画像のサイズ
LEN = 512
LENGTH = LEN*LEN

#RGB
DIM = 3

#LENGTHのRGB画像
im = cv2.imread('lena.bmp',flags=1)
if im is None:
	print "None"
	exit()
im.resize((LENGTH,DIM))

#代表点と画素値からそれぞれのクラスタ番号を得る
def get_len(repre,x):
	length = sys.maxint #代表点までの距離
	for i in range(0,len(repre)):
		if np.linalg.norm(repre[i]-x) < length:
			length = np.linalg.norm(repre[i]-x)
			num = i
	return num

#重心を計算し、新たな代表点を計算
def get_new(im,num):
	repre = np.zeros(K*DIM).reshape((K,DIM))
	count = np.zeros(K)
	for i in range(0,LENGTH):
			repre[num[i]] += im[i]
			count[num[i]] += 1
	for i in range(0,K):
		if count[i]!= 0:
			repre[i] /= count[i]
		else:
			pass
	return repre

#代表点の誤差を計算
def get_error(repre,new_repre):
	error = 0
	for i in range(0,len(repre)):
		error += np.linalg.norm(repre[i]-new_repre[i])
	return error

#初期クラス
num = np.zeros(LENGTH)
for i in range(0,LENGTH):
		num[i] = np.random.randint(0,K)

#現在代表点と誤差の初期化
repre = np.zeros(K*DIM).reshape((K,DIM))
new_repre = np.zeros(K*DIM).reshape((K,DIM))
error = get_error(repre,new_repre)

#k-meansの処理
while 1:
	#クラスタ番号更新
	repre = get_new(im,num)
	error = get_error(repre,new_repre)
	new_repre = repre
	print "error=",error
	if error < 0.1:
		break
	for i in range(0,LENGTH):
			num[i] = get_len(repre,im[i])

#符号なし8ビット整数を使う
im_kmean = np.zeros(LENGTH*DIM, dtype=np.uint8).reshape((LENGTH,DIM))
for i in range(0,LENGTH):
		im_kmean[i][:] = repre[num[i]]

im_kmean.resize((LEN,LEN,DIM))
cv2.imshow("test",im_kmean)
cv2.imwrite('result.bmp', im_kmean)
cv2.waitKey(0)
cv2.destroyAllWindows()

結果は
入力画像

出力画像(3色に減色)

まあちゃんと動いているみたいだけど遅い…
手順3のそれぞれの点から一番近い代表点を割り出す所がとにかく遅い。
これを解消するには木構造を使用したりすればいいらしい。
他にも早くするために最初の代表点を上手く選んでやるなんて方法もあるとのこと。