Python
numpy
Fortran90
情報処理

情報処理の問題をfortranとpythonで書く!#03 - 3次元空間上の2点間の距離

はじめに

大学の授業でfortranを使った情報処理の授業があり、その課題をpythonとfortranの両方で書いていきたいと思います。

プログラミングについては専門で習っているわけではなく、初心者なので誤りなどあると思います。もし、よかったらコメントで指摘をお願いします。

以下にそれぞれの環境を記しておきます

fortran90

windows10
cygwin2.8.1

python

windows10
python3.6.2

課題内容 問5-15

問題 15
3 次元空間の N 個の点を考える.これらの位置を座標(x, y, z)の形で読み込み,2 点間の距離をすべて計算せよ.さらに,距離の小さい順に並べ替えて,座標とともに出力せよ.ただし,𝑁 ≤ 10としてよい。

fortran

出力結果

  0.100000001      0.200000003      0.300000012
   1.10000002       2.20000005       3.29999995
   4.40000010       5.50000000       6.59999990
   7.69999981       8.80000019       9.89999962
distance=    3.74   ( 0.10 0.20 0.30)   ( 1.10 2.20 3.30)
distance=    5.72   ( 1.10 2.20 3.30)   ( 4.40 5.50 6.60)
distance=    5.72   ( 4.40 5.50 6.60)   ( 7.70 8.80 9.90)
distance=    9.29   ( 0.10 0.20 0.30)   ( 4.40 5.50 6.60)
distance=   11.43   ( 1.10 2.20 3.30)   ( 7.70 8.80 9.90)
distance=   14.96   ( 0.10 0.20 0.30)   ( 7.70 8.80 9.90)

入力ファイル

coordinate.txt
0.1 0.2 0.3
1.1 2.2 3.3
4.4 5.5 6.6
7.7 8.8 9.9
toi5_15.f90
program toi5_15
implicit none
integer:: n,i,k1,k2,k3,k4,k5,j,buf2
real::a(3,10),l2(100,100),l1(100,100),buf1,fl(100,3)

    open(12,file="coordinate.txt",status="old",err=100)

    DO i=1,10
        read(12,*,end=200) a(1,i),a(2,i),a(3,i)

        write(*,*) a(1,i),a(2,i),a(3,i)
    end do

200 i=i-1
    j=1
    DO k1=1,i
        DO k2=k1+1,i
            l2(k1,k2)=(a(1,k1)-a(1,k2))**2+(a(2,k1)-a(2,k2))**2+(a(3,k1)-a(3,k2))**2
            l1(k1,k2)=sqrt(l2(k1,k2))
            fl(j,1)=l1(k1,k2)
            fl(j,2)=k1
            fl(j,3)=K2
            j=j+1
        end do
    end do

    DO k3=1,j-1
        DO k4=1,j-2
            if(fl(k4,1)>fl(k4+1,1)) then
                buf1=fl(k4,1)
                fl(k4,1)=fl(k4+1,1)
                fl(k4+1,1)=buf1
                buf2=fl(k4,2)
                fl(k4,2)=fl(k4+1,2)
                fl(k4+1,2)=buf2
                buf2=fl(k4,3)
                fl(k4,3)=fl(k4+1,3)
                fl(k4+1,3)=buf2
            end if
        end do
    end do

    DO k5=1,j-1
        write(*,'(A,f8.2)',advance="no") "distance=",fl(k5,1)
        write(*,'(A,f5.2,f5.2,f5.2,A)',advance="no")"   (",a(1,int(fl(k5,2))),a(2,int(fl(k5,2))),a(3,int(fl(k5,2))),")"
        write(*,'(A,f5.2,f5.2,f5.2,A)')"   (",a(1,int(fl(k5,3))),a(2,int(fl(k5,3))),a(3,int(fl(k5,3))),")"
    end do

stop
100 write(*,'(A)') "no file!"
    stop

end program toi5_15

かなり助長になってしまった気がしますが,とりあえずできました。ソートをもとの2次元配列のままできないかなと考えましたが、あきらめて配列をもう一度作ってバブルソートで並び替えています。

今回学んだこと
fortran90には1文39(132)文字の制限があり,それをこえるとtruncated lineというエラーが出る。→@septcolor さんの指摘で132文字が正しいようです
→自分の場合最後のwrite文を1行で書くとエラーになりました

