Skip to content

Instantly share code, notes, and snippets.

@DSCF-1224

DSCF-1224/20200505_1430.f08

Created May 5, 2020 05:39
Show Gist options
  • Save DSCF-1224/0463fc73c05cd356b30a1fd5c0fa6953 to your computer and use it in GitHub Desktop.
Save DSCF-1224/0463fc73c05cd356b30a1fd5c0fa6953 to your computer and use it in GitHub Desktop.
Download ZIP
ガウス過程と機械学習 第1章 lm.py version 01 (Fortran)
Raw
20200505_1430.f08
! ==================================================================================================================================
! ISBN 978-4-06-152926-7
! ガウス過程と機械学習
!
! [reference]
! http://chasen.org/~daiti-m/gpbook/python/lm.py
! http://chasen.org/~daiti-m/gpbook/data/nonlinear.dat
! https://nbviewer.jupyter.org/gist/genkuroki/a37894d5669ad13b4cd27da16096bfd2
! https://gihyo.jp/book/2015/978-4-7741-7506-5
! ==================================================================================================================================
module mod_ISBN9784061529267
! <module>s to import
use, intrinsic :: iso_fortran_env
use, non_intrinsic :: ISBN9784774175065
! require all variables to be explicitly declared
implicit none
! accessibility of the <subroutine>s and <function>s in this <module>
public :: type_data_estimated ! type
public :: type_data_observed ! type
public :: type_data_regression ! type
public :: type_data_regression_degree ! type
public :: type_data_regression_weight ! type
public :: compute_weight ! subroutine
public :: estimate_output ! subroutine
public :: read_data ! subroutine
public :: save_result ! subroutine
private :: compute_estimated_output ! function
! <type>s for this <module>
type type_data_observed
integer(INT32), public :: size
real(REAL64), allocatable, public :: inpt(:)
real(REAL64), allocatable, public :: otpt(:)
end type type_data_observed
type type_data_regression_degree
integer(INT32), public :: poly
integer(INT32), public :: trig
end type type_data_regression_degree
type type_data_regression_weight
real(REAL64), allocatable, public :: poly(:)
real(REAL64), allocatable, public :: trig(:)
end type type_data_regression_weight
type type_data_regression
type(type_data_regression_degree), public :: degree
type(type_data_regression_weight), public :: weight
end type type_data_regression
type type_data_estimated
integer(INT32), public :: size
real(REAL64), allocatable, public :: inpt(:)
real(REAL64), allocatable, public :: otpt(:)
end type type_data_estimated
! contained <subroutine>s and <function>s are below
contains
pure function compute_estimated_output (data_regression, input) result(output)
! arguments for this <function>
type(type_data_regression), intent(in) :: data_regression
real(REAL64), intent(in) :: input
! return value of this <function>
real(REAL64) :: output
! variables for this <function>
integer(INT32) :: itr_dgr
itr_dgr = data_regression%degree%poly
output = data_regression%weight%poly(itr_dgr - 1) + data_regression%weight%poly(itr_dgr) * input
itr_dgr = itr_dgr - 2_INT32
do while (itr_dgr .ge. 0_INT32)
output = data_regression%weight%poly(itr_dgr) + input * output
itr_dgr = itr_dgr - 1_INT32
end do
do itr_dgr = 1_INT32, data_regression%degree%trig, 1_INT32
output = &!
output + &!
data_regression%weight%trig(2 * itr_dgr - 1) * cos(itr_dgr * input) + &!
data_regression%weight%trig(2 * itr_dgr ) * sin(itr_dgr * input)
end do
end function compute_estimated_output
subroutine compute_weight (data_degree, data_observed, data_regression)
! arguments for this <subroutine>
type(type_data_regression_degree), intent(in) :: data_degree
type(type_data_observed), intent(in) :: data_observed
type(type_data_regression), intent(inout) :: data_regression
! variables for this <subroutine>
integer(INT32) :: bffr_dgr
integer(INT32) :: bffr_stat
integer(INT32) :: itr_cl, itr_rw
integer(INT32), allocatable :: pivot(:)
real(REAL64), allocatable :: bffr_vector(:)
real(REAL64), allocatable :: bffr_matrix(:,:)
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
data_regression%degree = data_degree
allocate(data_regression%weight%poly(0:data_degree%poly), stat= bffr_stat)
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `data_regression%weight%poly(:)` !'
allocate(data_regression%weight%trig(1:2*data_degree%trig), stat= bffr_stat)
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `data_regression%weight%trig(:)` !'
! STEP.02
! allocate the arrays to store the design matrix
bffr_dgr = data_degree%poly + 2 * data_degree%trig
allocate(matrix_design_unit(1:data_observed%size, 0:bffr_dgr), stat= bffr_stat)
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `matrix_design_unit(:,:)` !'
allocate(matrix_design_prod(0:bffr_dgr, 0:bffr_dgr), stat= bffr_stat)
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `matrix_design_prod(:,:)` !'
allocate(pivot(0:bffr_dgr), stat= bffr_stat)
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `pivot(:)` !'
allocate(bffr_vector(0:bffr_dgr), stat= bffr_stat)
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `bffr_vector(:,:)` !'
allocate(bffr_matrix(0:bffr_dgr, 0:bffr_dgr), stat= bffr_stat)
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `bffr_matrix(:,:)` !'
! STEP.03
! compute the design matrix
matrix_design_unit(1:data_observed%size, 0) = 1
do itr_cl = 1_INT32, data_degree%poly, 1_INT32
do concurrent (itr_rw = 1:data_observed%size:1)
matrix_design_unit(itr_rw, itr_cl ) = &!
matrix_design_unit(itr_rw, itr_cl - 1) * data_observed%inpt(itr_rw)
end do
end do
do itr_cl = 1_INT32, data_degree%trig, 1_INT32
do concurrent (itr_rw = 1:data_observed%size:1)
matrix_design_unit(itr_rw, data_degree%poly + 2 * itr_cl - 1) = cos( itr_cl * data_observed%inpt(itr_rw) )
matrix_design_unit(itr_rw, data_degree%poly + 2 * itr_cl ) = sin( itr_cl * data_observed%inpt(itr_rw) )
end do
end do
! STEP.04
! compute the weight
matrix_design_prod = matmul( transpose(matrix_design_unit), matrix_design_unit )
call lu_decomposition_pivot(&!
dim = bffr_dgr + 1_INT32, &!
matrix = matrix_design_prod(0:bffr_dgr, 0:bffr_dgr), &!
pivot = pivot(0:bffr_dgr) &!
)
call lu_inverse_pivot(&!
dim = bffr_dgr + 1_INT32, &!
mtrx_trg = matrix_design_prod(0:bffr_dgr, 0:bffr_dgr), &!
pivot = pivot(0:bffr_dgr), &!
mtrx_tmp = bffr_matrix(0:bffr_dgr, 0:bffr_dgr), &!
vctr_tmp = bffr_vector(0:bffr_dgr) &!
)
bffr_vector(0:bffr_dgr) = matmul( matrix_design_prod, matmul( transpose(matrix_design_unit), data_observed%otpt ) )
data_regression%weight%poly(0: data_degree%poly) = bffr_vector(0:data_degree%poly)
data_regression%weight%trig(1:2*data_degree%trig) = bffr_vector(data_degree%poly + 1:bffr_dgr)
! STEP.04
! deallocate arrays used in this subroutine
deallocate(bffr_vector)
deallocate(bffr_matrix)
deallocate(matrix_design_prod)
deallocate(matrix_design_unit)
deallocate(pivot)
end subroutine compute_weight
subroutine estimate_output (size, data_observed, data_regression, data_estimated)
! arguments for this <subroutine>
integer(INT32), intent(in) :: size
type(type_data_observed), intent(in) :: data_observed
type(type_data_regression), intent(in) :: data_regression
type(type_data_estimated), intent(inout) :: data_estimated
! variables for this <subroutine>
integer(INT32) :: itr_pnt
integer(INT32) :: bffr_stat
real(REAL64) :: interval_inpt
real(REAL64) :: maxval_inpt
real(REAL64) :: minval_inpt
! STEP.01
! store the given number of points
data_estimated%size = size
! STEP.02
! allocate the arrays to store the estimated output
allocate( data_estimated%inpt(1:size), stat= bffr_stat )
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `data_estimated%inpt` !'
allocate( data_estimated%otpt(1:size), stat= bffr_stat )
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate `data_estimated%otpt` !'
! STEP.03
! prepare the values to compute input
minval_inpt = minval(array= data_observed%inpt(:))
maxval_inpt = maxval(array= data_observed%inpt(:))
interval_inpt = maxval_inpt - minval_inpt
minval_inpt = minval_inpt - 0.1_REAL64 * interval_inpt
maxval_inpt = maxval_inpt + 0.1_REAL64 * interval_inpt
interval_inpt = (maxval_inpt - minval_inpt) / ( size - 1_INT32 )
! STEP.04
! compute the input and estimated output
do concurrent (itr_pnt = 1_INT32 : size : 1_INT32)
data_estimated%inpt(itr_pnt) = minval_inpt + (itr_pnt - 1_INT32) * interval_inpt
data_estimated%otpt(itr_pnt) = compute_estimated_output( data_regression, data_estimated%inpt(itr_pnt) )
end do
end subroutine estimate_output
subroutine read_data (file, data_observed)
! arguments for this <subroutine>
character(len=*), intent(in) :: file ! file path to read
type(type_data_observed), intent(inout) :: data_observed ! 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
do
read(unit=unit_read, fmt='()', iostat=bffr_stat)
if ( IS_IOSTAT_END(bffr_stat) ) exit
bffr_size = bffr_size + 1_INT32
end do
data_observed%size = bffr_size
! STEP.03
! allocate the array to store the data
allocate( data_observed%inpt(1:bffr_size), stat= bffr_stat )
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate an array `data_observed%inpt` !'
allocate( data_observed%otpt(1:bffr_size), stat= bffr_stat )
if (bffr_stat .gt. 0_INT32) stop 'Failed to allocate an array `data_observed%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=*) &!
data_observed%inpt(itr_line), &!
data_observed%otpt(itr_line)
end do
! STEP.06
! close the file
close(unit= unit_read, status='keep')
end subroutine read_data
subroutine save_result (data_observed, data_regression, data_estimated)
! arguments for this <subroutine>
type(type_data_observed), intent(in) :: data_observed
type(type_data_regression), intent(in) :: data_regression
type(type_data_estimated), intent(in) :: data_estimated
! 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, data_observed%size, 1_INT32
write(unit= unit_save, fmt='(ES24.16E3,",",ES24.16E3)', advance='yes') &!
data_observed%inpt(itr), &!
data_observed%otpt(itr)
end do
write(unit= unit_save, fmt='(/,A,",",A)', advance='yes') &!
'degree', &!
'weight (poly)'
do itr = 0_INT32, data_regression%degree%poly, 1_INT32
write(unit= unit_save, fmt='(I2,",",ES24.16E3)', advance='yes') &!
itr, data_regression%weight%poly(itr)
end do
write(unit= unit_save, fmt='(/,A,",",A)', advance='yes') &!
'degree', &!
'weight (trig)'
do itr = 1_INT32, 2 * data_regression%degree%trig, 1_INT32
write(unit= unit_save, fmt='(I2,",",ES24.16E3)', advance='yes') &!
itr, data_regression%weight%trig(itr)
end do
write(unit= unit_save, fmt='(/,A,",",A)', advance='yes') &!
'input', &!
'output (estimated)'
do itr = 1_INT32, data_estimated%size, 1_INT32
write(unit= unit_save, fmt='(ES24.16E3,",",ES24.16E3)', advance='yes') &!
data_estimated%inpt(itr), &!
data_estimated%otpt(itr)
end do
close(unit= unit_save, status='keep')
end subroutine save_result
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_data_estimated) :: obj_data_estimated
type(type_data_observed) :: obj_data_observed
type(type_data_regression) :: obj_data_regression
type(type_data_regression_degree) :: obj_data_degree
obj_data_degree%poly = 2
obj_data_degree%trig = 1
call read_data('downloaded/chap01/nonlinear.dat', obj_data_observed)
call compute_weight(obj_data_degree, obj_data_observed, obj_data_regression)
call estimate_output(101, obj_data_observed, obj_data_regression, obj_data_estimated)
call save_result(obj_data_observed, obj_data_regression, obj_data_estimated)
end program linear
! ==================================================================================================================================
! EOF
! ==================================================================================================================================
Raw
20200505_1430.plt
# ==================================================================================================================================
# ISBN 978-4-06-152926-7
# ガウス過程と機械学習
#
# [reference]
# http://chasen.org/~daiti-m/gpbook/python/lm.py
# http://chasen.org/~daiti-m/gpbook/data/nonlinear.dat
# https://nbviewer.jupyter.org/gist/genkuroki/a37894d5669ad13b4cd27da16096bfd2
# https://gihyo.jp/book/2015/978-4-7741-7506-5
# ==================================================================================================================================
set datafile separator comma
set key on outside right center vertical Left reverse
set format x '%3.1f'
set format y '%3.1f'
set grid x
set grid y
set xlabel 'input'
set ylabel 'output'
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:3::3 with lines title 'regression'
# ==================================================================================================================================
# EOF
# ==================================================================================================================================
Raw
20200505_1430.sh
# GNU Fortran (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0
gfortran -c -O2 -std=f2008 -Wall -fcheck=all -pedantic-errors -ffree-line-length-none ../ISBN9784774175065/Module/ISBN9784774175065.f08
gfortran -c -O2 -std=f2008 -Wall -fcheck=all -pedantic-errors -ffree-line-length-none fortran2008/chap01/20200505_1430.f08
gfortran -o 20200505_1430.exe 20200505_1430.o ISBN9784774175065.o
Raw
result.csv
input output (observed)
2.7613622000000002E+000 7.8126940000000000E-001
-2.5020370999999999E+000 5.7840239999999998E-001
-6.5341980000000000E-001 -8.3648389999999995E-001
-5.0937080000000001E-001 -1.0659939000000000E+000
1.0698122999999999E+000 -5.0531780000000004E-001
1.2444854000000001E+000 6.5693199999999993E-002
-1.5838630000000001E-001 -1.7132896000000000E+000
-1.8188962000000000E+000 1.5392700000000001E-001
-3.6129367000000001E+000 1.0440645000000000E+000
-2.8263055000000001E+000 7.7416410000000002E-001
-1.0204458000000001E+000 -8.3045159999999996E-001
2.6113039999999998E-001 -1.3202885000000000E+000
-9.6677950000000001E-001 -1.0839772000000001E+000
4.7017170000000003E-001 -4.5600799999999997E-002
-7.2799290000000005E-001 -3.6915999999999997E-002
3.0501329999999999E-001 -1.1703207000000000E+000
3.9014330000000003E-001 -7.7879779999999998E-001
3.0430492000000000E+000 1.2799982999999999E+000
-1.6307559000000000E+000 -2.2101540000000000E-001
-3.5160841999999999E+000 5.7242400000000004E-001
-6.4898370999999999E+000 -4.1269260000000002E-001
1.8201852000000001E+000 9.1200530000000002E-001
-7.7676290000000003E-001 -8.3715919999999999E-001
-3.4106993999999999E+000 6.3177320000000003E-001
-9.1959919999999995E-001 -1.2476601000000000E+000
-1.1248267000000001E+000 -4.3716270000000002E-001
-1.5116411999999999E+000 -7.7745160000000002E-001
1.8203301999999999E+000 1.2516906999999999E+000
-1.1185539000000000E+000 -1.6227400999999999E+000
degree weight (poly)
0 -9.5426604458734299E-002
1 7.2322955596856808E-002
2 2.4200714833955781E-002
degree weight (trig)
1 -9.0875410400920353E-001
2 2.9939963775084627E-001
input output (estimated)
-7.4431257300000002E+000 6.9552597830008667E-002
-7.2143364588000001E+000 -1.4027833846537097E-001
-6.9855471875999999E+000 -3.0678636857615749E-001
-6.7567579163999998E+000 -4.2452885579149924E-001
-6.5279686452000005E+000 -4.9047261765132050E-001
-6.2991793740000004E+000 -5.0415200510174263E-001
-6.0703901028000002E+000 -4.6769317054032072E-001
-5.8416008316000001E+000 -3.8570325998311028E-001
-5.6128115604000000E+000 -2.6503030428100305E-001
-5.3840222892000007E+000 -1.1440632303864626E-001
-5.1552330180000006E+000 5.6007758553942744E-002
-4.9264437468000004E+000 2.3515166258984493E-001
-4.6976544756000003E+000 4.1164218251905055E-001
-4.4688652044000001E+000 5.7436643716501989E-001
-4.2400759332000000E+000 7.1306103630202688E-001
-4.0112866619999998E+000 8.1884697365285619E-001
-3.7824973908000001E+000 8.8469237050273142E-001
-3.5537081196000000E+000 9.0577895330588232E-001
-3.3249188484000003E+000 8.7975316571811313E-001
-3.0961295772000001E+000 8.0684882795014645E-001
-2.8673403060000000E+000 6.8987495085335016E-001
-2.6385510347999999E+000 5.3406933982626303E-001
-2.4097617635999997E+000 3.4682561820848634E-001
-2.1809724924000005E+000 1.3730789677439117E-001
-1.9521832212000003E+000 -8.4026828557700384E-002
-1.7233939500000002E+000 -3.0597365979354668E-001
-1.4946046788000000E+000 -5.1716375804396542E-001
-1.2658154075999999E+000 -7.0665684874777057E-001
-1.0370261363999997E+000 -8.6451139157927481E-001
-8.0823686520000049E-001 -9.8230269972365425E-001
-5.7944759400000034E-001 -1.0535620059054203E+000
-3.5065832280000020E-001 -1.0741135925416412E+000
-1.2186905160000006E-001 -1.0422924159427531E+000
1.0692021960000009E-001 -9.5903088272764547E-001
3.3570949080000023E-001 -8.2781025596135172E-001
5.6449876199999949E-001 -6.5447922356220645E-001
7.9328803319999963E-001 -4.4694908457612104E-001
1.0220773043999998E+000 -2.1478143917284798E-001
1.2508665755999999E+000 3.1310129440453860E-002
1.4796558468000001E+000 2.8001841054421850E-001
1.7084451180000002E+000 5.2003185927877627E-001
1.9372343892000004E+000 7.4062411589730281E-001
2.1660236604000005E+000 9.3221302716061816E-001
2.3948129316000006E+000 1.0868600355797073E+000
2.6236022028000008E+000 1.1986839100368205E+000
2.8523914739999991E+000 1.2641672555520851E+000
3.0811807451999993E+000 1.2823398279488851E+000
3.3099700163999994E+000 1.2548290996810072E+000
3.5387592875999996E+000 1.1857754415026061E+000
3.7675485587999997E+000 1.0816163404239116E+000
3.9963378299999999E+000 9.5075089977906180E-001
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment