2017年
11月
1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30

setup diary

2007|12|
2008|01|02|03|04|05|06|07|08|09|10|11|12|
2009|01|02|03|04|05|06|07|08|09|10|11|12|
2010|01|02|03|04|05|06|07|08|09|10|11|12|
2011|01|02|03|04|05|06|07|08|09|10|11|12|
2012|01|02|03|04|05|06|07|08|10|11|12|
2013|01|02|03|04|05|06|07|08|09|10|11|12|
2014|01|02|03|04|06|08|11|
2015|01|02|03|04|05|06|07|08|10|11|12|
2016|01|02|03|04|05|06|07|08|09|10|11|12|
2017|01|02|03|04|05|06|07|08|09|10|11|12|
2018|01|02|03|04|05|06|07|08|09|10|11|12|
2019|01|02|03|04|05|06|07|08|09|10|11|12|
2020|01|02|03|04|05|06|07|08|09|10|11|12|
2021|01|02|03|04|05|06|07|08|09|10|11|12|
2022|01|02|03|04|05|06|07|08|09|10|11|12|
2023|01|02|03|04|05|06|07|08|09|10|11|

2017-11-23 Rで式量計算

_ 元素記号処理と原子量抽出

最近は必要なデータ処理の多くの部分をrubyでやっている。しかし、様々なデータ処理を行う上で、rubyには二つの欠点があると思う。一つは配列同士の四則演算などが面倒である点で、もう一つはすぐにグラフが書けない点である。配列にもいくつかの演算が定義されているが、+は配列を結合だし、*は繰り返しだし、通常の四則演算では無い。mapをつかってゴリゴリ書くか、Vectorを変換すれば、この問題は一部解決する。しかし、グラフを書くことは出来ない。というわけで、今回はRでデータ処理を行うことにした。 Rのスキルも少しは上がってきたので、主要なデータ処理の部分は無事に書き終えた。計算には式量の値が必要なのだが、式量は手で入力している。できれば化学式から式量を自動で計算するようにしたいというのが、次の課題である。rubyを使って式量を計算するプログラムは随分前に書いて、よく使っている。これをRで書くに当たって、二つのプロセスが必要になる。一つは化学式から元素記号と数字を探して取り出すところで、もう一つは元素記号を原子量に変換する部分である。

元素記号と数字を探すのには正規表現を使うのだが、gregexprを使うと、これが可能であることが分かった。しかし、gregexprは文字列の場所と長さを返すので、そこから文字列を取り出すためにsubstrを使ったり、長さを取り出すのにattrを使ったりと、文字列を取り出すまでが非常に面倒である。そこで、stringrというlibraryを使うことにした。library(stringr)としてこれを読み込むと、便利な関数が使えるようになる。最終的に使ったのがstr_match_allというもので、マッチした文字列と、括弧内の文字列を取り出してくれる。これで元素記号と数値を容易に取り出すことができた。

あとは、元素記号から原子量に変換すれば良い。原子量はファイルから読み込むが、read.tableで読み込むと、元素記号のベクトルと原子量のベクトルを合わせたdataframeとなる。ハッシュを使えれば、aw["Cl"]のように原子量を取り出せるのだが、Rではそう簡単にはいかない。まず考えた方法は、aw$V2[which(aw$V1=="Cl")]とか、subset(aw,V1=="Cl")$V2として、原子量を取り出す方法である。しかしこれだと、長い上に元素記号のベクトルを使って、原子量のベクトルを一度に取り出すということができなくなってしまう。

そこで考えたのがnamesを使う方法である。

aww<-aw$V2
names(aww)<-aw$V1

として、原子量のベクトルであるawwにnames属性をつけてやると、aww["Cl"]で原子量を取り出すことができる。さらに、aww[c("Na","Cl")]といくつかの原子量を一度に取り出すこともできる。しかし、dataframe以外にawwという余分なベクトルを作らないといけないのが欠点である。

最終的には、dataframeの列に元素記号の名前をつけると良いことが分かった。read.tableでrow.namesオプションを使うと列の名前を取り込むことができ、列の名前を使うとaw["Cl",]で原子量を取り出すことができる。カンマとかの意味がよく分かっていないが、動けば良いだろう。

そして、最終的にできた関数が次のようなものである。

library(stringr)
fw<-function(fml){
 aw<-read.table("formula.dat",row.names=1)
 el<-str_match_all(fml,"([A-Z][a-z]?)([0-9.]*)")[[1]]
 return(sum( aw[el[,2],] * ifelse(el[,3]=="",1,as.numeric(el[,3])) ))
}

fw("H2O")などとすると、無事に原子量が計算できる。rubyに比べたら、長くなったし、書くのに苦労したが、これでまた少しRのスキルが上がっていたら良いな。