#!/usr/bin/env python

import fileinput 
import pylab
import matplotlib.pyplot as plt
from matplotlib import rc
from math import *
fig_width_pt = 246.0  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inch
golden_mean = (sqrt(5)-1.0)/2.0         # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt  # width in inches
fig_height = fig_width*golden_mean      # height in inches
fig_size =  [fig_width,fig_height]
params = {'backend': 'ps',
           'axes.labelsize': 20,
           'text.fontsize': 20,
           'legend.fontsize': 20,
           'xtick.labelsize': 20,
           'ytick.labelsize': 20,
           'text.usetex': True,
           'figure.figsize': fig_size}
pylab.rcParams.update(params)

#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
## for Palatino and other serif fonts use:
rc('text', usetex=True)
rc('font', family='serif')


time_s=[]
spot_p=[]
forward_p=[]
b_theta=[]
kappa=[]
for line in fileinput.input():
    items=line[:-1].split()
    spot_p.append(items[0])
    forward_p.append(items[1])
    b_theta.append(items[2])
    time_s.append(items[4])
    kappa.append(items[5])
    print items[0],items[1],items[2],items[3],items[4],items[5]


fig = pylab.figure(1, figsize=(13,10))
plt.xlabel('Time to maturity')
plt.ylabel('Prices')
plt.plot(time_s,spot_p,label='Spot price')
plt.plot(time_s,forward_p,label='Forward price')
plt.legend()
plt.xlim([0, 1.6])
#plt.show()
plt.savefig('price_his.eps')


fig = pylab.figure(2,figsize=(13,10))
plt.xlabel('Time to maturity')
plt.ylabel(r'$\kappa$')
plt.plot(time_s,kappa,label=r'$\kappa$')
plt.legend()
plt.xlim([0, 1.6])
plt.ylim([-10, 5])
#plt.show()
plt.savefig('kappa.eps')


fig = pylab.figure(3,figsize=(13,10))
plt.xlabel('Time to maturity')
plt.ylabel(r'$\bar{\theta}$')
plt.plot(time_s,b_theta,label=r'$\bar{\theta}$')
plt.legend()
#plt.xlim([0, 1.6])
#plt.ylim([-10, 5])
#plt.show()
plt.savefig('b_theta.eps')

 

ブルーベル

日本野鳥の会のメンバーで丹頂鶴を見たこともある英国人と遭遇。奇遇ですね。ブルーベルを見てきた。

Pythonでグラフを描く。
pylab.plot(xx,yy)
pylab.show()

UbuntuでPythonのモジュールをインストールする。

sudo apt-get install python-numpy python-scipy python-matplotlib python-scientific python-stats

 

テートブリテン

