Skip to content

Instantly share code, notes, and snippets.

@DSCF-1224

DSCF-1224/20200504_0007.f08

Created May 3, 2020 15:13
Show Gist options
  • Save DSCF-1224/64aee379283c12022f00764f5129c091 to your computer and use it in GitHub Desktop.
Save DSCF-1224/64aee379283c12022f00764f5129c091 to your computer and use it in GitHub Desktop.
Download ZIP
ガウス過程と機械学習 第1章 linear.py version 01
Raw
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
! ==================================================================================================================================
Raw
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
# ==================================================================================================================================
Raw
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment