ガウス過程と機械学習 第1章 linear.py version 01

20200504_0007.f08

! ==================================================================================================================================
! ISBN 978-4-06-152926-7
! ガウス過程と機械学習
! 
! [reference]
! http://chasen.org/~daiti-m/gpbook/python/linear.py
! http://chasen.org/~daiti-m/gpbook/data/linear.dat
! https://nbviewer.jupyter.org/gist/genkuroki/a37894d5669ad13b4cd27da16096bfd2
! http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/
! ==================================================================================================================================

module mod_ISBN9784061529267

    ! <module>s to import
    use, intrinsic :: iso_fortran_env

    ! require all variables to be explicitly declared
    implicit none


    ! accessibility of the <subroutine>s and <function>s in this <module>
    public  :: type_learned_data        ! type
    public  :: type_obsereved_data      ! type
    public  :: compute_estimated_output ! function
    public  :: compute_weight           ! subroutine
    public  :: estimate_output          ! subroutine
    public  :: read_data                ! subroutine
    public  :: save_result              ! subroutine
    public  :: show_learned_data        ! subroutine
    public  :: show_observed_data       ! subroutine
    private :: inverse_matrix           ! subroutine


    ! <type>s for this <module>
    type type_obsereved_data
        integer(INT32),             public :: size
        real(REAL64),  allocatable, public :: inpt(:)
        real(REAL64),  allocatable, public :: otpt(:)
    end type type_obsereved_data

    type type_learned_data
        integer(INT32),             public :: degree     ! the degree of the regression
        integer(INT32),             public :: num_points
        real(REAL64),  allocatable, public :: weight(:)
        real(REAL64),  allocatable, public :: inpt(:)
        real(REAL64),  allocatable, public :: otpt(:)
    end type type_learned_data


    ! contained <subroutine>s and <function>s are below
    contains


    pure function compute_estimated_output (obj_learned_data, val_input) result(val_output)

        ! arguments for this <function>
        type(type_learned_data), intent(in) :: obj_learned_data
        real(REAL64),            intent(in) :: val_input

        ! return value of this <function>
        real(REAL64) :: val_output

        ! variables for this <function>
        integer(INT32) :: itr_dgr


        itr_dgr    = obj_learned_data%degree
        val_output = obj_learned_data%weight(itr_dgr - 1) + obj_learned_data%weight(itr_dgr) * val_input
        itr_dgr    = itr_dgr - 2_INT32

        do while (itr_dgr .ge. 0_INT32)
            val_output = obj_learned_data%weight(itr_dgr) + val_input * val_output
            itr_dgr    = itr_dgr - 1_INT32
        end do

    end function compute_estimated_output


    subroutine compute_weight (degree, obj_observed_data, obj_learned_data)

        ! arguments for this <subroutine>
        integer(INT32),            intent(in)    :: degree
        type(type_obsereved_data), intent(in)    :: obj_observed_data
        type(type_learned_data),   intent(inout) :: obj_learned_data


        ! variables for this <subroutine>
        integer(INT32)             :: bffr_stat
        integer(INT32)             :: itr_cl, itr_rw
        real(REAL64),  allocatable :: matrix_design_unit(:,:)
        real(REAL64),  allocatable :: matrix_design_prod(:,:)


        ! STEP.01
        ! store the given degree of regression
        ! allocate the array to store the weight of regression
        obj_learned_data%degree = degree

        allocate(obj_learned_data%weight(0:degree), stat= bffr_stat)
        if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `obj_learned_data%weight` !'


        ! STEP.02
        ! allocate the arrays to store the design matrix
        allocate(matrix_design_unit(1:obj_observed_data%size, 0:degree), stat= bffr_stat)
        if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `matrix_design_unit` !'

        allocate(matrix_design_prod(0:degree, 0:degree), stat= bffr_stat)
        if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `matrix_design_prod` !'


        ! STEP.03
        ! compute the design matrix
        matrix_design_unit(1:obj_observed_data%size, 0) = 1.0_REAL64

        do itr_cl = 1_INT32, degree, 1_INT32
            do concurrent (itr_rw = 1:obj_observed_data%size:1)
                matrix_design_unit(itr_rw, itr_cl    ) = &!
                matrix_design_unit(itr_rw, itr_cl - 1) * obj_observed_data%inpt(itr_rw)
            end do
        end do


        ! STEP.04
        ! compute the weight
        matrix_design_prod = matmul( transpose(matrix_design_unit), matrix_design_unit )
        call inverse_matrix(degree + 1_INT32, matrix_design_prod)
        obj_learned_data%weight = matmul( matrix_design_prod, matmul( transpose(matrix_design_unit), obj_observed_data%otpt ) )


        ! STEP.04
        ! deallocate arrays used in this subroutine
        deallocate(matrix_design_prod)
        deallocate(matrix_design_unit)

    end subroutine compute_weight


    subroutine estimate_output (num_points, obj_observed_data, obj_learned_data)

        ! arguments for this <subroutine>
        integer(INT32),            intent(in)    :: num_points
        type(type_obsereved_data), intent(in)    :: obj_observed_data
        type(type_learned_data),   intent(inout) :: obj_learned_data

        ! variables for this <subroutine>
        integer(INT32) :: itr_pnt
        integer(INT32) :: bffr_stat
        real(REAL64)   :: interval_inpt
        real(REAL64)   :: minval_inpt


        ! STEP.01
        ! store the given number of points
        obj_learned_data%num_points = num_points


        ! STEP.02
        ! allocate the arrays to store the estimated output
        allocate( obj_learned_data%inpt(1:num_points), stat= bffr_stat )
        if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `obj_learned_data%inpt` !'

        allocate( obj_learned_data%otpt(1:num_points), stat= bffr_stat )
        if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `obj_learned_data%otpt` !'


        ! STEP.03
        ! prepare the values to compute input
        minval_inpt = minval(array= obj_observed_data%inpt(:))

        interval_inpt &!
        = ( maxval(array= obj_observed_data%inpt(:)) - minval_inpt ) &!
        / ( num_points - 1_INT32 )


        ! STEP.04
        ! compute the input and estimated output
        do concurrent (itr_pnt = 1_INT32 : num_points : 1_INT32)
            obj_learned_data%inpt(itr_pnt) = minval_inpt + (itr_pnt - 1_INT32) * interval_inpt
            obj_learned_data%otpt(itr_pnt) = compute_estimated_output( obj_learned_data, obj_learned_data%inpt(itr_pnt) )
        end do


    end subroutine estimate_output


    subroutine inverse_matrix (dim, matrix)

        ! arguments for this <subroutine>
        integer(INT32), intent(in)    :: dim
        real(REAL64),   intent(inout) :: matrix(1:dim, 1:dim)

        ! variables for this <subroutine>
        integer(INT32)              :: row_pv
        integer(INT32)              :: bffr_stat
        integer(INT32)              :: bffr_val_swap_i
        integer(INT32)              :: itr_dm, itr_cl, itr_rw
        integer(INT32), allocatable :: idx_rw(:)
        real(REAL64)                :: val_pv
        real(REAL64)                :: bffr_val_swap_r


        ! STEP.01
        ! allocate the array to store the index of column
        allocate(idx_rw(1:dim), stat= bffr_stat)
        if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate array `idx_rw` !'


        ! STEP.02
        ! initialize the array to store the index of column
        do itr_rw = 1_INT32, dim, 1_INT32
            idx_rw(itr_rw) = itr_rw
        end do


        ! STEP.03
        ! compute the invertible matrix
        do itr_dm = 1_INT32, dim, 1_INT32

            ! STEP.01
            ! select the pivot
            row_pv = maxloc(array= abs(matrix(itr_dm:dim, itr_dm)), dim= 1) + itr_dm - 1_INT32
            val_pv = 1 / matrix(row_pv, itr_dm)

            ! STEP.02
            ! swap the rows of the coefficient matrix for pivot selection
            if (itr_dm .ne. row_pv) then

                do concurrent (itr_cl = itr_dm : dim)
                    bffr_val_swap_r        = matrix(itr_dm, itr_cl)
                    matrix(itr_dm, itr_cl) = matrix(row_pv, itr_cl)
                    matrix(row_pv, itr_cl) = bffr_val_swap_r
                end do

                bffr_val_swap_i = idx_rw(itr_dm)
                idx_rw(itr_dm)  = idx_rw(row_pv)
                idx_rw(row_pv)  = bffr_val_swap_i

            end if

            ! STEP.03
            ! elimination
            matrix(itr_dm, itr_dm) = 1
            matrix(itr_dm, 1:dim)  = matrix(itr_dm, 1:dim) * val_pv

            do itr_rw = 1_INT32, dim, 1_INT32

                if (itr_rw .ne. itr_dm) then

                    val_pv                 = matrix(itr_rw, itr_dm)
                    matrix(itr_rw, itr_dm) = 0

                    do concurrent (itr_cl = 1_INT32 : dim : 1_INT32)
                        matrix(itr_rw, itr_cl) = &!
                        matrix(itr_rw, itr_cl) - &!
                        matrix(itr_dm, itr_cl) * val_pv
                    end do

                end if

            end do

        end do


        ! STEP.04
        ! restore the order of rows
        do itr_cl = 1_INT32, dim, 1_INT32

            if (idx_rw(itr_cl) .eq. itr_cl) cycle

            itr_dm = itr_cl

            do while (idx_rw(itr_dm) .ne. itr_cl)
                itr_dm = itr_dm + 1_INT32
            end do

            do concurrent (itr_rw = 1_INT32 : dim : 1_INT32)
                bffr_val_swap_r        = matrix(itr_rw, itr_dm)
                matrix(itr_rw, itr_dm) = matrix(itr_rw, itr_cl)
                matrix(itr_rw, itr_cl) = bffr_val_swap_r
            end do

            idx_rw(itr_dm) = idx_rw(itr_cl)

        end do


        ! STEP.05
        ! deallocate the array
        deallocate(idx_rw)

    end subroutine inverse_matrix


    subroutine read_data (file, obj_observed_data)

        ! arguments for this <subroutine>
        character(len=*),          intent(in)    :: file              ! file path to read
        type(type_obsereved_data), intent(inout) :: obj_observed_data ! object to store the read data


        ! variables for this <subroutine>
        integer(INT32) :: bffr_size
        integer(INT32) :: bffr_stat
        integer(INT32) :: unit_read
        integer(INT32) :: itr_line


        ! STEP.01
        ! open the target file
        open(&!
            action  = 'read',      &!
            file    = file,        &!
            form    = 'formatted', &!
            iostat  = bffr_stat,   &!
            newunit = unit_read,   &!
            status  = 'old'        &!
        )

        if (bffr_stat .gt. 0_INT32) then

            write(unit=OUTPUT_UNIT, fmt='(2(A,/),A,I0)') &!
                'Failed to open file !', &!
                'PATH   > ' // file, &!
                'IOSTAT > ', bffr_stat

            stop

        end if


        ! STEP.02
        ! check the given data size
        bffr_size = 0_INT32

        loop_count_size: do
            read(unit=unit_read, fmt='()', iostat=bffr_stat)
            if ( IS_IOSTAT_END(bffr_stat) ) exit loop_count_size
            bffr_size = bffr_size + 1_INT32
        end do loop_count_size

        obj_observed_data%size = bffr_size


        ! STEP.03
        ! allocate the array to store the data
        allocate( obj_observed_data%inpt(1:bffr_size), stat= bffr_stat )
        if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate an array `inpt` !'

        allocate( obj_observed_data%otpt(1:bffr_size), stat= bffr_stat )
        if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate an array `otpt` !'


        ! STEP.04
        ! store the read data
        rewind(unit= unit_read, iostat= bffr_stat)
        if (bffr_stat .gt. 0_INT32) stop 'Failed to rewind !'


        ! STEP.05
        ! read out the data
        do itr_line = 1_INT32, bffr_size, 1_INT32
            read(unit=unit_read, fmt=*) &!
                obj_observed_data%inpt(itr_line), &!
                obj_observed_data%otpt(itr_line)
        end do


        ! STEP.06
        ! close the file
        close(unit= unit_read, status='keep')


    end subroutine read_data


    subroutine save_result (obj_observed_data, obj_learned_data)

        ! arguments for this <subroutine>
        type(type_obsereved_data), intent(in) :: obj_observed_data
        type(type_learned_data),   intent(in) :: obj_learned_data

        ! variables for this <subroutine>
        integer(INT32) :: itr
        integer(INT32) :: bffr_stat
        integer(INT32) :: unit_save


        open(&!
            action  = 'write',                         &!
            file    = 'fortran2008/chap01/result.csv', &!
            form    = 'formatted',                     &!
            iostat  = bffr_stat,                       &!
            newunit = unit_save,                       &!
            status  = 'replace'                        &!
        )


        write(unit= unit_save, fmt='(A,",",A)', advance='yes') &!
            'input', &!
            'output (observed)'

        do itr = 1_INT32, obj_observed_data%size, 1_INT32
            write(unit= unit_save, fmt='(ES24.16E3,",",ES24.16E3)', advance='yes') &!
                obj_observed_data%inpt(itr), &!
                obj_observed_data%otpt(itr)
        end do


        write(unit= unit_save, fmt='(/,A,",",A)', advance='yes') &!
            'input', &!
            'output (estimated)'

        do itr = 1_INT32, obj_learned_data%num_points, 1_INT32
            write(unit= unit_save, fmt='(ES24.16E3,",",ES24.16E3)', advance='yes') &!
                obj_learned_data%inpt(itr), &!
                obj_learned_data%otpt(itr)
        end do


        close(unit= unit_save, status='keep')

    end subroutine save_result


    subroutine show_learned_data (obj_learned_data)

        ! arguments for this <subroutine>
        type(type_learned_data), intent(in) :: obj_learned_data

        ! variables for this <subroutine>
        integer(INT32) :: itr

        write(unit= OUTPUT_UNIT, fmt='(/,A)')       'weight'
        write(unit= OUTPUT_UNIT, fmt='(A6,1X,A23)') 'degree', 'value'

        do itr = 0_INT32, obj_learned_data%degree, 1_INT32
            write(unit= OUTPUT_UNIT, fmt='(I6,1X,ES23.15E3)', advance='yes') &!
                itr, obj_learned_data%weight(itr)
        end do

    end subroutine show_learned_data


    subroutine show_observed_data (obj_observed_data)

        ! arguments for this <subroutine>
        type(type_obsereved_data), intent(in) :: obj_observed_data

        ! support variables for this <subroutine>
        integer(INT32) :: itr

        write(unit= OUTPUT_UNIT, fmt='(A10,2(1X,A23))', advance='yes') 'No.', 'input', 'output'

        do itr = 1_INT32, obj_observed_data%size, 1_INT32
            write(unit= OUTPUT_UNIT, fmt= '(I10,2(1X,ES23.15E3))', advance= 'yes') &!
                itr, &!
                obj_observed_data%inpt(itr), &!
                obj_observed_data%otpt(itr)
        end do

    end subroutine show_observed_data
    
