多角形の内外判定

多角形の内外判定をする必要があって、pythonで作ってみた。

参考にしたのがmatlabのinpolygon関数（http://www.mathworks.co.jp/jp/help/matlab/ref/inpolygon.html)

だけど、諸事情でend pointが閉じていないものにした。

使ったのはCrossing Number Algorithmという一番単純なもので、多角形の辺が、ある点より上にある数が奇数であれば、その点は多角形の内側にあるというもの。

import matplotlib.pyplot as plt
import numpy as np

def main():
    n = 2000
    x = np.ndarray(n, dtype=float)
    y = np.ndarray(n, dtype=float)
    inside = np.ndarray(n, dtype=bool)
        
    x = (np.random.rand(n) - 0.5)*2.5
    y = (np.random.rand(n) - 0.5)*2.5
    inside[:] = False
        
    x1, y1 = [], []
    
    # define a polygon
    for i in range(5):
        x1.append(np.cos(i*2.*np.pi/5.))
        y1.append(np.sin(i*2.*np.pi/5.))
    
    # check the points inside the polygon
    for i, (sx, sy) in enumerate(zip(x, y)):
        inside[i] = inpolygon(sx, sy, x1, y1)
    
    # plot results
    x1.extend([x1[0]])
    y1.extend([y1[0]])
    plt.plot(x[inside], y[inside], 'ro')
    plt.plot(x[inside==False], y[inside==False], 'ko')
    plt.plot(x1, y1)
    plt.savefig("inpoly.png")
        
def inpolygon(sx, sy, x, y):
    ''' 
    x[:], y[:]: polygon
    sx, sy: point
    '''     
    np = len(x)
    inside = False
    for i1 in range(np): 
        i2 = (i1+1)%np
        if min(x[i1], x[i2]) < sx < max(x[i1], x[i2]):
            #a = (y[i2]-y[i1])/(x[i2]-x[i1])
            #b = y[i1] - a*x[i1]
            #dy = a*sx+b - sy
            #if dy >= 0:
            if (y[i1] + (y[i2]-y[i1])/(x[i2]-x[i1])*(sx-x[i1]) - sy) > 0:
                inside = not inside

    return inside

if __name__ == "__main__":
    main()

ためしてみると、うまくいった。こんな単純な計算でできるんだなあ。

pythonで構造体的なものを

しばらく何も書いていなかったけれど、久々に投稿。
最近はpythonでコードをよく書いている。だけどまだまだ勉強し始めたばかりで、使いこなせていないのが現状だけど、だんだん使い慣れて気持ちよくコードが書けるようになってきた。

pythonの良いところはいろいろある。例えば、気に入っているのが、

• perlと比べるときれいにコードを書ける（だけどperlは大好き）
• numpy, scipy：数値計算が簡単
• ipython：変数の中身を簡単に見れたりとか
• ctypes：cとかfortranを自由に呼び出せる。okadaプログラムDC3D.f（http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_J.html）も簡単にwrapperを作れる。
• generatorとか、変な機能があってわくわくする。

きりがないけれど、そもそもpythonを使い始めた理由が数値計算にも結構使えそうだったため。なので、c、fortranで良く使う構造体が使えないかな、と思い始めた。classで実装すればいいけれど、それでは遅そうだし（未確認）。

ということでnumpy.dtypeがそれっぽそう？なので試してみる。とりあえず国土地理院（terras.gsi.go.jp）から、GEONETの座標値を落としてきて、plotさせてみた。


なんと、data['date']、data['x']のようにしてリスト？が出せる。さらに、添字が連続しているようで、data[i][0]、data[i][1]のようにも使える。これはかなり良いかもしれない。速度はどれくらい違うんだろう。numpy、すごいなあ。

import numpy as np
import datetime as dt
import matplotlib.pyplot as plt

def main():
    posf = "950255.11.pos"

    # structure?
    pos_t = np.dtype([('date', dt.datetime),
                      ('x', np.float), 
                      ('y', np.float), 
                      ('z', np.float), 
                      ('lat', np.float), 
                      ('lon', np.float), 
                      ('h', np.float), 
                      ])

    with open(posf, 'rt') as f:
        for line in f:
            if line[1:6] == "EPOCH":
                nd = int(line[72:75]) # number of data
            if line[0:5] == "+DATA":
                break

        if nd > 0:
            data = np.ndarray([nd], dtype=pos_t)

        i = 0
        for line in f:
            if line[0:1] == "*":
                continue
            elif line[0:5] == "-DATA":
                break
            else:
                date_string    = line[1:20]
                x1,   y1,   z1 = line[21:38], line[39:56], line[57:74]
                lat1, lon1, h1 = line[75:92], line[93:110], line[111:128]

                #data[i][0] = dt.datetime.strptime(date_string, "%Y %m %d %H:%M:%S")
                #data[i][1], data[i][2], data[i][3] = np.float(x1), np.float(y1), np.float(z1)
                #data[i][4], data[i][5], data[i][6] = np.float(lat1), np.float(lon1), np.float(h1)
                data[i]['date'] = dt.datetime.strptime(date_string, "%Y %m %d %H:%M:%S")
                data[i]['x'], data[i]['y'], data[i]['z'] = np.float(x1), np.float(y1), np.float(z1)
                data[i]['lat'], data[i]['lon'], data[i]['h'] = np.float(lat1), np.float(lon1), np.float(h1)

                i+=1

    plt.plot_date(data['date'], data['x'])
    plt.savefig("test.png")
    plt.show()

    return data

if __name__ == '__main__':
    data = main()

ところで、プロットさせてみた座標値（x）がこちら。東日本大震災で大きく観測点が移動している。x座標値で10cmくらい。石川県なんだけどな、この観測点・・・。

fortran コマンドライン引数処理

今日はfortranでコマンドライン引数を処理するサブルーチンを作ってみたので貼っておく。

fortranは言語仕様ではコマンドライン引数へのアクセスができないけども、通常のcompilerではiargc関数、getargサブルーチンが拡張仕様として備わっているので、最近では問題ないのかなと思う。

* Fortran 2003に対応しているcompilerならば、下記の関数を使う。
count = command_argument_count(), CALL get_command_argument(i, arg)

インターフェースはperlの標準モジュール"Getopt::Std"を真似した。使い方も大体同じ。
http://perldoc.perl.org/Getopt/Std.html

例えば、"-f"とか、"-i"とかよくやるやり方でfortranに引数を渡せるというまぁまぁ便利なものになった。

$ ./a.out -f filename -i
 f:filename
 i:1
 s:0

ソースは下記のようなもの。例えば下記のように"f:is"というオプションリストを渡すと、fについては文字列が格納され、i、sにはフラグ（0:off,1:on）がoptsという文字列配列に返される。順番は"f:is"で指定した順に入る。

program getarg
    implicit none

    character(len=256) :: opts(3)

    call getopt("f:is", opts)

    print *, "f:", trim(opts(1))
    print *, "i:", trim(opts(2))
    print *, "s:", trim(opts(3))

    stop
end program getarg

! --------------------------------------------------
subroutine getopt(optlist, opts)
! --------------------------------------------------
! parse command line arguments specified with
! single-character swiches.
!
! Usage: 
! call getopt("f:is", opts)
!
! where, switches followed by ":" takes arguments, 
! and others work as flag.
!
! for example, "./a.out -f file -i" returns
!  opts(1): file
!  opts(2): 1
!  opts(3): 0
! --------------------------------------------------

    type opt
        character(len=1)   :: name
        logical            :: arg  ! flag opt takes arg or not.
        character(len=256) :: val
    end type opt

    character(len=*), intent(in)  :: optlist
    character(len=*), intent(out) :: opts(:)

    integer                :: i, j, nopts
    character(len=256)     :: argv
    character(len=1)       :: optname
    type(opt), allocatable :: opts_tmp(:)

    allocate(opts_tmp( size(opts) ))

    ! parse optlist
    nopts = 0
    do i =1,len_trim(optlist)
        optname = optlist(i:i)
        if (optname /= ":") then
            nopts = nopts+1
            if (nopts > size(opts)) then
                print *, "error: number of options exceeds the size of opts"
            endif

            opts_tmp(nopts)%name = optname
            opts_tmp(nopts)%arg  = .false.
            opts_tmp(nopts)%val  = "0"
        endif

        if (optname == ":") opts_tmp(nopts)%arg = .true.
    enddo

    ! check command line arguments
    do i = 1, iargc()
        call getarg(i, argv)

        if (argv(1:1) == "-") then
            ! serach option list
            do j = 1,nopts
                if (opts_tmp(j)%name == argv(2:2)) then
                    ! found option, set value
                    if (opts_tmp(j)%arg .eqv. .true.) then
                        call getarg(i+1, opts_tmp(j)%val)
                    else
                        opts_tmp(j)%val = "1"
                    endif
                endif
            enddo
        endif
    enddo

    do i = 1, nopts
        opts(i) = opts_tmp(i)%val
    enddo

    deallocate(opts_tmp)

    return
end subroutine getopt

fortran用行列関係の関数いくつかまとめ

fortranで行列関係の関数をいくつか作ったので、まとめておく。といっても、MATLABの良く使いそうなやつをいくつか真似して作ってみただけ。

fortranは数値計算は何も考えなくても早いけども、ちょっとしたことが不便すぎる。いろいろツールを作っておかないとそのちょっとしたことで時間を取られてしまうし、コードも汚くなってしまう。まだfortran77のソースも良く見るけど、fortran90, fortran95の機能を知りながら使う気には到底なれないなあ。

それにしてもMATLABで便利なのって、結局、転置行列X'とか行列積X*Yが表現できることだったり、変数がそのまま残ることだったりするような気がする（速度を出すにはテクニックが必要で、まだ全然使いこなせない）。

あとは、MATLABのsave関数的なものがあると便利そうではある。

! ---------------------------------------------------------
module matrix
! ---------------------------------------------------------
! inv_lapack
!   Ai = inv_lapack(A) : return inverse of A(m,n) using lapack.
!
! whos
!   call whos(X) : shows size, bytes, and type of X.
!                  only real(8) and integer are supported.
!   
! linspace
!   X = linspace(x0,x1,n) : return linearly spaced vector just same as MATLAB.
!
! eye, ones, zeros (same as MATLAB, Octave)
!   X = eye(n,[m])   : get identity matrix
!   X = ones(n,[m])  : get matrix of all ones
!   X = zeros(n,[m]) : get matrix of all zeros
!
! show_dmat
!   call show_dmat(X,[fmt]) : show real(8) X(n,m). 
!
! ---------------------------------------------------------
! Note:
!   build: gfortran matrix.f90 main.f90 -llapack
! ---------------------------------------------------------
! Created by sako  2013/05/14 
! ---------------------------------------------------------
    implicit none

    interface whos
        module procedure whos_d1, whos_d2, whos_i1, whos_i2
    end interface

contains
    ! --------------------------------
    ! estimate inverse of array A(m,n).
    ! Ai = inv_lapack(A)
    ! --------------------------------
    function inv_lapack(A) result (AI)
        real(8), intent(in) :: A(:,:)
        real(8) :: AI(size(A,1),size(A,2))
        integer, parameter :: NB = 64
        integer :: info
        integer :: m, n

        ! for lapack
        integer :: LWORK
        integer, allocatable :: IPIV(:) ! dimension N
        real(8), allocatable :: WORK(:) ! dimension LWORK

        ! get size of A(m,n)
        m = size(A,1)
        n = size(A,2)

        LWORK = n*NB
        allocate(IPIV(n))
        allocate(WORK(LWORK))

        AI = A
        call DGETRF(m, n, AI, m, IPIV(1:min(m,n)), info)  ! factorize
        call DGETRI(n, AI, m, IPIV, WORK, LWORK, info)  ! inverse

        deallocate(IPIV)
        deallocate(WORK)

        if (INFO==0) then
            !write(*,*) "succeeded."
        endif
        if (INFO/=0) then
            print *, 'error: failed to estimate inverse at "inv_lapack".'
        endif

        return
    end function

    ! --------------------------------
    ! show variable size, bytes, type
    ! like whos of MATLAB. (call whos(X))
    ! Size                     Bytes    Kind
    ! ====                     =====    ====
    ! 1x100                      800  double
    ! --------------------------------
    subroutine whos_d1(x)
        real(8), intent(in) :: x(:)
        call show_whos(1, size(x), kind(x), 'double')
        return
    end subroutine whos_d1

    subroutine whos_d2(x)
        real(8), intent(in) :: x(:,:)
        call show_whos(size(x,1), size(x,2), kind(x), 'double')
        return
    end subroutine whos_d2

    subroutine whos_i1(x)
        integer, intent(in) :: x(:)
        call show_whos(1, size(x), kind(x), 'integer')
        return
    end subroutine whos_i1

    subroutine whos_i2(x)
        integer, intent(in) :: x(:,:)
        call show_whos(size(x,1), size(x,2), kind(x), 'integer')
        return
    end subroutine whos_i2

    subroutine show_whos(n, m, b, valtype)
        integer, intent(in) :: n, m, b
        character(len=*), intent(in) :: valtype

        character(len=11) :: strn, strm
        character(len=16) :: strb

        write(strn,'(i11)') n
        write(strm,'(i11)') m
        write(strb,'(i16)') b*n*m  ! bytes

        write(*,'(a48)') "          Size                     Bytes    Kind"
        write(*,'(a48)') "          ====                     =====    ===="
        write(*,'(a11,"x",a11,1x,a16,1x,a7)') adjustr(strn), adjustl(strm), adjustr(strb), adjustr(valtype)
        return
    end subroutine

    ! --------------------------------
    ! return linearly spaced matrix. 
    ! just same as linspace of MATLAB.
    ! --------------------------------
    function linspace(x0,x1,nd) result(X)
        real(8), intent(in) :: x0, x1
        integer, intent(in), optional :: nd
        real(8), allocatable :: X(:)

        real(8) :: dx
        integer :: n, i, ierr
        
        if (present(nd)) then
            n = nd
        else
            n = 100
        endif

        ! error check
        if (n < 1) n = 1  ! force to be 1

        allocate(X(n), stat=ierr)
        if (ierr /= 0) then
            print *, 'error: failed to allocate memory at "linspace".'
            stop
        endif

        ! special case for n==1
        if (n == 1) then
            X(1) = x1
            return
        endif

        dx = (x1-x0)/dble(n-1)
        do i = 1,n
            X(i) = x0 + dx*dble(i-1)
        enddo

        return
    end function linspace

    ! --------------------------------
    ! return idendity matrix E(n[,m])
    ! --------------------------------
    function eye(n,m) result(E)
        !implicit none
        integer, intent(in) :: n
        integer, intent(in), optional :: m
        real(8), allocatable :: E(:,:)
        integer :: i, n2, ierr

        if (present(m)) then
            n2 = m
        else
            n2 = n
        endif

        allocate(E(n,n2), stat=ierr)
        if (ierr /= 0) then
            print *, 'error: failed to allocate memory at "matrix_eye".'
            stop
        endif


        ! initialize
        E(1:n,1:n2) = 0.d0
        do i = 1,n
            E(i,i) = 1.d0
        enddo
        return
    end function eye

    ! --------------------------------
    ! return matrix X(n[,m]) of all ones
    ! --------------------------------
    function ones(n,m) result(X)
        !implicit none
        integer, intent(in) :: n
        integer, intent(in), optional :: m
        real(8), allocatable :: X(:,:)
        integer :: i, n2, ierr

        if (present(m)) then
            n2 = m
        else
            n2 = n
        endif

        allocate(X(n,n2), stat=ierr)
        if (ierr /= 0) then
            print *, 'error: failed to allocate memory at "matrix_ones".'
            stop
        endif

        ! initialize
        X(1:n,1:n2) = 1.d0
        return
    end function ones

    ! --------------------------------
    ! return matrix X(n[,m]) of all zeros
    ! --------------------------------
    function zeros(n,m) result(X)
        !implicit none
        integer, intent(in) :: n
        integer, intent(in), optional :: m
        real(8), allocatable :: X(:,:)
        integer :: i, n2, ierr

        if (present(m)) then
            n2 = m
        else
            n2 = n
        endif

        allocate(X(n,n2), stat=ierr)
        if (ierr /= 0) then
            print *, 'error: failed to allocate memory at "matrix_zeros".'
            stop
        endif

        ! initialize
        X(1:n,1:n2) = 0.d0
        return
    end function zeros

    ! --------------------------------
    ! show matrix A(m,n)
    ! --------------------------------
    subroutine show_dmat(A,fmtf)
        !implicit none
        real(8), intent(in) :: A(:,:)
        character(len=*), intent(in), optional :: fmtf
        character(len=25) :: fmtc
        integer :: m, n, i, j
        m = size(A,1); n = size(A,2)

        if (present(fmtf)) then
            fmtc = fmtf
        else
            fmtc = '(f13.4," ")'
        endif
        do i=1,m
            do j=1,n
                write(*,fmtc, advance='no') A(i,j)
            enddo
            write(*,*) ""
        enddo
        return
    end subroutine show_dmat

end module matrix

fortranの副プログラムからの返り値は参照渡しかどうか

ふと、fortranのfunctionとsubroutineではどのように返り値が返ってきているのか気になった。fortranの引数は参照渡しなので、恐らくsubroutineに関しては配列をallocateして返したりしても同じアドレスに違いない。

だけど、functionに関してはどうなのだろう？というのも、function内でallocateして配列を返すようなものを結構書いていたから。fortran95からはスコープから外れるとメモリを解放してくれることになっているので、deallocateなどしていない。

ということで、試してみた。

program test_allocate
    implicit none

    integer, parameter :: n = 1000000
    integer :: array1(n)
    integer, allocatable :: array2(:)
    array1(1:n) = 1

    ! function
    array1 = allocate_array(n)
    print *, loc(array1)

    ! subroutine
    call allocate_array_sub(n, array2)
    print *, loc(array2)
    stop
contains

subroutine allocate_array_sub(n, array)
    integer, intent(in) :: n
    integer, allocatable, intent(out) :: array(:)

    allocate(array(n))
    array(1:n) = 1

    print *, loc(array)

    return
end subroutine allocate_array_sub

function allocate_array(n) result(array)
    integer, intent(in) :: n
    integer, allocatable :: array(:)

    allocate(array(n))
    array(1:n) = 1

    print *, loc(array)

    return
end function allocate_array
end program test_allocate

すると、こんな結果に。

      140120599949328  ! function内でallocate
              6299808  ! functionの外
             39948288  ! subroutine内でallocate
             39948288  ! subroutineの外

ということで、試したのはgfortranだけれども、function内でallocateした配列は値のコピーが返り、subroutineから返る配列は参照渡しで返るみたいだ。本当に言語仕様なのかはよく分からない。

fortranのファイルオープン関数

fortranのファイルオープンをしようとすると、いつも引数の指定方法を忘れてしまう。cとかperlみたいにopenするのに慣れすぎた。

いちいち外部入力装置番号を覚えておかなくてはいけないとか、オプションも自由気ままなんて、あまりにカオスすぎる。

ということで、perlなどのようにファイルを開くことができるようなwrapperを作ってみた。これ以外で使う時は、scratchでファイルを作ってOut of Coreのプログラムを書く時くらいだろうか。バイナリなんかめったにいじらないし。

こんな感じで使える。

    funit = open_file('>', "test.dat")   !書き込み専用
    funit = open_file('<', "test.dat")   !読み込み専用
    funit = open_file('>>', "test.dat")  !追加書き込み専用
    funit = open_file('+<', "test.dat")  !読み書き可能
    funit = open_file('+>', "test.dat")  !読み書き可能（ファイル新規作成）
    funit = open_file('+>>', "test.dat") !読み書き可能追加書き込み

本体のモジュールが次のような感じ。未使用の外部入力装置番号を検索してそれを自動で割り当てて、あとはただのwrapper。未使用の外部入力装置番号のインデックスは、モジュールのprivateメンバ変数にしてしまっても良いかもしれない。

module file
    implicit none

    ! define error codes
    integer, parameter :: ERR_FOPEN = -1
    integer, parameter :: ERR_FMODE = -2
    integer, parameter :: ERR_FUNIT = -8

contains
    ! ----------------------------------------
    ! open file like 'c' and 'perl'.
    ! if success, returns unit number.
    ! if error, returns minus number.
    ! 
    ! write only:
    !   funit = open_file('>', file_name)
    !   funit = open_file('w', file_name)
    ! ----------------------------------------
    integer function open_file(mode, file_name)
        character(len=*), intent(in) :: mode
        character(len=*), intent(in) :: file_name

        integer :: ios
        character(len=7) :: fstatus
        character(len=9) :: faction
        character(len=6) :: fpos
        integer :: funit
        funit = get_file_unit()

        if ((mode == '<') .or. (mode == 'r')) then
            ! read only
            fstatus='old'
            faction='read'
            fpos='rewind'
        elseif ((mode == '>') .or. (mode == 'w')) then
            ! write only
            fstatus='replace'
            faction='write'
            fpos='rewind'
        elseif ((mode == '>>') .or. (mode == 'a')) then
            ! append only
            fstatus='unknown'
            faction='write'
            fpos='append'
        elseif ((mode == '+<') .or. (mode == 'r+')) then
            fstatus='old'
            faction='readwrite'
            fpos='rewind'
        elseif ((mode == '+>') .or. (mode == 'w+')) then
            fstatus='replace'
            faction='readwrite'
            fpos='rewind'
        elseif ((mode == '+>>') .or. (mode == 'a+')) then
            fstatus='unknown'
            faction='readwrite'
            fpos='append'
        else
            ! invalid mode
            open_file = ERR_FMODE
            return
        endif

        ! file open
        open(unit=funit, file=file_name, iostat=ios, status=fstatus, &
             action=faction, position=fpos)

        ! check file open error
        if (ios==0) then
            open_file = funit
        else
            open_file = ERR_FOPEN
        endif
        return
    end function open_file

    ! ----------------------------------------
    ! get file unit which has not opened.
    ! return minus if error found.
    ! ----------------------------------------
    integer function get_file_unit()
        integer :: i
        integer, parameter :: MAX_UNIT = 999
        logical :: fopened
        
        do i = 100, MAX_UNIT
            inquire(unit=i, opened=fopened)
            if (fopened .neqv. .true.) then
                get_file_unit = i
                return
            endif
        enddo

        ! not found unopened file unit.
        get_file_unit = ERR_FUNIT
        return    
    end function get_file_unit
end module file

RTKLIB v2.4.2が公開された

RTKLIB v2.4.2が公開されていた。
毎日HPを確認していたんだが、ついに！

http://www.rtklib.com/rtklib.htm

ベータ版はほとんど見てないけど、v2.4.1からどれくらい変わったのだろうか？
PPP-ARの仕組みは全く勉強していないので、これを機にソースを読んで勉強したい。

Okada(1992) DC3Dの使い方

Okada (1992) [Bull. Seism. Soc. Am., 82, 1018-1040]の式のfortran77プログラムが下記で公開されている。
http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_J.html

ある場所（緯度経度）での変位を出力しようと思ったら、結構頭がこんがらがった。というわけで、wrapperを作ったので、プログラムをメモ代わりにアップしておく。断層パラメータの記載は一般的な定義のままで、位置は断層上端左端にした。ああ、ややこしかった。

注意としては、深さは正の値を入れること、長さLの場合はAL1=0, AL2=Lとすること、幅Wの場合は、AW1=-W, AW2=0とすること。
一般的な断層パラメータの定義が数式と少し違い（走向は北から反時計回りとか）、注意が必要ですね。すぐ忘れそうなので、やはりwrapperを作っておいた方が良いな。もう少しモジュールがたまってきたら、arでライブラリにしてしまいたい。

ちなみに、こんな感じで呼び出す。