@BlackDice さんがわざわざコードを書いてくださり、多次元配列のソートやファイル指定、変数宣言など様々なことを学ばせてもらいました。自分のコードを編集しなおそうとしたのですが、知識不足で@BlackDice さんのコードで学んだことを書くとコピペとほとんど変わらなくなってしまうため、以下に気を付けることをリストアップすることにしました。

  1. 行番号はfortranのみの機能であるため、ファイルの終端処理についてはiostatを用いる(iostat指定子はファイル読み取り終了の場合は負、通常の完了は0、エラー状態が発生すると正の整数値を返すようです) https://jp.xlsoft.com/documents/intel/cvf/vf-html/pg/pg21_02_02.htm
  2. DO文などのmaxを指定するときはパラメータを宣言してその変数に数字を代入することで、あとで変更するときの手間を省く
  3. 変数は定数として使用するときはparameterを宣言し、::の後に初期値を代入する(parameterのときは初期値を指定しないといけないみたいです。::は初期値を指定するときのみ必要であるとのこと、以下参考サイト) http://www.nag-j.co.jp/fortran/FI_4.html
  4. fortranのsubroutineは配列を受け取るときに次元を落として受け取ることができるため、2次元配列を1次元配列に次元を落としてソートすることが可能(こんなことができるとは知りませんでした。よく配列のソートを使うので、参考にさせてもらいます。)
  5. 書式指定は、同様の指定をする場合、f5.2,f5.2,f5.23f5.2などとまとめることができる(指定子の文字列がシングルクォートで区切られている場合は文字列はダブルクォートで区切る必要があるみたいです)http://www-space.eps.s.u-tokyo.ac.jp/~amano/education/fortran/chap06.html
  6. stop関数はstop "文字列"のように書くと文字列を出力してプログラムが終了する
  7. 最後に気を付けることはサブルーチンでは変数宣言でintentをつけないとエラーで返ってくることです(いままで、サブルーチンは使っていませんでした)

それと、fortranには構造体はないものと思っていたのですが、fortran90から構造体が実装されているみたいですね、今度配列の代わりに使ってみようと思います

python

だいぶ遅れましたが、pythonのコードを書いていきます。いままでの投稿でいくつかコメントをいただき、勉強してnumpyなどを使っていますが、まだ習得できているとまではいかず、ところどころlistなどを混ぜてしまいました:frowning2:

とりあえず、出力結果から

[[ 0.1  0.2  0.3]
 [ 1.1  2.2  3.3]
 [ 4.4  5.5  6.6]
 [ 7.7  8.8  9.9]]
[[3.7416573867739413, array([ 0.1,  0.2,  0.3]), array([ 1.1,  2.2,  3.3])],
 [5.715767664977295, array([ 4.4,  5.5,  6.6]), array([ 7.7,  8.8,  9.9])],
 [5.715767664977295, array([ 1.1,  2.2,  3.3]), array([ 4.4,  5.5,  6.6])],
 [9.2881645118936174, array([ 0.1,  0.2,  0.3]), array([ 4.4,  5.5,  6.6])],
 [11.43153532995459, array([ 1.1,  2.2,  3.3]), array([ 7.7,  8.8,  9.9])],
 [14.962620091414472, array([ 0.1,  0.2,  0.3]), array([ 7.7,  8.8,  9.9])]]

見た目が汚いのは、書式指定をできていないためです...

プログラム

nda5-15.py
import numpy as np
from operator import itemgetter
from pprint import pprint

p=np.loadtxt("coordinate.txt")
print(p)
lon=[]
for i in range(0,len(p)):
    for k in range(i+1,len(p)):
        h=p[i]-p[k]
        m=np.linalg.norm(h)
#        print(m)
        inte=[]
        inte.insert(0,m)
        inte.append(p[i])
        inte.append(p[k])
        lon.insert(1,inte)


pprint(sorted(lon,key=itemgetter(0)))

とりあえず、numpyを使うことでさまざまな関数を利用できてかなり短く書くことができたのですが、numpyの配列が追加などができないところで戸惑ってしまいました。
今回のように、for文で1回の処理で一つの要素が出てきて、それをnumpyの配列に入れる時はどうしたらいいのか悩んで...結局追加できるpythonのlistに入れ込んで、そのあとにソートしました。その結果,要素が複雑で書式指定も断念しています。

要素ごとに段落をpprintで変更できたのはよかったところです。(@wakamezake さんにコメントで教えていただきました)