「Perl」タグアーカイブ

これからgitを勉強しようとする人がはじめに読むべき1冊

最近、Rスクリプトとかpythonスクリプトの管理にGitを使っています。
勉強するのに困ったのは、日本語の書籍が全くといっていいほどないこと。最初は、ネットの情報をたよりにやっていたのですが、そもそもシステム自体が理路整然としたものであるバージョン管理システムなのに拾い読みしてもなかなか理解が進みません。そんな時に見つけたのが、WEB+DB PRESS Vol.50でのgit特集でした。

IMGP5974

著者

現在のメンテナであるJunio C Hamanoが書かれてます。

素晴らしい点

本特集は全6章に分けられていて、バージョン管理システムの予備知識がない人も読み進められるようになっています。そして、4章までは、個人でgitを使う場合、またはローカルリポジトリを使う場合の説明に割かれていて、残り2章はグループでgitを使う場合、または共用リポジトリを使う場合について説明がなされています。ただ、最後の2章も必ずしもグループでしか使えないという説明ではなくて、一人で使う場合の応用として説明がなされています。つまり、gitは一人だろうがグループだろうが有用なツールだということが、非常にうまく説明されています。

ちょっと感動してしまったこと

特集の中のコラムで、gitは階層的なのかP2Pなのか、といった話が出るのですがgitはどちらにもなり得るシステムであるそうな。それは至極当たり前のことで、トップレベルのリポジトリにも、例えば仮に僕が作ったパッチを反映することもできるわけで、そういった柔軟性がgitが指示される理由の一つなのかも。
ただ、視点を現実に替えてみると会社も役所もすごく階層的なんですよね。社長がいて、部長がいて、みたいな。うーん、なぜなんだろう。gitには階層構造もP2Pもあるのに、社会でP2Pが発達しにくいのってなぜ??その疑問に対するjunio氏なりのヒントがコラムに書いてあって思わず納得してしまいました。ネタバレになるので、これ以上は本誌でどぞ。

まとめ

本特集のおかげで、gitの使い方がだいぶわかりました。いちお今書いているちょこちょこしたものもgithubにホスティングするようにしてます。
myuhe’s auto-complete-acr.el at master – GitHub
今のところ、CUIがメインなのでなれてないとめんどいですが、tortoisGitなるものもできてきてるのでこれからは、かなり普及するかも。

WEB+DB PRESS Vol.50

WEB+DB PRESS Vol.50 WEB+DB PRESS編集部

技術評論社 2009-04-24
売り上げランキング : 7688

おすすめ平均 star
star森田創特集(?)
starperl, PHP, SQL

Amazonで詳しく見る

auto-complete.elを拡張するR用のelisp書いた

※20090608追記:以下に書いているelispはさらに拡張されバージョンアップしています。詳しくは、auto-complete.elを拡張するR用のelispをさらに拡張してみた。をご覧ください。

Rでコードを書くときは、もっぱらESS on emacs(meadow)で、すこぶる便利に使わせてもらっているのですが、唯一不満な点が補完でした。abbrev、dabbrevといちおあることはあるのですが、補完候補が何かコマンドを押さないと出てこないし、しかも補完候補がすごく見にくい。

実は、昔perl書くときに使っていたxyzzyはsiteinit.lにごにょごにょすると、この問題はあっさり解決してしまうわけで。

ミニバッファの補完候補リストをチップヘルプにする

この機能になれていたというのもあって、emacsの補完のしょぼさはかなり苦痛でした。といいながらもESSの便利さに負けてしまい、xyzzyはサブのエディタに降格してしまったわけですが・・・

んで、学会終わって余裕もできたので、たまにはemacsの環境整理でもしてみるかと、Rjpwiki見てたら谷村さんから次のようなコメントが・・・

ESSのメーリングリストにも投げたのですが、便利すぎるauto-complete.elのR用拡張って誰か書いていませんか。rubyとかc++用の拡張のみならず、octave用の拡張も書かれているのに、R用のものがまだありません。私家版ですでに書いている人がいれば、是非公開してみてください。auto-complete-modeが初耳の方はDemoを見てみてください。ESSのR-modeで使ってみましたが、便利すぎます。後はあらかじめキーワードなどをリストアップした拡張を作れば、さらに便利になるのですが。。。 — 谷村 2009-03-04 (水) 12:02:57

そもそもauto-complete.el自体を知らなかったので、調べてみるとなんかxyzzyの補完候補リストみたいなものを出してくれるんだそうな。しかも容易に拡張可能ですでに各種言語の拡張やoctaveの拡張まで出てるんだそうな。詳しくは動画デモを見てもらえるとその便利さがわかってもらえると思うです。後、以下のリンクも参考になります。

テキスト入力中に補完候補を自動的に表示してくれる auto-complete.el をリリースしました — ありえるえりあ
情報源による拡張が可能な auto-complete 0.1.0 をリリースしました — ありえるえりあ
auto-complete.el にやてふ拡張を作ってみた – がべこれログ – Garbage Collection Log

お~知らなかった。便利すぎるやないですか。しかし、octaveの拡張elispが出ていて、なぜRが出てないんだろう・・・

というわけで、すでにどなたか書かれているかもしれませんが、R用のauto-complete.el拡張elisp書いてみました。
ソースコードはこちらから。今回は初めてgithubにホスティングしてみました。
ちなみにスクリプトを書くにあたってauto-complete-octave.elを参考(というかほとんどまんま)をさせてもらいました。thanks Yen-Chin,Lee!!

使い方
1.auto-complete.elとauto-complete-acr.elをロードパスが通ったフォルダに起きます。

2..emacsに次の3行を追加します。

[cpp]
(require ‘auto-complete)
(global-auto-complete-mode t)
(require ‘auto-complete-acr)
[/cpp]

4.emacsでESSを起動します。
5.M-x auto-complete-modeします。

その他
・ESSモードとiESSモードのhookに引っかけています。

ソースコード見ていただくとわかるのですが、ただRの関数をリストにしてぶち込んでいるだけです。芸がなくてごめんなさい。
関数もnamespace:statsのものを取り込んだだけで、まだかなり不完全です。最初はNAMESPACEファイルを動的に読み込んでリスト化させようかと試みましたが、そこまでの芸当はできませんでした(当方lisp歴3日)。
手力でNAMESPACEファイルからリストを生成しています。NAMESPACEからlispのリスト生成に使ったpythonスクリプトも恥ずかしながらgithubに置いときます。

[python]
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# NAMESPACEにある関数名をlispで使えるように変換させます。
# NAMESPACEのexpoort()の括弧内にある関数名をコピーしてファイルを指定してから使ってください。