素晴らしい晴天(夕方からは雨だったけど）。テート・ブリテンへ。バスで5分乗り換えなしは魅力的な目的地。
財布を持っていくのを忘れてコーヒーも飲めなかった。

#!/usr/bin/python

#Project Euler Problem 27
#Start 01/May/2010
#End   01/May/2010

import math

def is_prm(i):
    if i<0:
        return 0
    else:
        sq_i=math.ceil(math.sqrt(i))
        k=2
        while k<=sq_i:
            if i%k==0:
                return 0
            k+=1
        return 1

def fun(a,b,n):
    return math.pow(n,2)+a*n+b

prm=[]
prm.append(1)
for i in range(2,1000):
    k=is_prm(i)
    if k==1:
        prm.append(i)
print prm
mx=0
for a in range(-1000,1001,1):
    for b in prm:
#        print a,b
        n=0
        while 1:
            k=is_prm(fun(a,b,n))
            if k==1:
                n+=1
            else:
                break
        if n>mx:
            mx=n
            print a,b,n,a*b

 

チェルシー王立病院

チェルシーにあるチェルシー王立病院。ここはチェルシー・フラワー・ショーの会場です。普段は戦争関連の年金受給者の病院です。ひどい表現だな。フラワー・ショーは既に有名過ぎてただのお金集めのイベントというのが私の感想なので行ったことはない。値段も法外にしか思えない。

普段の王立病院は散歩に最適です。隣のアーミー博物館は子供クラスが開かれています。さすがイギリス。

プロジェクト・オイラーは29が解けた。100の100乗を計算するなぞ昔のプログラマーとしてはびびるのだが、多倍長演算が使える言語は良いな。ソートもどうやろうかななんて考える必要が無かった。最後にリストとセットを入れ替えて、リストの長さを出力。ああ、これがウェブ2.0という奴か。おっさんはびっくりだよ。

#!/usr/bin/python

#Project Euler Problem29
#Start 30/Apr/2010
#Start 30/Apr/2010

seq=[]

for i in range(2,101):
    for j in range(2,101):
        seq.append(pow(i,j))



seq1=set(seq)
seq=list(seq1)
seq.sort()
#print seq
print len(seq)

離乳食2日目。今日は舌を出すことを発見したらしい。

 

木曜日はパブへ

Problem 23は45秒に。ついでに40分程かけて26を解いた。やっとレベル1に到達。

朝に電車を見せる。キーという甲高い音は苦手なようだが泣かない。次は買い物。バギーに乗ると大体寝る。
午後は学校へ。
チェコ人のR君は来週から働くらしい。結局博士号はどうなるのだろう？お金の問題は人生の問題。R君とドイツ人のFとその他諸々でパブへ。

 

国際化著しい

朝起きてミルクを飲ませ、風呂にいれ学校へ。11時ごろに着。ジョセフしかいねえ。昼ご飯を食べるついでに再び病院へ行き昨日だし忘れたパスポートを見せる。

午後は全員来た。我が部屋も5人で全員出身が違うというなかなかの国際化だ。ちなみにアジア人は私一人。後はイタリアとフランスとポルトガルとイスラエル。ポルトガル人はもうすぐ帰るらしいけど。

夕方はL君がやってきてビジュアルベーシックの使い方を教えてくれと来た。昔々、まだサラリーマンだった時に使う破目になっていて毎日不快だった。もうあまり覚えてない。

Project Eulerは23が解けた。でも１分以上かかっている。abundantを計算するだけで1分かかるのですけど…

 

Python memo

print sum( int(x) for x in str(pow(2,1000)) )

 

papas and mamas

500桁の数字が10×50の表示でファイルに書かれている。パイソンで500個の要素を持つ配列として読み込みその後の計算に用いる。

そんなシンプルでどこにでもありそうな問題。なのに私は上手でお洒落なコードを書くことができない。

お昼はタイ料理。日が高いうちに大学を出て、オックスフォード・サーカスにあるパパス・アンド・ママスへ。複数形を日本語で書くと全然違う単語を想像させるがこれはパパとママのことで子供、赤ん坊関連のショップである。来週には妻と子がやって来る。バギーを購入しなくてはいけないのだけれど、よく分からない。種類が多すぎる。4ヶ月の子供は横に寝かせるべきなのか？10ヶ月の子供を連れている他の客に話しかけてどんな感じが聞いてみた。バギーって買い換えるのか。その人はバギーを使わなくなるまでに3台くらいは買うとのこと。話しかけた相手が大金持ちなのか、それともそれが普通なのか。今度お姉様に電話で聞こう。

小雨が降ったり止んだりの中を近くのJohn Lewisへ。こちらでも子供用品を見る。こちらのほうが空いているのは値段設定が高いからだろうか。さすがイギリス一の高級デパート（多分）。地下にフードホールがあると書いてあるので覗いてみる。残念ながら日本のデパ地下とは違ってスーパーが入っているだけだ。王室御用達のウェイトローズである。折角なのでrump steakとサラダ、パンを買う。パンが大変美味しい。これはフランス産に違いない。イギリスにこんな美味しいパンがあるわけない。近所のウェイトローズに同じパンがあるかどうかで私の食生活は大きく変わるな。

PicasaからBlogThisするところが壊れている。新しい写真をアップロードすることが出来なくなっている。既にアップしてあるのが終わったら終わりかな。前は出来たのだけどな。何かしたっけ？よく聞く話の一つでウィンドウズが止まった、直してくれ、何したの？何もしてないよ、を思い出した。確かに何もしていない（笑）。chromeを入れたくらいかな。そういえばemacsのフォントがある時を境に（何をしたのかは覚えていない）変わって、これはこれで良いやと受身になっていた辺りからなんとなく悪いのかも。

 

Project Euler

Project Eulerを始めてみた。簡単な方から（正解者数が多い方から）4問だけ解いた。
これまでの所、数学を使う必要が全くない。5問目のProblem 4は少し頭を使うかも。コンピューターに計算させても良いか。主題はPythonの練習。

教務課（に相当する部署）からメールが来ておめーの書類が間違っているから受け取れないぞとのこと。先週もありましたね。大学がそこまで私のために書類をチェックしてくれるなんて涙が出そうです。でもその書類は前回の時に秘書の方と顔を合わせながら作成したから間違っているはずがないのです。と思っていたら早速秘書の方からちょっと教務課がバカだから待っていてとのメールが来た。結局教務課がバカだったようだ。日本に帰りたい。

 

PythonからCを呼ぶ

CからPythonを呼ぶこともできるが、今回はいらない。

---------------------
import ctypes

共有ライブラリをロードしておく。今回初めて共有ライブラリを作ってみた。それまでは.aしか知らなかった。時代は変わる、ボブディラン。

libmain=ctypes.CDLL("./test.so")

共有ライブラリの中にある、関数をPythonにあてる。

ptest1=libmain.test1
ptest2=libmain.test2


引数と戻り値のタイプを指定する。
ptest1.argtypes = []
ptest1.restype = ctypes.c_double

ptest2.argtypes = []
ptest2.restype = ctypes.c_double

これはウェブから
def test_test(lst):
cintarray=ctypes.c_double*len(lst)
clist=cintarray(*lst)
ptest1(clist,len(lst))
return [x for x in clist]


どちらでも動く。
print libmain.ptest2(None)
print ptest2(None)

print test_test([0.0,0.0])