test.f90:
program test
    use okada92_util
    implicit none

    type(disp) :: du
    type(flt)  :: f
    real(8) :: lat, lon, h

    integer :: i, j

    f%lambda = 1.d0   ! ずるしている
    f%mu     = 1.d0
    f%lat    = 32.7725d0
    f%lon    = 136.0282d0
    f%depth  = 5000.d0
    f%str    = 261.766d0
    f%dip    = 8.903d0
    f%rake   = 90.d0
    !f%rake   = 0.d0
    f%length = 100000.d0
    f%width  = 25263.533d0
    f%slip   = 5.d0

    do i = -50, 50
        do j = -50, 50
            lat = f%lat + dble(i)*0.1d0
            lon = f%lon + dble(j)*0.1d0
            h   = 0.d0

            du = flt_disp(f, lat, lon, h)
            write(*,'(F7.3,",",F7.3,3(",",F15.7))') lon, lat, du%e, du%n, du%u
        enddo
    enddo

    stop
end program test

モジュールは下記のようなもの。

okada92_util.f90:
module okada92_util
    type disp
        real(8) :: e, n, u
    end type disp

    type flt
        !real(8) :: alpha          ! (lambda+mu)/(lambda+2*mu) !αをインプットにした方が良かったかも
        real(8) :: lambda, mu    ! (lambda+mu)/(lambda+2*mu)
        real(8) :: lat, lon      ! latitude, longitude (deg) at the corner of upper edge
        real(8) :: depth         ! fault depth (positive is ...)
        real(8) :: str           ! strike angle (deg)
        real(8) :: dip           ! dip angle    (deg)
        real(8) :: rake          ! rake angle   (deg)
        real(8) :: length, width ! fault length, width (m)
        real(8) :: slip          ! slip (m)
    end type flt

    ! ellipsoid
    type ellip
        real(8) :: a, b, e2, f, fi
    end type ellip

    ! parameters
    real(8), parameter :: PI = 3.1415926535897932d0
    real(8), parameter :: D2R = PI/180.d0
    real(8), parameter :: R2D = 180.d0/PI
 
contains
    ! ----------------------------------------------------------
    type(disp) function flt_disp(f, lat, lon, h)
    ! ----------------------------------------------------------
    ! Calculates displacement at (lat, lon, h) due to a 
    ! rectangular fault. The equations are given by 
    ! Okada (1992) [Bull. Seism. Soc. Am., 82, 1018-1040].
    !
    ! This function is a wrapper of 'DC3D' on the website:
    ! http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_E.html
    !
    ! Input:
    !   f : type(flt), includes
    !    f%lambda, mu    ! (lambda+mu)/(lambda+2*mu)
    !    f%lat, lon      ! latitude, longitude (deg) at the corner of upper edge.
    !    f%depth         ! fault depth (positive is downward) 
    !    f%str           ! strike angle (deg)
    !    f%dip           ! dip angle    (deg)
    !    f%rake          ! rake angle   (deg)
    !    f%length, width ! fault length, width (m)
    !    f%slip          ! slip (m)
    !   lat, lon, h : latitude(deg), longitude(deg), h(m)
    !                 positive h is upward.
    ! Output:
    !   disp        : type(disp),
    !     disp%e         ! eastward disp. (m)
    !     disp%n         ! northward disp. (m)
    !     disp%u         ! vertical disp. (m)
    !
    ! Created:
    !   2013/04/29  Sako
    ! ----------------------------------------------------------
        type(flt), intent(in) :: f
        real(8),   intent(in) :: lat, lon, h

        ! inputs of DC3D
        real(8) :: a, x, y, z, d, dip
        real(8) :: al1, al2, aw1, aw2
        real(8) :: disl1, disl2, disl3

        ! output of DC3D
        real(8) :: ux, uy, uz
        real(8) :: uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz
        integer :: iret

        ! tmp
        real(8) :: dx, dy, t

        !real(8) :: alpha          ! (lambda+mu)/(lambda+2*mu)
        a     = (f%lambda+f%mu)/(f%lambda+2.d0*f%mu)
        d     = f%depth   ! 断層上端の深さ。深さは正の値にする。
        dip   = f%dip
        al1   = 0.d0
        al2   = f%length  ! 断層の長さ。
        aw1   = -f%width  ! 断層の幅。図を良く見るとAW1はマイナスということで結構紛らわしい。
        aw2   = 0.d0
        disl1 = f%slip * dcos(f%rake*D2R)  ! すべり量を各軸に分解。
        disl2 = f%slip * dsin(f%rake*D2R)
        disl3 = 0.d0

        ! 目的の緯度経度をローカル座標系のx, yに変換する。走向の定義が北からなので、ここもややこしい。
        call bl2xy(f%lat, f%lon, lat, lon, dx, dy)
        t = (-f%str+90.d0)*D2R;     ! rotate to local coordinate
        x = dx*dcos(t) + dy*dsin(t) ! t = -str+90
        y =-dx*dsin(t) + dy*dcos(t)
        z =-h
        
        !  SUBROUTINE  DC3D(ALPHA,X,Y,Z,DEPTH,DIP,
        !*                AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
        !*                UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET) 
        call DC3D(a, x, y, z, d, dip, al1, al2, aw1, aw2, disl1, disl2, disl3, &
                  ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz, iret)

        if (iret == 0) then
            flt_disp%e = ux*dsin(f%str*D2R) - uy*dcos(f%str*D2R)  ! 忘れずに座標変換。走向の向きは北からだからこうなる。
            flt_disp%n = ux*dcos(f%str*D2R) + uy*dsin(f%str*D2R)
            flt_disp%u = uz
        else
            ! singular location
            flt_disp%e = 0.d0
            flt_disp%n = 0.d0
            flt_disp%u = 0.d0
        endif

        return
    end function flt_disp

    ! ----------------------------------------------------------
    subroutine bl2xy(lat0, lon0, lat, lon, dx, dy)
    ! ----------------------------------------------------------
    ! Calculates distance between reference station(lat0, lon0) and 
    ! station2(lat, lon). 
    ! 
    ! Reference: Ohno S.(1987), Geodetic Methods (in Japanese).
    !   see equations [1.59] and [1.62]. 
    !
    ! Input:
    !   lat0, lon0 : position of station1 (reference)
    !   lat,  lon  : position of station2
    !
    ! Output:
    !   dx, dy     : distance between stations (sta1 -> sta2).
    ! 
    ! Created:
    !   2013/04/29 Sako
    ! ----------------------------------------------------------
        implicit none

        real(8), intent(in)  :: lat0, lon0, lat, lon
        real(8), intent(out) :: dx, dy

        type(ellip) :: ell
        real(8) :: dlat, dlon, Bm
        real(8) :: a, fi, f, b, e2


        ! WGS84
        ell = ellipsoid('WGS84')
        a = ell%a
        fi= ell%fi
        f = ell%f
        b = ell%b
        e2= ell%e2

        dlat = (lat - lat0)*D2R
        dlon = (lon - lon0)*D2R
        Bm = dsin(dlat)
        dx = a*dcos(lat0*D2R) / dsqrt(1.d0-e2*dsin(lat0*D2R)**2.d0)*dlon
        dy = a*(1.d0-e2)*dlat*(1.d0+3.d0/4.d0*e2-3.d0/4.d0*e2*dcos(2.d0*Bm)+1.d0/8.d0*e2*(dlat)**2.d0*dcos(2.d0*Bm))

        return
    end subroutine bl2xy

    ! ----------------------------------------------------------
    type(ellip) function ellipsoid(ellip_name)
    ! ----------------------------------------------------------
        implicit none
        character(len=*), intent(in), optional :: ellip_name
        real(8) :: a, fi, f
        
        ! default is GRS80
        a = 6378137.d0
        fi= 298.257222101d0

        if (present(ellip_name)) then
            if (ellip_name == 'GRS80') then
                a = 6378137.d0
                fi= 298.257222101d0
            elseif (ellip_name == 'WGS84') then
                a = 6378137.d0;
                fi= 298.257223563d0;
            else 
                ! ...
            endif
        endif

        f = 1.d0/fi

        ellipsoid%a = a
        ellipsoid%fi= fi
        ellipsoid%f = f
        ellipsoid%b = a*(1.d0-f)
        ellipsoid%e2= 1.d0 - (1.d0-f)**2.d0

        return
    end function ellipsoid
end module okada92_util

fortran90の構造体へのアクセス速度 メモ