end module mod_ISBN9784061529267


! ==================================================================================================================================


program linear

    ! <module>s to import
    use,     intrinsic :: iso_fortran_env
    use, non_intrinsic :: mod_ISBN9784061529267

    ! require all variables to be explicitly declared
    implicit none

    ! variables for this <program>
    type(type_learned_data)   :: obj_learned_data
    type(type_obsereved_data) :: obj_observed_data


    call read_data('downloaded/chap02/linear.dat', obj_observed_data)
    call compute_weight(2, obj_observed_data, obj_learned_data)
    call show_learned_data(obj_learned_data)
    call estimate_output(101, obj_observed_data, obj_learned_data)
    call save_result (obj_observed_data, obj_learned_data)

end program linear


! ==================================================================================================================================
! EOF
! ==================================================================================================================================

20200504_0007.plt

# ==================================================================================================================================
# ISBN 978-4-06-152926-7
# ガウス過程と機械学習
# 
# [reference]
# http://chasen.org/~daiti-m/gpbook/python/linear.py
# http://chasen.org/~daiti-m/gpbook/data/linear.dat
# https://nbviewer.jupyter.org/gist/genkuroki/a37894d5669ad13b4cd27da16096bfd2
# http://www.caero.mech.tohoku.ac.jp/publicData/Daiguji/
# ==================================================================================================================================

set datafile separator comma
set key on outside right center vertical Left reverse
set xrange [-2.5:3.5]
set xzeroaxis linetype -1
set yzeroaxis linetype -1

plot 'result.csv' using 1:2 every ::1:0::0 with points title 'observed'   , \
     'result.csv' using 1:2 every ::1:1::1 with lines  title 'regression'

# ==================================================================================================================================
# EOF
# ==================================================================================================================================

20200504_0007.sh

gfortran -c -O2 -std=f2008 -Wall -fcheck=all -pedantic-errors -ffree-line-length-none 20200504_0007.f08
gfortran -o 20200504_0007.exe 20200504_0007.o