作業は、2地点間の距離と角度を求めるというシンプルなもの
ただ、データ数が多く、組み合わせが全部で22,232,390もある。
まずRでプログラミング
特に変なことはしていないと思いますが...
- Tree<-read.csv("TreePlot",header=T)
- IncW3<-read.table("IncW3Th10S3",header=T,sep=" ")
- IncW5<-read.table("IncW5Th10S3",header=T,sep=" ")
- IncW11<-read.table("IncW11Th10S3",header=T,sep=" ")
- Tree<-as.matrix(Tree[,1:5])
- IncW3<-as.matrix(IncW3)
- IncW5<-as.matrix(IncW5)
- IncW11<-as.matrix(IncW11)
- DistW3<-matrix(rep(NA,nrow(Tree)*nrow(IncW3)*7),ncol=7)
- DistW5<-matrix(rep(NA,nrow(Tree)*nrow(IncW5)*7),ncol=7)
- DistW11<-matrix(rep(NA,nrow(Tree)*nrow(IncW11)*7),ncol=7)
- for (i in 1:nrow(Tree)){
- for (j in 1:nrow(IncW3)){
- DistW3[i*j,1:4]<-c(Tree[i,4:5],IncW3[j,1:2])
- DistW3[i*j,5]<-sqrt((Tree[i,4]-IncW3[j,1])^2+(Tree[i,5]-IncW3[j,2])^2)
- DistW3[i*j,6]<-atan((Tree[i,5]-IncW3[j,2])/(Tree[i,4]-IncW3[j,1]))
- DistW3[i*j,7]<-atan((Tree[i,5]-IncW3[j,2])/(Tree[i,4]-IncW3[j,1]))*180/pi
- }
- for (j in 1:nrow(IncW5)){
- DistW5[i*j,1:4]<-c(Tree[i,4:5],IncW5[j,1:2])
- DistW5[i*j,5]<-sqrt((Tree[i,4]-IncW5[j,1])^2+(Tree[i,5]-IncW5[j,2])^2)
- DistW5[i*j,6]<-atan((Tree[i,5]-IncW5[j,2])/(Tree[i,4]-IncW5[j,1]))
- DistW5[i*j,7]<-atan((Tree[i,5]-IncW5[j,2])/(Tree[i,4]-IncW5[j,1]))*180/pi
- }
- for (j in 1:nrow(IncW11)){
- DistW11[i*j,1:4]<-c(Tree[i,4:5],IncW11[j,1:2])
- DistW11[i*j,5]<-sqrt((Tree[i,4]-IncW11[j,1])^2+(Tree[i,5]-IncW11[j,2])^2)
- DistW11[i*j,6]<-atan((Tree[i,5]-IncW11[j,2])/(Tree[i,4]-IncW11[j,1]))
- DistW11[i*j,7]<-atan((Tree[i,5]-IncW11[j,2])/(Tree[i,4]-IncW11[j,1]))*180/pi
- }
- }
結果
ユーザ システム 経過 1522.730 1.510 1528.272
約25分ぐらいですか
次にFortranの場合
コードはこちら
- program DistCalc
- implicit none
- real,dimension(1565,5)::Tree
- real,dimension(7146,3)::IncW3
- real,dimension(4952,3)::IncW5
- real,dimension(2108,3)::IncW11
- real pi
- real DistW3(1565*7146),AngW3(1565*7146)
- real DistW5(1565*4952),AngW5(1565*4952)
- real DistW11(1565*2108),AngW11(1565*2108)
- integer i,j
- open(10, file='Tree')
- open(20, file='IncW3')
- open(30, file='IncW5')
- open(40, file='IncW11')
- read(10,*) ((Tree(i,j),j=1,5),i=1,1565) !iは行、jは列
- read(20,*) ((IncW3(i,j),j=1,3),i=1,7146)
- read(30,*) ((IncW5(i,j),j=1,3),i=1,4952)
- read(40,*) ((IncW11(i,j),j=1,3),i=1,2108)
- do i=1,1565
- do j=1,7146
- DistW3(i*j) = sqrt((Tree(i,4)-IncW3(j,1))**2 + (Tree(i,5)-IncW3(j,2))**2)
- AngW3(i*j) = atan2((Tree(i,5)-IncW3(j,2)),(Tree(i,4)-IncW3(j,1)))
- end do
- end do
- do i=1,1565
- do j=1,4952
- DistW5(i*j) = sqrt((Tree(i,4)-IncW5(j,1))**2 + (Tree(i,5)-IncW5(j,2))**2)
- AngW5(i*j) = atan2((Tree(i,5)-IncW5(j,2)),(Tree(i,4)-IncW5(j,1)))
- end do
- end do
- do i=1,1565
- do j=1,2108
- DistW11(i*j) = sqrt((Tree(i,4)-IncW11(j,1))**2 + (Tree(i,5)-IncW11(j,2))**2)
- AngW11(i*j) = atan2((Tree(i,5)-IncW11(j,2)),(Tree(i,4)-IncW11(j,1)))
- end do
- end do
- pi=4*atan(1.0)
- open(50,file="DistW3")
- do i=1,1565
- do j=1,7146
- write(50,*) Tree(i,4), Tree(i,5), IncW3(j,1), IncW3(j,2),&
- DistW3(i*j), AngW3(i*j), (AngW3(i*j)*180)/pi
- end do
- end do
- open(70,file="DistW5")
- do i=1,1565
- do j=1,4952
- write(70,*) Tree(i,4), Tree(i,5), IncW5(j,1), IncW5(j,2),&
- DistW5(i*j), AngW5(i*j), (AngW5(i*j)*180)/pi
- end do
- end do
- open(80,file="DistW11")
- do i=1,1565
- do j=1,2108
- write(80,*) Tree(i,4), Tree(i,5), IncW11(j,1), IncW11(j,2),&
- DistW11(i*j), AngW11(i*j), (AngW11(i*j)*180)/pi
- end do
- end do
- end program DistCalc
で結果はというと、一目瞭然ですね!!
$ f95 DistCalc.f90 $ ./a.out real 3m20.994s user 3m14.140s sys 0m5.330s
今回の場合には、Fortranのコードを書くだけで結構時間がかかった（まだ初心者なので）ので、Rだけで作業しても良かったかもしれません。
ですが、今後のことを考えてもFortranをお勉強したほうが良いですね〜
できればC++もお勉強したいのですが...
2 件のコメント:
よくわからないのにコメント付けて申し訳ないですが，
for (i in 1:nrow(Tree)){
for (j in 1:nrow(IncW3)){
DistW3[i*j,1:4]<-c(Tree[i,4:5],IncW3[j,1:2])
：
のようなところは，添え字が i, j のときに，対象 i と対象 j の距離と角度というつもりなのではないでしょうか?
結果を格納する場所の添え字に i*j を使ってしまっては，正しい場所に正しい結果が格納されないのではないでしょうか?
つまり，例えば i=3, j=4 のときも，i=2, j=6 のときも i*j=12 なので，i*j が 12 になる場合に同じ場所にどんどん上書きされていくことになりますね。必然，i=3, j=4 のときの結果は上書きされてなくなってしまいます。
R では，for を使わないでベクトル化演算を用いることで，格段の実行速度が得られますよ。
コメントありがとうございます。
たしかにご指摘の通りです。このコード自体も後に書き換えて，i*jを使わないように修正しています。
また，Rとの比較という点も，ご指摘のように結果はフェアでとは言いがたいですね。今後修正します。
ありがとうございました。
コメントを投稿