import re
mer = []
f = open("D:RscriptessNAMESPACE")
for line in f:
line = line.strip()
mer.append(line)
f.close()
out = "".join(mer)
out = re.sub(" ","",out)
out = re.sub(‘"’,"",out)
out = re.sub(",","" "",out)
out = re.sub("^""",’"’,out)
out = re.sub("$""",’"’,out)
[/python]

ということで、このコードはgithubにてhostして凄腕lisperの方のさらなるハックを期待したい・・・なる他力本願なことを考えています。
改良したelispを作られた方いらっしゃったら、是非githubのアカウントを取ってidを教えてください。collborateしてbrunchできるようにします。
まだまだ知名度の低いgithubですが、使い出すと面白いですのでガンガンbrunch、pushしてもらえると凄くうれしいです。

後、答えられるか甚だ不安ですが、動かんぞ!!とかこんな馬鹿な書き方あるか!!みたいなお叱りはコメントでお願いします。

python本読んでみた

来年からの新しいお仕事に向けて、そしてArcGISをごにょごにょするためのワザ取得に向けて、python本読んでみました。

 とりあえず、今回買ったのは4冊。うち1冊は洋書。pythonは、日本語リソースが少ないとか聞きますが、どうして結構あるもんです。ま、バージョンがちと古かったりもするけど。
 いきなり、4冊読破はつらいので、とりあえず「みんなのpython」あたりから読み始める。
 で、斜め読みして気づいたのは、とりあえずシンプル。Perlとかだと標準で入ってる正規表現もpythonだとモジュール扱いなんですな。
 それと、Rと比較するとなかなかに興味深い。

 Rと似ている点
 あまり構文糖がない。これはPerlと全く異なる訳で、誰が書いても似たようなスクリプトになるってこと。人によってこれを長所ととるか、短所ととるかは大きく分かれると思いますが、少なくとも人が書いたスクリプトをよく使う、または自分の過去のスクリプト資産に強く依存しているってな人はすごく便利な特長だと思う。
 大昔に自分が書いたPerlのスクリプトとかたまに見直すと、まぢで読めない。その時の自分の知識が浅はかで後で読むなんてこと考えてなかったっていうのもありますが、それさっ引いても、よめん。
 その点、Rで書いたやつとかは、だいぶ読めます。それと人が書いたやつもおおかたわかる。そういう点は言語体系としていかしますな。

 Rと違う点
 データ型がまぢ少ない。これは、多分極力シンプルに作りたいという意図と、純粋オブジェクト指向言語として、面倒なデータ扱いそうなときは、オブジェクト指向的に扱ってくれよっていうことなんだろうか。
 一方で、Rはベクトルに始まり、マトリックス、リスト、データフレームだの定義済のデータ型がたくさんあって、それらを柔軟に使いこなすのは結構大変だったりするわけで。逆にオブジェクト指向については、plot()だのsummary()だのオブジェクト指向的ポリモーフィズムをうまく使った便利な機能があるかわりに、クラスの定義だの継承だのあんまり厳格になってなさそうな気が。どだい、Rってのは、統計解析環境なわけで、開発環境ではないっていう割り切りがそういう風に進化させたのかもしれんですな。ま、エンドユーザーとしては、そこまでは必要ないと言えばない。

 いろいろと特徴があるわけですが、共存させるとすんごく面白いかもしれない。つまり、pythonをhubにして、GISとRの幸せな出会いをプロデュースとか。
 試しになんか作ってみよう。

 

結構厚い・・・ま、全部読むわけぢゃないし。

そして、賽・・・じゃなくてスクリプトは投げられた

昨日、作ったスクリプトの実戦投入。
 まず、ベクタのデータをラスタライズして、テキストに変換。最初は、メッシュサイズを10mにしようかと思ってましたが、テキストに吐き出したら、メモリが足らんことになりそうだったので、20mに。

 ここまでは、不本意ながらArcGISにお任せ。特段文句も言わず成功。

 そして、テキストをRに投げてみる。まずは、サンプルモードだったスクリプトを実戦モードにカスタム。
  
 スクリプト投げ込む。完了。

 何の文句も言わない!!
 偉い。偉すぎるぞ、僕スクリプト!!
 メモリもどうやらだいぢょぶそうなので、メイン妖精さんに計算お願いして帰宅。

 「Rプログラミング&グラフィックス」によると、Rはメモリリークの心配がないんだそうな。なぜなのかは知らないですけどね。ただ、perlの時のような理不尽な失敗は多分ないということだろう。
 
 まずは、無事に計算が終わってくれることを願うばかりですが、興味としては、どれくらい時間がかかるかな~ってことか。
 メモリも贅沢に使って、高速に計算させることに特化してますからね。
 
 これで、メモリが足らないと言われたら、データをチビチビ切り刻みながら計算させるしかないんだけど、これは速度とトレードオフでだいぶ遅くなりそう。
 Rjpwikiによると、チビチビ読み込み、計算・・・・を繰り返すより、全体のオブジェクトをどかんと読み込んでおいて、その中で計算させる方が高速なんだそうな。
 ということで、今回も全部のラスタデータをがつんとマトリックスオブジェクトにして、それを読み込んだまま、出力用のオブジェクトも先に読み込んでおくという、とても贅沢な仕様にしているので、だいぶ早く計算してくれると思うんだが・・・

 ま、明日の朝を待ちますか。

 

   

   

 

ラスタ演算もRでやってみる

結局、振り出しに戻った解析作業。
 とりあえず、データに戻って、ぶち込む変数を再度計算させようとしているわけですが、それにはラスタデータの集計が必要です。

 といっても、あんな腐れArcGIS使ってたら精神衛生上悪くってしょうがない。

 ということで、ラスタをテキストで吐き出して、そのデータでラスタ演算させるRスクリプト書いてみました。

 やっていることは、ターゲットセルから円形のバッファを発生させて、バッファ内の特定セル値をカウントさせるという、ラスタ演算入門編にでも出てきそうなネタ。
 んが、こんな単純処理でも腐れGISソフトってば、理不尽なエラー吐きまくりやがるのですな。多分。

 それに比べれば、Rで書けば、ミスしても「なんて僕はお馬鹿さんなんだろう」と自責の念に駆られるだけですむので気が楽。
 しかも、とても親切なエラーメッセージもついてくるので、見通しも立てやすいです。

 ってわけで、書いたスクリプトはこんな感じ。

in.raster <- sample(1:10,1000000,replace=TRUE)
dim(in.raster) <- c(1000,1000)

buffer.length <- 5

x.square.buffer <- rep((buffer.length *-1):buffer.length,times=(buffer.length*2+1))
y.square.buffer <- rep(c((buffer.length *-1):buffer.length),each=11)
distance.buffer <- sqrt(x.square.buffer^2 + y.square.buffer^2)

fun <- function(x,y){
  if(x <= buffer.length){
    y
  }
}

x.buffer <- unlist(mapply(fun,distance.buffer,x.square.buffer))
y.buffer <- unlist(mapply(fun,distance.buffer,y.square.buffer))

out.raster <- in.raster * 0

for(out.y in (buffer.length + 1):(nrow(in.raster) – buffer.length)){
  for(out.x in  (buffer.length + 1):(nrow(in.raster) – buffer.length)){
    check.vector <- in.raster[cbind(x.buffer + out.x,y.buffer + out.y)]
    out.raster[out.x,out.y] <-length(check.vector[check.vector == 5])
  }
}

細かくチェックはしとらんが、大枠はいけてるはず。

今回は、シンプルを心がけてみました。というのも、今回の処理は結構時間がかかりそうなので、無駄なループとかは極力排除してベクトル演算中心のコーディングです。

最初は、apply系をメイン処理に使おうと考えていたのですが、どうやらapply系は速度が遅いんだそうな。
perlのmapとかもそうみたいだけど、forループのラッパーな関数って、総じて遅いってことなんだろうか。
ちなみに、apply系関数の源流ってlispにあるみたいですな。つってもlispなんてエディタいぢるときくらいにしか見ないけど。

 たしかに、Cとかlispとかいろんな言語の影響を受けているであろうRの言語仕様って、非常に複雑なうえに構文糖的な仕様もあるんだけど、perlなんかと比べても見やすい気が立てやすい・・・気がするだけかも。あ、でも結構過去の遺産をうまく使えてたりするんですよね。

 さあ、後はほんちゃんのデータにこのスクリプトを突っ込むだけだ。ターゲットは6億セルだが、どんなもんだか。
 昔、perlで同じようなスクリプト書いたときは、謎のメモリリーク起こして、ODさんに泣きついたもんですが、今回は頼れるのは自分だけ。ま、メモリたらんくなったら、分散処理にするとか、メッシュサイズ大きくするとか、対策はある。

 ので、細かいことはそのとき考えるか。
   

   

ニシキヘビの餌食

 最初は、調査データを整理するために22歳頃にかじり始めたエクセルVBAが最初でした。うまく書けてボタンをポチッと押すと、ズババッとスプレッドシートに文字が埋まっていく姿を見たときは感動したものです。

 その後、24になると、エクセルの行数に制限があるという致命的な欠陥からperlを勉強し始めました。竹中さんのサイトとにらめっこし、ウンウンうなりながら勉強したperlは動いている中身は見えませんが、配列やハッシュといったデータ構造、人間がとてもできないような計算をやってくれるというプログラミングそのものにも興味が出始めた頃でしょうか。

 そのご2年のブランクを経て26くらいから本格的にRで計算をさせるようになってきました。Rをプログラミング言語というのかとなると、異論ありそうですが、統計に特化した言語体系は、そもそもの目的のデータ処理に絶妙にはまるわけで、その便利さを享受しつつ現在に至ります。

 と、ここで気づいたのですが、ここのところ新しいプログラミング言語って勉強してない。ま、どの言語も初級レベルの理解に留まってますがね。にしても、このまま、Rがメイン言語というのも、なんだか。
 例えば、次のような会話を想像してみる。

 「どういうやつで書いてるんですか?」

 「Rです。」

 「・・・」

 う~ん、いかん、いかんよ、それでは。そう、クワトロ・バジーナ曰く

 「まだだ、まだ終わらんよ!!」

 というわけで、現実逃避がてら最近新しい言語勉強しているとこです。どうせ仕事でアプリ作ることになるんだし、マイナスにはならんだろ。

 んで、問題は、何語を勉強するかということ。条件は優先度が高い順にこんな感じ。

 ・GUIアプリが楽に書ける。
 ・開発コストが高いのはちょっとつらい。
 ・拡張性が高い。
 ・教科書が豊富。
 ・マルチプラットフォーム。
 ・名前が格好良い。

 で、いろいろと開発言語を調べてみる。

C++
 現存のGUIアプリの開発言語としては、恐らく相当多い部類のはず。つまり、多くのサンプルコード、教科書がある。何やらji-sa氏もこそ勉している模様。ということで、勉強するには文句なしの言語なんだが、何せ高級すぐる。メモリ管理とかポインタとか多分理解できない気がする。ということで、覚えるのに時間かかるし、開発コストが相当高い気がする。

ruby
 国産プログラミング言語。純粋なオブジェクト指向で、開発コストが相当低いそうな。国産であるため、日本語のテキストも豊富。なんだけど、GUIアプリってなさそう。いずれは勉強したい言語ではあるが、今回は却下。

Perl
 ライトウェイト言語の代表格。何せ使ったことがあるだけに開発コストが相当低く済むはず。なんですが、これもあまりGUI開発に関するテキストがない。多分作れるだろうが途中で関門が相当多い気がする。

VS.net
 MSの開発環境パッケージ。とかく、GUI周りのめんどくさいことは全部お任せできそうで、良い感じ。んだけど当然マルチプラットフォームは断念せざるを得なくなるわけで、何よりMSに囲い込まれるというのが生理的に嫌。

JAVA
 何はさておき、マルチプラットフォーム。開発環境もEclipseとかnetbeansとかいろいろあるみたい。ま、ここまでは良いもののJAVAはどちらかというとネットワークよりって気が。スタンドアローンでの拡張性があんまりない気がする。

 とまあ、いろいろと考えてみて決めたのが、python(ニシキヘビ)です。その特徴ですが、まずインタプリタであるということ。つまり、僕のようなへっぽこが、トライアンドエラーでやるのに便利ってとこ。
 次にある程度、GUIアプリ開発に関するノウハウがありそうだということ。guiツールキットもデフォルトのTkinterとか、linuxではおなじみGTK、QTなどなどいろいろありそうです。
 そして、マルチプラットフォームというのもぐぐっと来ますな。MSに囲い込まれずに済みます。
 とどめが豊富なモジュール群と将来的拡張性。例えば、Rpyモジュールを使えば、なんと簡単にRを使役することができるし、gnuplotでグラフなんかも書けちゃいます。それと、ArcGISやQGISと言ったGISはpythonでマクロを書くことができるので、GIS -> python -> R なんていう複合ワザも可能になるかもしんない。

 とりあえず、使い始めたばかりなんですが、開発環境は、guiツールキットはwxpython、ビルダーにwxgladeを使って、起動画面程度はできました。後は毎木調査データを取り込めるスプレッドシートからグラフを作れるようにして、帳票とかも出力できるようにしたいとこ。

 Rと違うのは当然として、Perlなんかと比べてもだいぶ違いますな。テキスト類も買い揃えないと。
 ま、ぼちぼち気長に勉強していきます。

Rのテキスト処理

ご無沙汰です。
Mwクwリ a.k.a プレゼン作りはデザインから派 です。

 支部会で使うデータで地形の凹凸度を出す必要があったのですが、出しようがなかったのでRで書いてみることにしました。

 前に九州全部のDEMとかで組んだときはPerlだったんですが、コードが長くなっちゃうのと、ループのループのループみたいな構造になっちゃうのがちと嫌。

 んでどうにか完成。半日くらいかかっちゃいました。やっぱ、sum()とかで、ベクトルの合計とかができるのは便利。
 コードもだいぶ短く書けますし、無駄なループもないので処理も速い(多分)。
 デメリットは、前にも書きましたが、データ構造の行き来がイマイチわからん。

 例えば今回もデータフレームから取り出した数値をつなげていってベクトルにしても、ベクトルにならないという、とんちのようなトリックを理解できるのに2時間ほどかかってしまいました。

 これはどゆことかいうと、データフレームから数値を取ってきてもそれはベクトルでなく、データフレームになっていて、データフレームをつなげていくと、それはベクトルではなく、リストになってしまうのですね~~
 たいした違いないかと思うと、コレが恐ろしいくらい違うわけで、例えば・・・・

   hogehoge <- hogehoge + 1

 このオブジェクトがベクトルだと、ベクトルの要素全部に1が足されるのですが、リストだとエラーを出しやがります。
というわけで、今回のスクリプトの肝は、リストをベクトルに変換するunlist()。この関数で変換すれば、普通のベクトルのように使えます。
オブジェクト指向型の功罪というか、なんというか・・・
こういった類のトラップは往々にしてあるわけでそゆのを回避しながらコーディングするのがなかなか一苦労なわけでして・・・ま、こういうコードを財産にしていけばそういう悩みも減っていくと信じて、しばらくはRでいろいろ書いてみようかと思います。

##凹凸度計算スクリプト
##周囲10セルの高度を参照して計算
##最初6行はヘッダなのでスキップ
########################################
setwd(“d:/【研究】葉枯れ/地図データ”)
dat<-read.table(“dem.txt”,header=F,skip=6)

##最初の10行はバッファなので、すべてを0にする。
roughness <- dat[FALSE,]
for( i in 1:10){
  col <- dat[i,]
  col <- col*0
  roughness <-rbind(roughness,col)
}

##周囲10セルの値から凹凸度を計算
##行頭10と行末10はバッファなので0にする。
for( j in 11:(nrow(dat)-10)){
v.point.rough <- c(0,0,0,0,0,0,0,0,0,0)
  for( p in 11:(ncol(dat)-10)){
    v.rough <- NULL
    rel.rough <- NULL
    point.rough <- NULL
    v.rough <- c(dat[c(j-10),(p-10):(p+10)],dat[c(j-9),(p-10):(p+10)],dat[c(j-8),(p-10):(p+10)],dat[c(j-7),(p-10):(p+10)],dat[c(j-6),(p-10):(p+10)],dat[c(j-5),(p-10):(p+10)],dat[c(j-4),(p-10):(p+10)],dat[c(j-3),(p-10):(p+10)],dat[c(j-2),(p-10):(p+10)],dat[c(j-1),(p-10):(p+10)],dat[c(j),(p-10):(p-1)],dat[c(j),(p+1):(p+10)],dat[c(j+1),(p-10):(p+10)],dat[c(j+2),(p-10):(p+10)],dat[c(j+3),(p-10):(p+10)],dat[c(j+4),(p-10):(p+10)],dat[c(j+5),(p-10):(p+10)],dat[c(j+6),(p-10):(p+10)],dat[c(j+7),(p-10):(p+10)],dat[c(j+8),(p-10):(p+10)],dat[c(j+9),(p-10):(p+10)],dat[c(j+10),(p-10):(p+10)])
    v.rough <- unlist(v.rough)
    target <-unlist(dat[j,p])
    rel.rough <- v.rough – target
    point.rough <- sum(rel.rough)
    v.point.rough <- c(v.point.rough,point.rough)
  }
buff <- c(0,0,0,0,0,0,0,0,0,0)
v.point.rough <- c(v.point.rough,buff)
roughness <- rbind(roughness,v.point.rough)
}

##最初の10行はバッファなので、すべてを0にする
for( i in 1:10) {
  col <- dat[i,]
  col <- col*0
  roughness <-rbind(roughness,col)
}
 write.table(roughness, “d:/【研究】葉枯れ/地図データ/roughness.txt”, quote=F,col.names=F, append=T,row.names=F)

久しぶりにPerlで書いてみた

 職場の人からGISのテキストデータを編集して欲しいと頼まれ、久しぶりにコーディング。処理としては格納されているx座標とy座標の順番を変えるという単純なものですが、40万行くらいあるので表計算ソフトぢゃ無理。

 最初はRでちゃちゃっと終わらせようかな~と思ってたのですが、もらったテキストが所々にヘッダーのような行が入ってるんですな↓

#Data
#Region 1
#  16
-50852.785000000003 -47664.510000000002
-50837.785000000003 -47627.010000000002
#Region 1
#  15
-51565.279999999999 -46182.625
-51552.779999999999 -46228.875

元々Rは統計用言語なので、こういうのぶっ込むと途中途中のヘッダーを認識してくれなかったりすんので、うまくいきません。ま、強引にすればいくんでしょうけど・・・

ということで、久しぶりにperlで書くことにしてみました。実に4年ぶりの作業。とりあえずActive Perlwーwロwサwクして書き始めようとするものの・・・・完全に忘れてる!!

何$って??@って・・・初歩的なとこの復習を昔自分が書いたコード見ながら勉強。とりあえず、どうにかでけました。正味2時間くらい。思ってたよりさくっといけたです。

#mifの座標のxyを入れ替える。
#途中に文字列あるので、条件分岐の必要あり。
#################################################

open(FILE,”<vector.txt”) or die “open error!!”;
open(OUT,”>>output.mif”) or die “open error”;
while(<FILE>){
    chomp;
    if($_ =~ m/#.*/){
       print OUT “$_
“;
    }
    else{
        @array = split /s+/;
        print OUT “$array[1] $array[0]
“;
    }
}

close FILE;
close OUT;

 ヘッダー要素には頭に#が付くという規則性があるようなので、そこんとこで条件分岐。数秒で処理完了。う~んPerlもなかなか。

 ま、単純な処理だったってのもありますけど。

 改めて思いましたが、やっぱPerlって可読性低いですな。syntactic sugarが豊富っちゅーのもあるんでしょうが、昔自分が書いたコード見ても何書いてるかさっぱりですもの。その点Rは結構厳格というか・・・

 つっても、柔軟なテキスト処理はやっぱPerlの方が上。PerlとRをうまく使い分けてコーディングするのが賢い使い方なんでしょうな。

パーティナイト

 昨日の結婚式はなんだかんだで2時くらいまで飲んでいて、今日はそのダメージで研究室でグダーっと漫画見たり、ネット見たり、後輩の指導に忙しいbunsho君にちょっかい出したり。

 perl教えてたみたいなので、

  漏れ「フフフ。今はRの時代なのだよ。Rだったらsplitしてハッシュとか作らなくてもデータフレームからhogehoge$Hとかで一撃なのだよ。」

   とか偉そうなこと言ってみたら

  bunsho「perl教えたの、あなたでしょうが!!」

と、逆に怒られてしまいました。

 んなことしてたら、外はいつの間にか日暮れ間近・・・・

 

台風過ぎて・・・

 2号館からの夕日。ここから見る夕日はなぜか引きつけられてしまいます。

 実は今日のビッグイベントはこれから。
 ちょくさんの仕事のお得意先さんが主催で親不孝のクラブでイベントがあるのですな!!親不孝のクラブとか何年ぶりですかね・・・

  21時に親不孝の長浜公園に集合。なんか、親不孝も寂しくなっちゃったもんです。こでちゃんは、お通夜に出てきて、遅れて集合。
 23時くらいになると、人がわんさか集まりだしてグッチャグチャ。その後しばらくしてヘアカットショーも始まり(主催は美容室さん)、完全に酸素不足・・・・

 その後、だいぶ人もはけて、良い雰囲気に。なんつっても、曲の選曲が良い!クラブジャズものからブラジリアン、定番レアグルーヴものとかうまくつないでくれます。miceteethとかもかかって、おおっと思ったり。

 これまでの飲み疲れも貯まっていたようで、今回は2時にドロップアウト。もうちょっといたかったですけど、体がね・・・

 でも、久々に楽しかったダンスゥィンナイト!!でした。

Rでプログラムでけた!!

 なかなか難産でしたが、Rで組んでたプログラムどうにか出来ました。

 機能としては、株の大きさからべき乗式使って胸高直径を推定。うんで、その胸高直径からさらに測っていない樹高を今度は拡張相対成長式使って推定するという合わせ技。しかもどっちも非線形回帰。非線形なのでPerlでもできるかもしれんが、多分めんどくさそう。といって、プロットごとに回帰させるのもちょっと。というわけで、Rで組んでみようということになったんですな~

 使うデータセットはDBH(間伐木の分は当然欠損)、根株の直径、樹高(欠損値あり)、それと各プロットのID。

 Rで組むときの大きなメリットは再利用のしやすさでしょうね。今回のやり方も全てデータフレームで管理するようにしてるので、例えば、ヘッダーのIDをそろえると(今回の場合はname DBH GBH,Hとか)他のデータでも応用が簡単にできます。そこは、基本的に1次元のデータ構造しか持たないPerlからすると、大きなメリットです。ま、その分データ構造で挙動が違ってくるので、そこの扱いに慣れるのは大変ですけど・・・

 それと、解析結果の二次利用が簡単ということでしょうか。今回の場合は非線形回帰で得られたパラメータを使って推定に使うというまさしく非線形回帰の結果を二次的に利用しているのですが、パラメータを取り出すのもcoef()で簡単に実現してしまいます。

 もしかすると、PerlでもCPANとかで探すと見つかるのかもしれませんが(ちなみに僕が探したときは見つかりませんでした)、マニアックな解析してそその結果を他の処理に渡す・・・なんてやつはRプログラムの独壇場でしょう。

 ということで、書いたソースを追いときます。ちなみにヘッダーは以下のような感じです。

    name,Subplot,色,No,sp,GBH,G 0.3,H,under.c,備考

 汚いので晒せるもんでもないですが、つっこみ大歓迎ということで・・・・
 

setwd(“d:”)
library(nls)
dat<-read.csv("testdata.csv",header=T) ##毎木データをプロットごとのデータフレームに分ける ############################################################# #name列はプロットID列。summaryの結果をname.matrixにぶっ込む name.matrix <- summary (dat$name) name.flame <- data.frame(name.matrix) ##プロットIDのそれぞれをベクトルにしてname.datにぶっ込む name.dat <- rownames(name.flame) ##各プロットの毎木データをそれぞれのデータフレームに分ける x <- 0 for( i in 1:length(name.dat)){ x <- x+1 each.name <- name.dat[x] each.dataframe <- subset(dat,dat$name == name.dat[x] ) ##assign()でプロットIDをオブジェクト名に当てる。 assign(each.name,each.dataframe) } ##一番外側のループ。各プロットデータに充てたデータフレームをループ ##今回は各プロットでべき乗式、拡張相対成長式を当てはめ、欠損値を補完する。 ############################ y <- 0 for( k in 1:length(name.dat)){ y <- y + 1 ##間伐している木の胸高を推定し、欠損値を補完 ##################################################### res.nls <- nls(GBH~A*G.0.3^B,data=get(name.dat[y]),start=c(A=1,B=1)) ##回帰係数をとる res.coef<-coef(res.nls) ##ラベルをとる res.coef <- as.vector(res.coef) par(ask=TRUE) plot(get(name.dat[y])$GBH~get(name.dat[y])$G.0.3) line <- function(x) {res.coef[1]*x^res.coef[2]} x <- seq(0,,length=100) lines(x,line(x),col="red") x <- 0 mod.GBH.vec<- NULL for( mod.GBH in 1:length(get(name.dat[y])$GBH)){ x <- x+1 if (is.na(get(name.dat[y])$GBH[x])==TRUE){ mod.GBH.vec <- append(mod.GBH.vec,res.coef[1]*get(name.dat[y])$G.0.3[x]^res.coef[2]) } else{ mod.GBH.vec <- append(mod.GBH.vec,get(name.dat[y])$GBH[x]) } } temp.data <- cbind(get(name.dat[y]),mod.GBH.vec) temp.data$DBH <- temp.data$mod.GBH.vec/pi assign(name.dat[y],temp.data) ##樹高未測定木の樹高を推定し、欠損値を補完 ##################################################### sub.forH <- subset(get(name.dat[y]),is.na(H)==FALSE) edit(sub.forH) inv.H <- 1/sub.forH$H resH.nls <- nls(inv.H ~ 1/(A * sub.forH$DBH^B)+ 1/C ,data = sub.forH,start=c(A=1,B=1,C=25)) resH.nls ##回帰係数をとる res.coef<-coef(resH.nls) ##ラベルをとる res.coef <- as.vector(res.coef) res.coef ##散布図と回帰線を描画 par(ask=TRUE) plot(sub.forH$H~sub.forH$DBH,xlim=c(0,50),ylim=c(0,30)) x <- seq(0,50,length=100) line <- function(x) {1/(1/(res.coef[1]*x^res.coef[2])+1/res.coef[3])} lines(x,line(x),col="red") x <- 0 mod.H <- NULL for( j in 1:length(get(name.dat[y])$name)){ x <- x+1 if (is.na(get(name.dat[y])$H[x])==TRUE){ mod.H <- append(mod.H,1/(1/(res.coef[1]*get(name.dat[y])$DBH[x]^res.coef[2])+1/res.coef[3])) } else{ mod.H <- append(mod.H,get(name.dat[y])$H[x]) } } temp.data <- cbind(get(name.dat[y]),mod.H) assign(name.dat[y],temp.data) edit(get(name.dat[y])) }