Fortranで機械学習がしたい：FortranからのPyTorchの利用

この記事[FortranからPythonを使いたい！:forpy](https://qiita.com/cometscome_phys/items/5b4347a48fb25b52e472)の続編です。
Fortranで計算して結果をそのまま機械学習したい、という場合に、forpyを利用することでFortranからPythonを呼び出してそのまま計算することが可能です。
forpyについては前の記事をみてください。

また、PyTorchについては[PyTorchで関数フィッティング](https://qiita.com/cometscome_phys/items/61331dd2eedd0c846336)を参考にしてください。この記事にあるPythonコードをFortranから実行するのが目的です。

# コード
まず、結論から述べます。Fortranコードのメイン部分はこんな感じになりました。

```fortran
subroutine test()
    use forpy_mod
    use pytorch
    implicit none
    interface 
        subroutine make_phi(phi,x0,n,k)
            implicit none
            real(4),intent(in)::x0(:)
            real(4),intent(out)::phi(:,:)
            integer,intent(in)::n,k
        end subroutine
    end interface
    integer::n
    real(8)::a0,a1,b0
    real(4),allocatable,asynchronous::x0(:),y0(:,:)
    real(8)::dx
    integer::d_input,d_middle,d_out
    type(object) :: net
    real(4),allocatable::phi(:,:)
    real(8)::learning_rate
    integer::nt
    integer::i,k
    integer::ierror

    ierror = forpy_initialize()

    n = 10
    a0 = 3d0
    a1 = 2d0
    b0 = 1d0
    allocate(y0(n,1),x0(n))

    dx = 4d0/dble(n-1)
    do i=1,n
        x0(i) = dble(i-1)*dx -2d0
    end do 

    y0(:,1) = a0*x0+a1*x0**2 + b0 + 3*cos(20*x0) 

    k = 4
    allocate(phi(n,k))
    call make_phi(phi,x0,n,k)

    d_input = k
    d_middle = 10
    d_out = 1

    call construct_neuralnet(net,d_input,d_middle,d_out)

    learning_rate = 0.001d0
    nt = 20000
    call train(net,phi,y0,learning_rate,nt)


    call forpy_finalize

end subroutine

subroutine make_phi(phi,x0,n,k)
    implicit none
    real(4),intent(in)::x0(:)
    real(4),intent(out)::phi(:,:)
    integer,intent(in)::n,k
    integer::i

    do i=0,k-1
        phi(:,i+1) = x0**i
    end do

    return

end subroutine

```
これを見ると、ただFortranで配列を作成しているだけなのがわかります。Pythonを呼び出しているのは```construct_neuralnet```と```train```です。この二つは```pytorch```というmoduleに定義されていまして、

```fortran
module pytorch
    use forpy_mod
    implicit none

    contains

    subroutine construct_neuralnet(net,d_inp,d_middle,d_out)
        integer,intent(in)::d_inp,d_middle,d_out
        type(object),intent(out) :: net
        integer::ierror
        type(tuple) :: args
        type(module_py) :: torch_sample
        type(list) :: paths

        ierror = get_sys_path(paths)
        ierror = paths%append(".")


        ierror = tuple_create(args, 3) !PythonのNetの引数が3つなので3
        ierror = args%setitem(0, d_inp)
        ierror = args%setitem(1, d_middle)
        ierror = args%setitem(2, d_out)
        ierror = import_py(torch_sample, "torchsample") 
        ierror = call_py(net, torch_sample, "Net", args) !Pythonの関数の引数はargsにまとめてある

        return
    end subroutine

    subroutine train(net,phi,y0,learning_rate,nt)
        type(object)::net
        integer::nt
        type(list) :: paths
        integer::ierror
        type(tuple) :: args
        type(module_py) :: trainfunc
        type(ndarray) ::y0_py,phi_py
        real(8)::learning_rate
        real(4),asynchronous,intent(in)::phi(:,:),y0(:,:)

        ierror = get_sys_path(paths)
        ierror = paths%append(".")

        ierror = import_py(trainfunc, "torchsample")

        ierror = ndarray_create_nocopy(phi_py, phi)
        ierror = ndarray_create_nocopy(y0_py, y0)

        ierror = tuple_create(args, 5) !Pythonの関数trainが5個の引数を持つので5
        ierror = args%setitem(0, net)
        ierror = args%setitem(1, learning_rate)
        ierror = args%setitem(2, phi_py)
        ierror = args%setitem(3, y0_py)
        ierror = args%setitem(4, nt)
        write(*,*) "training"

        !train(net,learning_rate,phi,y0,nt):
        ierror = call_py_noret(trainfunc,"train",args)

        return
    end subroutine
end module
```
となります。```construct_neuralnet```はPythonのオブジェクトである```net```を返します。そして、```train```はPythonコード```train```を実行します。この```net```と```train```は別のファイルで

```python:torchsample.py
import torch 
import torch.nn as nn 
import torch.nn.functional as F
from torch.autograd import Variable
import torch.optim as optim

class Net(nn.Module):
    def __init__(self,d_inp,d_middle,d_out):
        super(Net, self).__init__()
        self.fc1 = nn.Linear(d_inp,d_middle)
        self.fc2 = nn.Linear(d_middle, d_out)

    def forward(self, x):
        x = F.relu(self.fc1(x))
        x = self.fc2(x)        
        return x

def train(net,learning_rate,phi,y0,nt):
    print(net)
    params = list(net.parameters())
    print(params)

    inputdata = torch.from_numpy(phi)
    target = torch.from_numpy(y0)
    criterion = nn.MSELoss()
    optimizer = optim.Adam(net.parameters(), lr=learning_rate)

    print(nt)
    for i in range(nt):
        
        optimizer.zero_grad() #勾配の初期化？　必要らしい
        output = net(inputdata) #アウトプットの生成
        
        loss = criterion(output, target) #loss関数
        loss.backward() 
        optimizer.step()
        if i %1000 is 0:
            print(i,loss.item())
```
で定義されています。ですので、このPythonのファイルとFortranのファイルを用意すれば実行できます。

# コードの全体

Fortranコードの全体は

```fortran:fortranpytorch.f90
module pytorch
    use forpy_mod
    implicit none

    contains

    subroutine construct_neuralnet(net,d_inp,d_middle,d_out)
        integer,intent(in)::d_inp,d_middle,d_out
        type(object),intent(out) :: net
        integer::ierror
        type(tuple) :: args
        type(module_py) :: torch_sample
        type(list) :: paths

        ierror = get_sys_path(paths)
        ierror = paths%append(".")


        ierror = tuple_create(args, 3)
        ierror = args%setitem(0, d_inp)
        ierror = args%setitem(1, d_middle)
        ierror = args%setitem(2, d_out)
        ierror = import_py(torch_sample, "torchsample")
        ierror = call_py(net, torch_sample, "Net", args)

        return
    end subroutine

    subroutine train(net,phi,y0,learning_rate,nt)
        type(object)::net
        integer::nt
        type(list) :: paths
        integer::ierror
        type(tuple) :: args
        type(module_py) :: trainfunc
        type(ndarray) ::y0_py,phi_py
        real(8)::learning_rate
        real(4),asynchronous,intent(in)::phi(:,:),y0(:,:)

        ierror = get_sys_path(paths)
        ierror = paths%append(".")

        ierror = import_py(trainfunc, "torchsample")

        ierror = ndarray_create_nocopy(phi_py, phi)
        ierror = ndarray_create_nocopy(y0_py, y0)

        ierror = tuple_create(args, 5)
        ierror = args%setitem(0, net)
        ierror = args%setitem(1, learning_rate)
        ierror = args%setitem(2, phi_py)
        ierror = args%setitem(3, y0_py)
        ierror = args%setitem(4, nt)
        write(*,*) "training"

        !train(net,learning_rate,phi,y0,nt):

        ierror = call_py_noret(trainfunc,"train",args)

        !train(net,nt,optimizer,criterion,input,target)
        return
    end subroutine
end module

subroutine make_phi(phi,x0,n,k)
    implicit none
    real(4),intent(in)::x0(:)
    real(4),intent(out)::phi(:,:)
    integer,intent(in)::n,k
    integer::i

    do i=0,k-1
        phi(:,i+1) = x0**i
    end do

    return

end subroutine

subroutine test()
    use forpy_mod
    use pytorch
    implicit none
    interface 
        subroutine make_phi(phi,x0,n,k)
            implicit none
            real(4),intent(in)::x0(:)
            real(4),intent(out)::phi(:,:)
            integer,intent(in)::n,k
        end subroutine
    end interface
    integer::n
    real(8)::a0,a1,b0
    real(4),allocatable,asynchronous::x0(:),y0(:,:)
    real(8)::dx
    integer::d_input,d_middle,d_out
    type(object) :: net
    real(4),allocatable::phi(:,:)
    real(8)::learning_rate
    integer::nt
    integer::i,k
    integer::ierror

    ierror = forpy_initialize()

    n = 10
    a0 = 3d0
    a1 = 2d0
    b0 = 1d0
    allocate(y0(n,1),x0(n))

    dx = 4d0/dble(n-1)
    do i=1,n
        x0(i) = dble(i-1)*dx -2d0
    end do 

    y0(:,1) = a0*x0+a1*x0**2 + b0 + 3*cos(20*x0) 

    k = 4
    allocate(phi(n,k))
    call make_phi(phi,x0,n,k)

    d_input = k
    d_middle = 10
    d_out = 1

    call construct_neuralnet(net,d_input,d_middle,d_out)

    learning_rate = 0.001d0
    nt = 20000
    call train(net,phi,y0,learning_rate,nt)


    call forpy_finalize

end subroutine

program main
    implicit none
    call test()
end program main
```
でして、Pythonコードは上にある通りです。

コンパイルは

```
gfortran -c forpy_mod.F90 
gfortran fortranpytorch.f90 forpy_mod.o `python3-config --ldflags`  -L/usr/local/Cellar/gettext/0.21/lib
```
でできます。詳しくは[前の記事](https://qiita.com/cometscome_phys/items/5b4347a48fb25b52e472)をみてください。

# 実行結果
実行すると、

```
 training
Net(
  (fc1): Linear(in_features=4, out_features=10, bias=True)
  (fc2): Linear(in_features=10, out_features=1, bias=True)
)
[Parameter containing:
tensor([[-0.4106, -0.3108, -0.4385, -0.2647],
        [ 0.0301,  0.4899, -0.2945, -0.1717],
        [ 0.4867, -0.2890,  0.4036,  0.1898],
        [ 0.2753, -0.2922,  0.0601,  0.2919],
        [-0.2123,  0.3228,  0.0421, -0.3196],
        [ 0.1464, -0.4406, -0.0427,  0.2106],
        [-0.4179,  0.2699, -0.1437,  0.3141],
        [-0.0273, -0.4634, -0.4223, -0.4201],
        [-0.3101,  0.4756, -0.1964, -0.3231],
        [-0.3890,  0.1310, -0.2351,  0.1514]], requires_grad=True), Parameter containing:
tensor([ 0.0883,  0.1841,  0.4362,  0.2545,  0.0261,  0.3895, -0.1527,  0.0191,
         0.1739,  0.0740], requires_grad=True), Parameter containing:
tensor([[-0.1873,  0.1575,  0.1746, -0.1483,  0.0863,  0.1394,  0.3049,  0.0600,
         -0.1203, -0.1196]], requires_grad=True), Parameter containing:
tensor([0.2089], requires_grad=True)]
20000
0 37.118621826171875
1000 4.937760353088379
2000 3.480103015899658
3000 3.293708324432373
4000 3.213000774383545
5000 3.121199607849121
6000 2.9794585704803467
7000 2.6786561012268066
8000 1.734819769859314
9000 0.9694552421569824
10000 0.6194493770599365
11000 0.4148186147212982
12000 0.2840990424156189
13000 0.17756792902946472
14000 0.09583604335784912
15000 0.04160896688699722
16000 0.012834017165005207
17000 0.002260125009343028
18000 0.0001653230719966814
19000 7.0691885412088595e-06
```
こんな結果が得られます。