2014年9月11日木曜日

Perl

最近PDLを使おうかと思ってます。

いや、仕事ではないですが、せっかくperlの勉強をしたいと思っているので。

自由度が高いけど、書き手によってコードのエレガント差が変わる部分が大きい言語は結構好きかもしれません。

学習レベルによってコードが変わるのは、自分の成長度がわかっていいですね。
業務では嫌われる要因でしょうが・・・。

モダンPerlの新版はいつ出るんでしょうか。
結構楽しみにしてます。

あと、Perl入学式に参加して、
「書き方は自由です」
というのは、勇気づけられました。
もちろん、技巧をこらしても良い、ということではなく、
「汚く書いても動くし、技巧をこらしても動くけど、ちゃんと考えて書いてね」
と言われてる気がします。

それはいわば、自然言語に近いのかなと。

文章の流れを考えて、くどく「私は」「私は」と書くのは、
読みづらいし、大事な所で主語が抜けるのは、本当に意味が分からない。

話し言葉と文章では異なりますし。

そういった意味では、使い手のレベルに依存し、
高度なコードも、低度なコードも可能なので、自然言語に似ていると感じます。

さすが、ラリーさんです。

今後も発展して欲しいと思いますし、無くなることはないでしょう。
Linuxにはほとんど入ってますし。

死ぬまでにどれだけ理解出来るか、楽しみです。

2014年6月18日水曜日

カルマンフィルタ

足立先生の「カルマンフィルタの基礎」を勉強させていただきます。
初心者にとっては非常に良い本だと思ってます。

その中で、例題プログラムがあり、MATLABなのですが、
ビンボーなため、持っておりません。

そこでpythonで書いてみました。
最初の例題は、スカラー変数なので、若干弄っております。
コメントはしません。
色々おこられそうなので。。。

本当に勉強になります。

import numpy as np
import pylab as pl
import random

def kf(A,B,Bu,C,Q,R,u,y,xhat,P):
     xhatm = A*xhat+Bu*u
     Am = A;Bm=B;Cm=C;I =1
     Pm = A*P*Am + B*Q*Bm
     G = Pm*C/(C*Pm*Cm+R)
     xhat_new = xhatm + G*(y-C*xhatm)
     P_new = (I-G*C)*Pm
     return xhat_new,P_new,G

A = 1.
b = 1.
c = 1.
Q = 1.
R = 5.
N = 300
x = np.zeros(N)
y = np.zeros(N)
v = np.zeros(N)
w = np.zeros(N)
xhat = np.zeros(N)
t = np.zeros(N)
for i in range(N):
v[i] = random.normalvariate(0,Q)
w[i] = random.normalvariate(0,R)
t[i] = i
#print v[i],w[i]
y[0] = c*x[0] + w[0]
for k in range(1,N):
     x[k] = A*x[k-1]+b*v[k-1]
     y[k] = c*x[k] + w[k]
P = 0
xhat[0] = 0.0
for k in range(1,N):
      xhat[k],P,G = kf(A,b,0,c,Q,R,0,y[k],xhat[k-1],P)

pl.plot(t,y,label='measured')
pl.plot(t,x,label='ture')
pl.plot(t,xhat,label='estimate')
pl.grid(True)
pl.legend()
pl.savefig('kf')

pl.show()


とりあえずうまくいってるようです。


2014年6月15日日曜日

Harten-Yeeをちょっと修正

以前のものを多少修正しました。
もうちょっと美しく出来れば良いんですけど・・・。
1次風上も残してますが、20行くらい少なくなってます。


#advection #################################################
# Harten Yee Upwind ########################################
############################################################
#import#####################################################
from pylab import *
def cal_fai(z):
    eps = 0.0000000001
    if abs(z) > eps:
        return abs(z)
    else:
        return (z*z+eps*eps)/(2*eps)
def cal_sigma(z,dt,dx):
    return 0.5*(cal_fai(z)-dt/dx*z*z)
def minmod2(x,y):
    if abs(x)<=abs(y) and x*y >=0.0:
        return x
    elif abs(x)>abs(y) and x*y >= 0.0:
        return y
    elif x*y<0:
        return 0.0
    else:
        print "Error minmode2"