前のブログ(http://sak12.blogspot.com/2013/04/matlab.html)で、MATLABの構造体の扱いには注意しなければいけないということをメモしたけれど、これまでfortranで使っていたものは大丈夫なのか不安になった。

 で、同じように２乗を計算してみた（1000万回）。
どうやら、fortranでは構造体へのアクセス速度は、特に何も考えなくても大丈夫そう。
構造体をよく使っていたのでドキッとしたのだけれど、一安心。

あと、MATLABが速いことを改めて実感。MATLABでドット演算子を使った場合、fortranのforallと速度はほぼ変わらない。
MATLABは遅いと良く言われたりしてるけど、うまく使えばかなり高速になりそうだし、きちんとMATLABの勉強しようかな。

program main
    implicit none

    type struct1
        real(8), allocatable :: x(:)
    end type struct1

    type struct2
        real(8) :: x
    end type struct2

    integer, parameter :: n = 10000000
    type(struct1) :: s1
    type(struct2) :: s2(n)
    real(8)       :: x(n), y(n)
    integer       :: stime(8), etime(8)
    integer :: i

    allocate(s1%x(1:n))

    forall (i=1:n)
        x(i)    = 2.d0
        s1%x(i) = 2.d0
        s2(i)%x = 2.d0
    end forall

    ! x ただの配列同士
    forall (i=1:n) y(i) = 0.d0
    call date_and_time(values=stime)
    forall (i=1:n)
        y(i) = x(i) * x(i)
    end forall
    call date_and_time(values=etime)
    write(*,'("elapsed time (x):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0

    ! s1　構造体が配列を持つ
    forall (i=1:n) y(i) = 0.d0
    call date_and_time(values=stime)
    forall (i=1:n)
        y(i) = s1%x(i) * s1%x(i)
    end forall
    call date_and_time(values=etime)
    write(*,'("elapsed time (s1):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0

    ! s2　構造体配列
    forall (i=1:n) y(i) = 0.d0
    call date_and_time(values=stime)
    forall (i=1:n)
        y(i) = s2(i)%x * s2(i)%x
    end forall
    call date_and_time(values=etime)
    write(*,'("elapsed time (s2):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0

    stop
end program main

結果はこのとおり。

elapsed time (x):    0.0250(s)
elapsed time (s1):    0.0260(s)
elapsed time (s2):    0.0260(s)

----- おまけ -----
perlでやってみると、こんなに時間がかかる・・・。まあ数値計算はやっちゃだめということですかね。そりゃそうか。

use strict;
use Time::HiRes qw( gettimeofday tv_interval );

my $N = 10_000_000;
my ($i, @x, $y);
foreach $i (1 .. $N) {
    push(@x, 2);
}

# estimate square of @x.
my $t0 = [gettimeofday];
for ($i=0;$i<$N;$i++) {
    $y = $x[$i]*$x[$i];
}
my ($seconds, $microseconds) = gettimeofday;
my $elapsed = tv_interval($t0, [$seconds, $microseconds]);
print "elapsed time: $elapsed\n";

$ perl test.pl 
elapsed time: 1.291174

MATLABのループ メモ

今日は久々にMATLABでコードを書こうと思いたった。
fortran77でよく使われているプログラムがあるのだけど、それをもっと楽に使用するため。

しかし、適当にMATLABでコードを書くと、変なloopを書いてしまいそうで怖い。
で、行列の掛け算のやり方と速度を確認したのでメモしておく。

単純なfor loopと、MATLAB特有のドット演算子を使ったものと、さらに添字にしたらどうなるかな、というもの。

結果は、当然だけどドット演算子を使わないとめちゃくちゃ遅い。あと、添字をつけたらfor loopより遅くなってしまった。使い方がまずいんだろうか。
あと、構造体を絡めて計算する場合は、構造体配列にすると遅くなるので、構造体に配列を入れてしまうのが良いみたい。

for loop: y(i) = x(i) * x(i);
Elapsed time is 0.015437 seconds.

vectorized loop: y = x .* x
Elapsed time is 0.002330 seconds.

vectorized loop2: y(1:n) = x(1:n) .* x(1:n)
Elapsed time is 0.021473 seconds.

using struct: y = s.x .* x
Elapsed time is 0.003305 seconds.

構造体のフィールドに配列を入れた場合と、構造体配列を使った場合を比較すると、構造体配列はアクセスがとんでもなく遅いよう。
ということで、構造体を作るときは注意しなければいけないな。使い方が悪いだけかもしれないけども。しかし、こんなに遅くなるとは思わなかった。

using struct: y = s.x .* x
Elapsed time is 0.000005 seconds.

using struct (for loop): y = s2.x .* x
Elapsed time is 0.000239 seconds.

using struct (index): y = s2.x .* x
Elapsed time is 0.001868 seconds.

使ったコードは下記のもの。

clear all;
n = 1000000;

s = struct('x', ones(1,n));
x = 2 * ones(1,n);
y = zeros(1,n);

% for loop
disp('for loop: y(i) = x(i) * x(i);');
tic
for i = 1:n
    y(i) = x(i) * x(i);
end
toc

% vectorized loop
disp('vectorized loop: y = x .* x');
tic
y = x .* x;
toc

% vectorized loop?
disp('vectorized loop2: y(1:n) = x(1:n) .* x(1:n)');
tic
y(1:n) = x(1:n) .* x(1:n);
toc

% vectorized loop using structure
disp('using struct: y = s.x .* x');
tic 
y = s.x .* x;
toc

clear all;
n = 1000;

s = struct('x', ones(1,n));  % s.x(1:n) のようにs.xに配列を指定
s2(1:n) = struct('x', 1);    % s(1:n).x のように構造体の配列
y = zeros(1,n);

disp('using struct: y = s.x .* x');
tic 
y = s.x .* x;
toc

disp('using struct (for loop): y = s2.x .* x');
tic 
s2_x = zeros(1,n);
for i = 1:n
    s2_x(i) = s2(i).x;
end
y = s2_x .* x;
toc

disp('using struct (index): y = s2.x .* x');
tic 
s2_x = zeros(1,n);
s2_x = s2(1:n).x;
y = s2_x .* x;
toc