# Harten Yee ##############################################
class Cell:
    def __init__(self,dxl,u,f,fn,fl,fr):
        self.dxl   = dxl
        self.u     = u
        self.f     = f
        self.fn    = fn
        self.fl    = fl
        self.fr    = fr
        self.df    = 0.0
        self.gi    = 0.0
        self.gamma = 0.0
        self.fai   = 0.0
    def upwind(self,dt,cellup,celldn):
        ultmp = 0.5*(u + cellup.u)
        urtmp = 0.5*(celldn.u + u)
        if self.u > 0:
            self.fl = ultmp*cellup.f
            self.fr = urtmp*self.f
        else:
            self.fl = ultmp*self.f
            self.fr = urtmp*celldn.f
    def caldf(self,celldn):
        self.df = celldn.f - self.f
    def calgi(self,cellup):
        self.gi = minmod2(self.df,cellup.df)
    def calgamma(self,dt,celldn):
        if self.df != 0:
            self.gamma = cal_sigma(self.u,dt,self.dxl)*(celldn.gi - self.gi)/self.df
        else:
            self.gamma = 0.0
    def calfai(self,celldn):
        self.fai = cal_sigma(self.u,dt,self.dxl)*(self.gi+celldn.gi)-cal_fai(self.u+self.gamma)*self.df
    def yee(self,dt,cellup,celldn):
        csdn = 0.5*(celldn.u*celldn.f+self.u*self.f)
        csup = 0.5*(self.u*self.f + cellup.u*cellup.f)
        cfdn = 0.5*self.fai
        cfup = 0.5*cellup.fai
        self.fr = csdn+cfdn
        self.fl = csup+cfup
    def euler(self,dt):
        self.fn = self.f - dt/self.dxl*(self.fr - self.fl)
    def shift(self):
        self.f = self.fn
#set parameters ############################################
in_ff = open('test.txt')
imx=100
dx =  1.0
dt =  0.2
u  =  0.2
print "-"*40
print "u  =",u
print "dx =",dx
print "dt =",dt
print "-"*40
############################################################
#initial condition ################################
cells = []
fini = []
i = 0
for x in in_ff:
    print i
    a = x.split("\t")
    for k in range(0,1):
        print a[k]
    cells = cells + [Cell(float(dx),float(u),float(a[1]),float(a[1]),0.0,0.0)]
    fini = fini + [float(a[1])]
    i=i+1
imx=i
print "imx=",imx
print "*"*40
print "xl,f"
for i in range(0,imx):
    print i,"  ",cells[i].f,cells[i].u
# Time integration ################################
for n in range(0,1000):
    #Harten Yee##############################
    for i in range(0,imx-1):
        cells[i].caldf(cells[i+1])
    for i in range(2,imx-1):
        cells[i].calgi(cells[i-1])
    for i in range(3,imx-2):
        cells[i].calgamma(dt,cells[i+1])
    for i in range(2,imx-1):
        cells[i].calfai(cells[i+1])
    for i in range(1,imx-3):
        cells[i].yee(dt,cells[i-1],cells[i+1])
    for i in range(2,imx-2):
        cells[i].euler(dt)
    #shift#########################################
    for i in range(2,imx-2):
        cells[i].shift()
# Exact ###########################################
fout = [0]
fex = [0]
xl  = [0]
for i in range(0,imx-1):
    if i>=49 and i<=59:
        fex = fex +[1.0]
    else:
        fex = fex +[0.0]
    fout = fout + [float(cells[i].f)]
    xl = xl + [float(i)]
#output############################################
ff=open("output.dat","w")
for i in range(0,imx-1):
    ff.write(str(i))
    ff.write("    ")
    ff.write(str(fout[i]))
    ff.write("\n")
ff.close()
#####################graph#########################
y1=fout
y2=fini
y3=fex
print len(fout),len(fini),len(fex)
plot(xl,y1,label = 'Harten-Yee')
plot(xl,y2,label = 'INITIAL')
plot(xl,y3,label = 'EXACT')
xlabel('x')
ylabel('f')
title('Hatren-Yee')
axis([-0,100.0,-0.2,1.2])
legend()
grid(True)
savefig('Harten-Yee')
show()
###################################################


2014年5月11日日曜日

Java

今、仕事でJavaを使ってます。

正直、FortranやC的な書き方しか出来てません・・・。

今更ですが、自分自身オブジェクト指向が出来ていないと思ってます。

結城先生のJavaの本を読んで、多少わかった気にはなってますが、
まだまだだなと。

実装自体は、数値解析的な所かつライブラリを使う所なので難しい部分はありますが・・・。


2014年3月22日土曜日

Fortran

新しい会社では全く使わなくなったFortran+本の話です。

「数値計算のためのFrotran90/95プログラム入門」は以前の会社に入って良く読みました。京都大学の牛島先生の本です。

特定の分野では77は使われています。
ゼロから書く場合は90/95へのシフトがなされているようです。

 大学の時は77が多くて、
・Windowsかつ小さなプログラムならfcpad+salfordで解く
・Macならg77で解く
だった気がします。

Fortranは科学技術計算に特化していますが、文法は分っても、「自分が書きたいプログラム」の書き方の参考になるものは少なかった気がします。77だとそんなこと気にすることもないですが・・・。common文は使わない、とかimplicit realは使わないくらいでしょうか。

当時、90/95の書籍はあったのですが、流体の計算をする場合の書き方みたいなものを提示してある書籍かつインストール方法やgnuplotを使った可視化まで書いてある本はあまりなかった記憶があります。この本は自習用として良いのではないかとすぐ購入した覚えがあります。

結果、モジュールやサブルーチンの使い方、プログラム規模に応じた使い分け、浅水流方程式の特性曲線法を使ったプログラムまで掲載されており、手垢塗れになりました。付箋もいっぱい貼ってあります。

77はGOTO文を使わないとLOOPを抜け出せない、文番号をつかったFORMAT文、1行72文字まで等の制限がありますが、90/95は自由形式になり、それなりに使いやすくなってます。配列も使いやすくなってます。

Fortranを使う人間は、あまりプログラムがメインでない技術者が多く、文法とか書き方にこだわりすぎても仕方ありません。Cだとメモリに気を使わないとダメですし。
その点、Fortranはちょうど良いかと思います。フリーのライブラリも多くありますし。

フリーのコンパイラならgfortran、商用ならインテルコンパイラやPGIコンパイラがあります。Linuxかつ非商用目的なら、インテルコンパイラはフリーで使えます。

あとFortran一番の特徴は計算が速いことです。これが利用される目的として一番大きいですかね。 時間の数値積分のための繰り返し回数が多くなると、計算終わるまでに時間がかかります。計算が速いのは重要です。




2014年2月10日月曜日

関数プログラミング

プログラミングの書き方はいくつかあると思いますが、
私のこれまでの感覚だと、手続き型とオブジェクト指向型がメインだったように思います。

最近の書籍、雑誌、ネットでは「関数プログラミング」というキーワードが多くあるとおもいます。

私が関数プログラミングに触れ始めたのはEmacsからです。Emacsを使いこなそうとすると、Elispを使わなければなりません。いや、使わなくてもEmacsは使えるんですが、これがないなら他のエディタを使った方が良いです。

Elispを見て、それまでの命令型と全く違って面食らいました・・・。
かっこまみれは良いとして、関数型の思想を勉強してないので、頭が大混乱・・・。

そこから関数型言語としてCommonLispとか、Schemeとかあるんだなと。
ハッカーと画家も買って読みました。読み物として面白いです。
が、その威力が未だに見えてきません。変数書き換え不可!なぜ?

ある程度大規模な科学技術計算(私は流体解析メイン)はFortran、C/C++で書かれることが多く、空間離散点が数百万、時間積分も長時間やる場合は、Intel、PGIのコンパイラが多いようです。

中規模であれば、javaやC#なんかも使われているようです。

小規模であったり、統計解析や定常計算等であれば、Excel、Python、Ruby、R等が使われることが多いようです。

さて、自分のフィールドである科学技術計算に関数プログラミングを使えないかと考えたのですが、なかなか難しいです・・・。

その理由の第一は関数プログラミングを十分理解出来ていないことです。これは本読んで、コーディングすれば少しずつ解決しそうです。

第二に関数プログラミングが適用されている他のフィールドで求められている機能について理解出来ていないことです。これはなかなか難しそうですが、第一の理由を解決する過程で何か答えの断片を得られそうです。

その上で、第一、第二を踏まえて科学技術計算にどのように適用するのかを考えようと思います。

一つの機能としては、並列化への適用は可能性がありそうです。科学技術計算をやる上で、「解析時間の短縮」は大きなウェートです。
可読性向上については、厳しいような気がします。計算屋は命令型しか勉強しないですし、オブジェクト指向すら知らない人が多いです。実際、オブジェクト指向で書くと、遅くなることが多いです。
あとは、保守管理、過去の遺産の利用等々、いくつかの側面についても検討する必要があると思います。

ちなみに、関数プログラミングについては、Twitterで、おそらくその道のプロ(NY先生)であろう方からの返信があって、考え始めました。資料もご提示いただき、少し頭が晴れてきました。このページを見ているとは思いませんが、お礼申し上げます。ありがとうございます。

また、考えがまとまったら書こうと思います。



2014年2月9日日曜日

KKスキーム

河村・桑原スキームについてです。

1984年のLESに関する論文だったかと思います。

数値解析のうち、移流項の解法は以下の流れで勉強しました(私の場合)。

・中央差分(必ず発散)
・1次精度風上差分(解けるけど解が鈍る)
→中央差分の方が精度が良さそうだが解けない。なぜなら移流項(対流項)は移流速度に乗って情報が上流から下流に伝搬するからである。1次精度風上差分の場合、意図しない所で、拡散的な項(数値粘性)が入っているため解が鈍る。

・LaxWendroff法(2次精度だが振動が入る)
・MacCormack法(2次精度だが振動が入る)
→1次精度風上差分では数値粘性が入っていたため、安定的に解けた。同様の考え方で人工粘性を入れることで、振動を小さくする。

・KKスキーム
→3次精度を実現するために4次精度の中央差分+4次の拡散項を入れる。振動は残るが、河村先生の論文以降、実現象解明へ適用されてきた実績がある。

以降、QUICK、TVD(MUSCL系、Non-MUSCL系)、CIP(オリジナル、有理関数補間、CSL系)の勉強をしました。個人的にはCIP-CSLR1が好きですが、既存のコードの書き換えは結構大変です。マルチモーメントなので・・・。
その後は、小松先生の6Point Scheme、牛島先生のQSIスキーム、ENO、WENO等もかじってます。

今回、気が向いて、KKスキームについてちょっと勉強し直しました。
ベンチマークで使われる2D-Cavity Flowに適用しようと考えてます。
非圧縮性流体の解析なので、MAC系の解法を使って、移流項にKKスキームをと・・。

KKスキームは非保存系のスキームでした。ということで、ある論文の方法を使って保存系に変形しました。線形移流方程式に対して、保存系、非保存系で解いた結果、同様の結果を得られたので、適用は出来そうです。

うまくいったら、結果を載せようと思います。

2014年2月4日火曜日

xojo



xojoの2月の課題があったのでやってみた。
全然きれいじゃないし分ってないことを認識しているが、
とりあえずやってみみた。

下記ページを参考にさせていただいた。

http://www.geocities.co.jp/SiliconValley/2413/tips.html#16

なんかプリントスクリーンが巧くいかない・・・。

2014年1月24日金曜日

河川のシミュレーション!

河村先生の書籍です。
現在は絶版になっております。

大学の時に購入しました。
解析の基本的な説明の後、1次元、2次元、3次元の解析について記述されています。

偏微分方程式については移流拡散方程式、ポアソン方程式の解法を、
また一般座標系の説明があります。

河川に焦点を当てた(浅水流方程式に焦点を当てた)書籍は他にないような気がします。

サンプルプログラムもありますが、実用性は低いように思います。
また、流砂量公式については記述がありません。

ただし、勉強し始めた人にとっては有用な本の一つであることは間違いないと感じました。最近は大学レベルでも1からプログラムを書くことは少なくなっているようですし、
あくまでソフロウェアのユーザーであっても、知っておいた方が良いように思います。


個人的にはKKスキームについて記述されていないのが、ちょっと残念です。
日本人が作った数少ない移流項の解法ですからね。
(河村先生のその他の本には記述されています。河川だから省いたのかなと邪推しています)


日本の技術書は薄い本が多いように感じます。
ゼロから勉強しようにも、技術のバックグラウンドや基礎方程式について記述されているだけの本が多いです。基本的には基礎方程式は複雑なことが多く、複雑な方程式を解くには数値計算に頼らざるを得ないことが多いです。そのためには数値解析の基本的事項およびそれを実用的に使う方法まで書くには厚くすべきだと思ってます。
その点、海外の書籍は厚く、読むのが大変ですが、独学に向いています。
日本でもそのような流れにしていただきたいです。

ちなみに河村先生は多くの書籍を出されています。数値計算の基本的な所を記述した教科書になるものが多いです。私はファンなので買ってしまいます・・・。
エクセルVBAを使った解析について書かれているものもありますので、興味のある人は是非!



履歴

年明け一発目です。
せっかくなので、これまでのプログラミング履歴をまとめておこうと思います。
順番は適当です。


  • C/C++
Cは大学のときに最初に学んだ言語。その後、数値解析に目覚めて、CIP法を書いてみて、数式を展開して、思いどおりの結果が出てハマる。研究室に所属して、修論ではC++となるが、当時はあまりオブジェクト指向を理解していなかった。社会人になってからはチマチマ勉強して多少理解が進んだ。業務でも扱うこと(VC++)が多く、とりあえず、一番利用していると思う。


  • Fortran
科学技術計算をやるならFortranは外せない。もちろんゼロから書くならC/C++でも良いが、如何せん、研究室の遺産を利用しようとすると、書けなくても読めなければならなかった。まあ、C/C++に比べれば大したことはない。研究室では77しかなかったが、社会人になってから京都大学の牛島先生の本が出て、90/95を勉強した。当時は配列の動的確保やintent文が使えるなという印象だった。牛島先生の本は90/95の文法だけでなく、科学技術計算のプログラムの書き方について記述してあるのが優れていると思う。


  • python
初めて存在を知ったのは、研究室に入ってから。共有のパソコンにヘビのアイコンがあったのを記憶している。しばらく忘れていたが、社会人になってみんpyをちょくちょく読みながら少しずつ勉強した。高速な計算でなければ使いやすいと思う。数学者が作成した言語というのも好感度アップ。matplotlibによるグラフ作成や、ベクトルやコンター図も簡単にできるのが良い。また、マルチパラダイムなのも良い。特にFortranユーザーにも使いやすいと思う。海外研修時に簡単な温度応力解析のプログラムをこいつで書いた。


  • ruby
日本人が作ったとのことで、たのしいRubyを手に取ったのが始まり。pythonに比べて科学技術計算用途のmoduleは少ないと思う。ただし、オブジェクト指向の勉強としては優れていると思う。Railsは使ったことがない・・・。ごめんなさい。データ整理なら特に問題なく使える。業務で使ったことはない。ただし、日本人ならruby使って発展させた方が良い気がする。黒魔術も使えるし。


  • perl
仕事で海外研修時にO'Reillyの例のやつ(初めてのPerl)を読んでいた。出張中はWindowsしか触れなかったのでPadreを入れていたと記憶している。言語学者が作っているのであまりピンとこなかった。ただし、上記書籍は読んだ方が良いと思う。普通に面白い。業務では科学技術計算用Linuxの設定のついでにLAMPを入れた時に多少触ったくらい。python、ruby、perlについては用途が被るが好みで選択すれば良いと思う。どれが優れているとかは感じなかった。


  • xojo
real basicだったのxojoになった。勉強会も開かれているようだ。もちろん会社では一切使わない。移流方程式を解いて、図化する程度のプログラムは書いてみた。これも海外出張時に多少勉強した。日本語の情報は少ないが、英語の情報はたくさんあるし、リファレンスやサンプルはしっかりしている。ただ、サードパーティなので一抹の不安はある。Linux、Mac、Windowsの環境で使えるし、ネイティブコードを吐き出せるのは良い所。

  • VBA
いや、プログラム言語じゃない!と言われると思いますが、ごめんなさい。会社ではデータ整理やちょっとした統計解析なんかはVBAでやることが多いです。VBAが優れている点は、Excelがインストールしてあれば細かい設定は関係無しにみんな使える所です。通信関係は触ったことはないです。


  • VB.NET
メールを読み込んで添付ファイルを抽出するプログラムを業務で書いた。好んで使った訳でなく、遺産がVBだったから使ったまで。文法については「Must Inherits」とかカッコいいと思った。6.0から.NETで強烈に変化してしまったため、ユーザーは減ってるらしい。私の周り?そもそもプログラム書くことに情熱を持ってるやつなんていないっす。


  • java
javaによるCIP法の本が出てちょっとだけ勉強した。解析結果を表示しながら解析を進めるというもので、結構面白かった。ただ、Fortran的な書き方しかしてないので、javaを勉強した気はしなかった。ほんとはちゃんと勉強したい。Androidアプリも弄ることが出来そうだし。。。


  • ELisp
Emacsを使い始めると必然的に使わなきゃいけない。今だに使えるようになったとは思わないが、Emacsは楽しい。C/C++から入った私としては、文法に面食らった。何かに使うわけでなくEmacsのカスタマイズのみである。なおWindows環境でのEmacsライクなエディタにxyzzyがあるがこちらはCommon Lispでカスタマイズする。


  • AWK
読んだ本はプログラミング言語AWK。スクリプト言語はsh、sed、awk等が基本的な所にあると思う。利用としてはテキスト処理しかないが、結構簡単に書けるしパイプラインで繋いで処理していく過程は楽しい。完全な主観だが、ワンラインならはやりのスクリプトよりこっちの方が良いかもしれない。