diff --git a/CHANGELOG b/CHANGELOG index 64f9109..2dfcd34 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -15,6 +15,43 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Removed +## [1.3.2] - 2024-06-28 + +### Added +- Added conv1d_layer_type +- Added avgpool1d_layer_type +- Added maxpool1d_layer_type +- Added input2d_layer_type +- Added statement of need to README + +### Changed + +### Fixed +- Fix reference to fpm install in README +- Fixed unallocated lr_decay reference +- Fixed ifort and ifx compiler reference to temporary arrays +- Fixed temporary array reference in tests +- Fixed epsilon undefined value in mod_loss +- Fixed project link in CMakeLists.txt + +### Removed + +## [1.3.1] - 2024-05-25 + +### Added + +### Changed +- Add reference to gcc 12.3 compatibility in README + +### Fixed +- Fix attempt to assign size of random seed in random_setup +- Fix attempt to assign size of random seed in test_network +- Fix attempt to assing size of random seed in example/simple +- Fix abstract interface use in base_layer_type +- Add appropriate references to neural-fortran in CONTRIBUTING.md and associated procedures, and modules + +### Removed + ## [1.3.0] - 2024-03-13 ### Added diff --git a/CMakeLists.txt b/CMakeLists.txt index 54a260d..cae0c77 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,7 +15,7 @@ project(athena NONE) set( LIB_NAME ${PROJECT_NAME} ) set( PROJECT_DESCRIPTION "Fortran neural network" ) -set( PROJECT_URL "https://https://git.exeter.ac.uk/hepplestone/athena" ) +set( PROJECT_URL "https://github.com/nedtaylor/athena" ) set( CMAKE_CONFIGURATION_TYPES "Release" "Parallel" "Serial" "Dev" "Debug" "Parallel_Dev" CACHE STRING "List of configurations types." ) set( CMAKE_BUILD_TYPE "Release" @@ -97,13 +97,16 @@ set(LIB_FILES mod_batchnorm1d_layer.f90 mod_batchnorm2d_layer.f90 mod_batchnorm3d_layer.f90 + mod_conv1d_layer.f90 mod_conv2d_layer.f90 mod_conv3d_layer.f90 mod_dropout_layer.f90 mod_dropblock2d_layer.f90 mod_dropblock3d_layer.f90 + mod_avgpool1d_layer.f90 mod_avgpool2d_layer.f90 mod_avgpool3d_layer.f90 + mod_maxpool1d_layer.f90 mod_maxpool2d_layer.f90 mod_maxpool3d_layer.f90 mod_full_layer.f90 @@ -112,6 +115,7 @@ set(LIB_FILES mod_flatten3d_layer.f90 mod_flatten4d_layer.f90 mod_input1d_layer.f90 + mod_input2d_layer.f90 mod_input3d_layer.f90 mod_input4d_layer.f90 mod_container_layer.f90 @@ -221,6 +225,7 @@ if(BUILD_TESTS) endif() # add coverage compiler flags -if (CMAKE_BUILD_TYPE MATCHES "Debug*" OR CMAKE_BUILD_TYPE MATCHES "Dev*") +if ( ( CMAKE_Fortran_COMPILER MATCHES ".*gfortran.*" OR CMAKE_Fortran_COMPILER MATCHES ".*gcc.*" ) AND + ( CMAKE_BUILD_TYPE MATCHES "Debug*" OR CMAKE_BUILD_TYPE MATCHES "Dev*" ) ) append_coverage_compiler_flags() endif() \ No newline at end of file diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 8ff54a7..3937360 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -2,6 +2,8 @@ This document outlines the organisation of the athena codebase to help guide code contributors. +This document has been copied from the neural-fortran repository and used as a template from which to design this version. The original document can be found here: https://github.com/modern-fortran/neural-fortran/blob/main/CONTRIBUTING.md + ## Overall code organization The source code organisation follows the usual `fpm` convention: diff --git a/README.md b/README.md index 930b60f..789c112 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,28 @@ It was decided that this project should be migrated to allow for better communit --- +## Statement of need + +The ATHENA library leverages Fortran's strong support of array arithmatics, and its compatibility with parallel and high-performance computing resources. +Additionally, there exist many improvements made available since Fortran 95, specifically in Fortran 2018 (Reid 2018) (and upcoming ones in Fortran 2023), as well as continued development by the Fortran Standards committee. +All of this provides a clear incentive to develop further libraries and frameworks focused on providing machine learning capabilities to the Fortran community. + +While existing Fortran-based libraries, such as neural-fortran (Curcic 2019), address many aspects of neural networks, +ATHENA provides implementation of some well-known features not currently available within other libraries; these features include batchnormalisation, regularisation layers (such as dropout and dropblock), and average pooling layers. +Additionally, the library provides support for 1, 2, and 3D input data for most features currently implemented; this includes 1, 2, and 3D data for convolutional layers. +Finally, the ATHENA library supports many convolutional techniques, including various data padding types, and stride. + +One of the primary intended applications of ATHENA is in materials science, which heavily utilises convolutional and graph neural networks for learning based on charge densities and atomic structures. +Given the unique data structure of atomic configurations, specifically their graph-based nature, a specialised API must be developed to accommodate these needs. + +### References +- Reid, J. (2018). The new features of fortran 2018. SIGPLAN Fortran Forum, 37(1), 5–43. https://doi.org/10.1145/3206214.3206215 +- Curcic, M. (2019). A parallel fortran framework for neural networks and deep learning. SIGPLAN Fortran Forum, 38(1), 4–21. https://doi.org/10.1145/3323057.3323059 + + +Documentation +----- + ATHENA is distributed with the following directories: | Directory | Description | @@ -35,9 +57,6 @@ ATHENA is distributed with the following directories: | _test/_ | A set of test programs to check functionality of the library works after compilation | -Documentation ------ - For extended details on the functionality of this library, please check out the [wiki](https://github.com/nedtaylor/athena/wiki) **NOTE: There currently exists no manual document. This will be included at a later date** @@ -63,6 +82,10 @@ The library has been developed and tested using the following compilers: - ifort -- Intel 2021.10.0.20230609 - ifx -- IntelLLVM 2023.2.0 +#### Tested compilers +The library has also been tested with and found to support the following libraries: +- gfortran -- gcc 12.3 + ### Building with fpm The library is set up to work with the Fortran Package Manager (fpm). @@ -76,7 +99,7 @@ Run the following command in the repository main directory: To check whether ATHENA has installed correctly and that the compilation works as expected, the following command can be run: ``` - fpm test + fpm test --profile release ``` This runs a set of test programs (found within the test/ directory) to ensure the expected output occurs when layers and networks are set up. @@ -133,7 +156,7 @@ Using fpm, the examples are built alongside the library. To list all available e To run a particular example, execute the following command: ``` - fpm run --example [NAME] + fpm run --example [NAME] --profile release ``` where [_NAME_] is the name of the example found in the list. diff --git a/example/mnist/src/main.f90 b/example/mnist/src/main.f90 index 906f861..1df53e0 100644 --- a/example/mnist/src/main.f90 +++ b/example/mnist/src/main.f90 @@ -71,7 +71,7 @@ program mnist_example !!!----------------------------------------------------------------------------- !!! initialise random seed !!!----------------------------------------------------------------------------- - call random_setup(seed, num_seed=1, restart=.false.) + call random_setup(seed, restart=.false.) !!!----------------------------------------------------------------------------- diff --git a/example/mnist_3D/src/main.f90 b/example/mnist_3D/src/main.f90 index 64288d1..8de9d7e 100644 --- a/example/mnist_3D/src/main.f90 +++ b/example/mnist_3D/src/main.f90 @@ -72,7 +72,7 @@ program mnist_test !!!----------------------------------------------------------------------------- !!! initialise random seed !!!----------------------------------------------------------------------------- - call random_setup(seed, num_seed=1, restart=.false.) + call random_setup(seed, restart=.false.) !!!----------------------------------------------------------------------------- diff --git a/example/mnist_bn/src/main.f90 b/example/mnist_bn/src/main.f90 index cfb5386..dfb2e49 100644 --- a/example/mnist_bn/src/main.f90 +++ b/example/mnist_bn/src/main.f90 @@ -71,7 +71,7 @@ program mnist_example !!!----------------------------------------------------------------------------- !!! initialise random seed !!!----------------------------------------------------------------------------- - call random_setup(seed, num_seed=1, restart=.false.) + call random_setup(seed, restart=.false.) !!!----------------------------------------------------------------------------- diff --git a/example/mnist_drop/src/main.f90 b/example/mnist_drop/src/main.f90 index 158b055..2603bb5 100644 --- a/example/mnist_drop/src/main.f90 +++ b/example/mnist_drop/src/main.f90 @@ -71,7 +71,7 @@ program mnist_test !!!----------------------------------------------------------------------------- !!! initialise random seed !!!----------------------------------------------------------------------------- - call random_setup(seed, num_seed=1, restart=.false.) + call random_setup(seed, restart=.false.) !!!----------------------------------------------------------------------------- diff --git a/example/simple/src/main.f90 b/example/simple/src/main.f90 index c797f25..1e1e469 100644 --- a/example/simple/src/main.f90 +++ b/example/simple/src/main.f90 @@ -1,3 +1,6 @@ +!! This file contains a modified version of the "simple" example found in ... +!! ... neural fortran: +!! https://github.com/modern-fortran/neural-fortran/blob/main/example/simple.f90 program simple use athena use constants_mnist, only: real12, pi @@ -16,9 +19,8 @@ program simple !! set random seed - seed_size = 8 call random_seed(size=seed_size) - seed = [1,1,1,1,1,1,1,1] + allocate(seed(seed_size), source = 1) call random_seed(put=seed) write(*,*) "Simple function approximation using a fully-connected neural network" diff --git a/example/sine/src/main.f90 b/example/sine/src/main.f90 index fa0997d..dbb6e6b 100644 --- a/example/sine/src/main.f90 +++ b/example/sine/src/main.f90 @@ -1,3 +1,6 @@ +!! This file contains a modified version of the "sine" example found in ... +!! ... neural fortran: +!! https://github.com/modern-fortran/neural-fortran/blob/main/example/sine.f90 program sine use athena use constants_mnist, only: real12, pi diff --git a/fpm.toml b/fpm.toml index da601aa..4645072 100644 --- a/fpm.toml +++ b/fpm.toml @@ -1,5 +1,5 @@ name = "athena" -version = "1.3.0" +version = "1.3.2" license = "MIT" author = "Ned Thaddeus Taylor" maintainer = "n.t.taylor@exeter.ac.uk" diff --git a/src/athena.f90 b/src/athena.f90 index 76a5750..9f06bbd 100644 --- a/src/athena.f90 +++ b/src/athena.f90 @@ -49,6 +49,7 @@ module athena !! input layer types use input1d_layer, only: input1d_layer_type + use input2d_layer, only: input2d_layer_type use input3d_layer, only: input3d_layer_type use input4d_layer, only: input4d_layer_type @@ -58,6 +59,7 @@ module athena use batchnorm3d_layer, only: batchnorm3d_layer_type, read_batchnorm3d_layer !! convolution layer types + use conv1d_layer, only: conv1d_layer_type, read_conv1d_layer use conv2d_layer, only: conv2d_layer_type, read_conv2d_layer use conv3d_layer, only: conv3d_layer_type, read_conv3d_layer @@ -67,8 +69,10 @@ module athena use dropblock3d_layer, only: dropblock3d_layer_type, read_dropblock3d_layer !! pooling layer types + use avgpool1d_layer, only: avgpool1d_layer_type, read_avgpool1d_layer use avgpool2d_layer, only: avgpool2d_layer_type, read_avgpool2d_layer use avgpool3d_layer, only: avgpool3d_layer_type, read_avgpool3d_layer + use maxpool1d_layer, only: maxpool1d_layer_type, read_maxpool1d_layer use maxpool2d_layer, only: maxpool2d_layer_type, read_maxpool2d_layer use maxpool3d_layer, only: maxpool3d_layer_type, read_maxpool3d_layer diff --git a/src/lib/mod_avgpool1d_layer.f90 b/src/lib/mod_avgpool1d_layer.f90 new file mode 100644 index 0000000..0013409 --- /dev/null +++ b/src/lib/mod_avgpool1d_layer.f90 @@ -0,0 +1,463 @@ +!!!############################################################################# +!!! Code written by Ned Thaddeus Taylor +!!! Code part of the ATHENA library - a feedforward neural network library +!!!############################################################################# +!!! module contains implementation of a 1D average pooling layer +!!!############################################################################# +module avgpool1d_layer + use constants, only: real12 + use base_layer, only: pool_layer_type + implicit none + + + type, extends(pool_layer_type) :: avgpool1d_layer_type + real(real12), allocatable, dimension(:,:,:) :: output + real(real12), allocatable, dimension(:,:,:) :: di ! gradient of input (i.e. delta) + contains + procedure, pass(this) :: get_output => get_output_avgpool1d + procedure, pass(this) :: init => init_avgpool1d + procedure, pass(this) :: set_batch_size => set_batch_size_avgpool1d + procedure, pass(this) :: print => print_avgpool1d + procedure, pass(this) :: forward => forward_rank + procedure, pass(this) :: backward => backward_rank + procedure, private, pass(this) :: forward_3d + procedure, private, pass(this) :: backward_3d + end type avgpool1d_layer_type + + + interface avgpool1d_layer_type + module function layer_setup( & + input_shape, batch_size, & + pool_size, stride) result(layer) + integer, dimension(:), optional, intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, dimension(..), optional, intent(in) :: pool_size + integer, dimension(..), optional, intent(in) :: stride + type(avgpool1d_layer_type) :: layer + end function layer_setup + end interface avgpool1d_layer_type + + + private + public :: avgpool1d_layer_type + public :: read_avgpool1d_layer + + +contains + +!!!############################################################################# +!!! get layer outputs +!!!############################################################################# + pure subroutine get_output_avgpool1d(this, output) + implicit none + class(avgpool1d_layer_type), intent(in) :: this + real(real12), allocatable, dimension(..), intent(out) :: output + + select rank(output) + rank(1) + output = reshape(this%output, [size(this%output)]) + rank(2) + output = & + reshape(this%output, [product(this%output_shape),this%batch_size]) + rank(3) + output = this%output + end select + + end subroutine get_output_avgpool1d +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! forward propagation assumed rank handler +!!!############################################################################# + pure subroutine forward_rank(this, input) + implicit none + class(avgpool1d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: input + + select rank(input); rank(3) + call forward_3d(this, input) + end select + end subroutine forward_rank +!!!############################################################################# + + +!!!############################################################################# +!!! backward propagation assumed rank handler +!!!############################################################################# + pure subroutine backward_rank(this, input, gradient) + implicit none + class(avgpool1d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: input + real(real12), dimension(..), intent(in) :: gradient + + select rank(input); rank(3) + select rank(gradient); rank(3) + call backward_3d(this, input, gradient) + end select + end select + end subroutine backward_rank +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! set up layer +!!!############################################################################# +#if defined(GFORTRAN) + module function layer_setup( & + input_shape, batch_size, & + pool_size, stride) result(layer) + implicit none + integer, dimension(:), optional, intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, dimension(..), optional, intent(in) :: pool_size + integer, dimension(..), optional, intent(in) :: stride + + type(avgpool1d_layer_type) :: layer +#else + module procedure layer_setup + implicit none +#endif + + + layer%name = "avgpool1d" + layer%input_rank = 2 + allocate( & + layer%pool(layer%input_rank-1), & + layer%strd(layer%input_rank-1) ) + !!----------------------------------------------------------------------- + !! set up pool size + !!----------------------------------------------------------------------- + if(present(pool_size))then + select rank(pool_size) + rank(0) + layer%pool = pool_size + rank(1) + layer%pool = pool_size + end select + else + layer%pool = 2 + end if + + + !!----------------------------------------------------------------------- + !! set up stride + !!----------------------------------------------------------------------- + if(present(stride))then + select rank(stride) + rank(0) + layer%strd = stride + rank(1) + layer%strd = stride + end select + else + layer%strd = 2 + end if + + + !!-------------------------------------------------------------------------- + !! initialise batch size + !!-------------------------------------------------------------------------- + if(present(batch_size)) layer%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! initialise layer shape + !!-------------------------------------------------------------------------- + if(present(input_shape)) call layer%init(input_shape=input_shape) + +#if defined(GFORTRAN) + end function layer_setup +#else + end procedure layer_setup +#endif +!!!############################################################################# + + +!!!############################################################################# +!!! initialise layer +!!!############################################################################# + subroutine init_avgpool1d(this, input_shape, batch_size, verbose) + implicit none + class(avgpool1d_layer_type), intent(inout) :: this + integer, dimension(:), intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, optional, intent(in) :: verbose + + integer :: verbose_ = 0 + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + if(present(batch_size)) this%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! initialise input shape + !!-------------------------------------------------------------------------- + if(.not.allocated(this%input_shape)) call this%set_shape(input_shape) + + + !!----------------------------------------------------------------------- + !! set up number of channels, width, height + !!----------------------------------------------------------------------- + this%num_channels = this%input_shape(2) + allocate(this%output_shape(2)) + this%output_shape(2) = this%input_shape(2) + this%output_shape(1) = & + floor( (this%input_shape(1) - this%pool(1))/real(this%strd(1))) + 1 + + + !!-------------------------------------------------------------------------- + !! initialise batch size-dependent arrays + !!-------------------------------------------------------------------------- + if(this%batch_size.gt.0) call this%set_batch_size(this%batch_size) + + end subroutine init_avgpool1d +!!!############################################################################# + + +!!!############################################################################# +!!! set batch size +!!!############################################################################# + subroutine set_batch_size_avgpool1d(this, batch_size, verbose) + implicit none + class(avgpool1d_layer_type), intent(inout) :: this + integer, intent(in) :: batch_size + integer, optional, intent(in) :: verbose + + integer :: verbose_ = 0 + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + this%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! allocate arrays + !!-------------------------------------------------------------------------- + if(allocated(this%input_shape))then + if(allocated(this%output)) deallocate(this%output) + allocate(this%output( & + this%output_shape(1), this%num_channels, & + this%batch_size), & + source=0._real12) + if(allocated(this%di)) deallocate(this%di) + allocate(this%di( & + this%input_shape(1), & + this%input_shape(2), & + this%batch_size), & + source=0._real12) + end if + + end subroutine set_batch_size_avgpool1d +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! print layer to file +!!!############################################################################# + subroutine print_avgpool1d(this, file) + implicit none + class(avgpool1d_layer_type), intent(in) :: this + character(*), intent(in) :: file + + integer :: unit + + !! open file with new unit + !!-------------------------------------------------------------------------- + open(newunit=unit, file=trim(file), access='append') + + !! write convolution initial parameters + !!-------------------------------------------------------------------------- + write(unit,'("AVGPOOL1D")') + write(unit,'(3X,"INPUT_SHAPE = ",3(1X,I0))') this%input_shape + write(unit,'(3X,"POOL_SIZE =",1X,I0)') this%pool(1) + write(unit,'(3X,"STRIDE =",1(1X,I0))') this%strd(1) + write(unit,'("END AVGPOOL1D")') + + !! close unit + !!-------------------------------------------------------------------------- + close(unit) + + end subroutine print_avgpool1d +!!!############################################################################# + + +!!!############################################################################# +!!! read layer from file +!!!############################################################################# + function read_avgpool1d_layer(unit) result(layer) + use infile_tools, only: assign_val, assign_vec + use misc, only: to_lower, icount + implicit none + integer, intent(in) :: unit + + class(avgpool1d_layer_type), allocatable :: layer + + integer :: stat + integer :: itmp1 + integer, dimension(1) :: pool_size, stride + integer, dimension(2) :: input_shape + character(256) :: buffer, tag + + + !! loop over tags in layer card + tag_loop: do + + !! check for end of file + read(unit,'(A)',iostat=stat) buffer + if(stat.ne.0)then + write(0,*) "ERROR: file encountered error (EoF?) before END AVGPOOL1D" + stop "Exiting..." + end if + if(trim(adjustl(buffer)).eq."") cycle tag_loop + + !! check for end of convolution card + if(trim(adjustl(buffer)).eq."END AVGPOOL1D")then + backspace(unit) + exit tag_loop + end if + + tag=trim(adjustl(buffer)) + if(scan(buffer,"=").ne.0) tag=trim(tag(:scan(tag,"=")-1)) + + !! read parameters from save file + select case(trim(tag)) + case("INPUT_SHAPE") + call assign_vec(buffer, input_shape, itmp1) + case("POOL_SIZE") + call assign_vec(buffer, pool_size, itmp1) + case("STRIDE") + call assign_vec(buffer, stride, itmp1) + case default + !! don't look for "e" due to scientific notation of numbers + !! ... i.e. exponent (E+00) + if(scan(to_lower(trim(adjustl(buffer))),& + 'abcdfghijklmnopqrstuvwxyz').eq.0)then + cycle tag_loop + elseif(tag(:3).eq.'END')then + cycle tag_loop + end if + stop "Unrecognised line in input file: "//trim(adjustl(buffer)) + end select + end do tag_loop + + !! set transfer activation function + + layer = avgpool1d_layer_type( input_shape=input_shape, & + pool_size = pool_size, stride = stride & + ) + + !! check for end of layer card + read(unit,'(A)') buffer + if(trim(adjustl(buffer)).ne."END AVGPOOL1D")then + write(*,*) trim(adjustl(buffer)) + stop "ERROR: END AVGPOOL1D not where expected" + end if + + end function read_avgpool1d_layer +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! forward propagation +!!!############################################################################# + pure subroutine forward_3d(this, input) + implicit none + class(avgpool1d_layer_type), intent(inout) :: this + real(real12), dimension( & + this%input_shape(1), & + this%num_channels, & + this%batch_size), & + intent(in) :: input + + integer :: i, m, s + integer :: stride_idx + + + this%output = 0._real12 + !! perform the pooling operation + do concurrent(& + s = 1:this%batch_size, & + m = 1:this%num_channels, & + i = 1:this%output_shape(1)) + stride_idx = (i - 1) * this%strd(1) + 1 + this%output(i, m, s) = sum(& + input( & + stride_idx:stride_idx+this%pool(1)-1, m, s)) / & + product(this%pool) + end do + + end subroutine forward_3d +!!!############################################################################# + + +!!!############################################################################# +!!! backward propagation +!!!############################################################################# + pure subroutine backward_3d(this, input, gradient) + implicit none + class(avgpool1d_layer_type), intent(inout) :: this + real(real12), dimension( & + this%input_shape(1), & + this%num_channels, & + this%batch_size), & + intent(in) :: input + real(real12), & + dimension(& + this%output_shape(1), & + this%num_channels, & + this%batch_size), & + intent(in) :: gradient + + integer :: i, m, s + integer :: stride_idx + + + this%di = 0._real12 + !! compute gradients for input feature map + do concurrent( & + s = 1:this%batch_size, & + m = 1:this%num_channels, & + i = 1:this%output_shape(1)) + stride_idx = (i-1) * this%strd(1) + !! compute gradients for input feature map + this%di( & + stride_idx+1:stride_idx+this%pool(1), m, s) = & + this%di( & + stride_idx+1:stride_idx+this%pool(1), m, s) + & + gradient(i, m, s) / product(this%pool) + + end do + + end subroutine backward_3d +!!!############################################################################# + +end module avgpool1d_layer +!!!############################################################################# diff --git a/src/lib/mod_base_layer.f90 b/src/lib/mod_base_layer.f90 index cda77c3..ec14d28 100644 --- a/src/lib/mod_base_layer.f90 +++ b/src/lib/mod_base_layer.f90 @@ -35,6 +35,16 @@ !!! - get_gradients - get the gradients of the layer !!! - set_gradients - set the gradients of the layer !!!############################################################################# +!!! Attribution statement: +!!! The following procedures are based on code from the neural-fortran library +!!! https://github.com/modern-fortran/neural-fortran/blob/main/src/nf/nf_layer.f90 +!!! procedures: +!!! - get_num_params* +!!! - get_params* +!!! - set_params* +!!! - get_gradients* +!!! - set_gradients* +!!!############################################################################# module base_layer use constants, only: real12 use clipper, only: clip_type @@ -106,7 +116,7 @@ end subroutine set_shape_base end interface - abstract interface + interface !!-------------------------------------------------------------------------- !! initialise layer !!-------------------------------------------------------------------------- @@ -134,9 +144,10 @@ module subroutine set_batch_size(this, batch_size, verbose) end subroutine set_batch_size end interface - abstract interface + interface !!-------------------------------------------------------------------------- !! get number of parameters in layer + !! procedure modified from neural-fortran library !!-------------------------------------------------------------------------- !! this = (T, in) layer_type pure module function get_num_params(this) result(num_params) @@ -280,6 +291,7 @@ end subroutine layer_merge !!-------------------------------------------------------------------------- !! get learnable parameters of layer + !! procedure modified from neural-fortran library !!-------------------------------------------------------------------------- !! this = (T, in) layer_type !! param = (R, out) learnable parameters @@ -291,6 +303,7 @@ end function get_params !!-------------------------------------------------------------------------- !! set learnable parameters of layer + !! procedure modified from neural-fortran library !!-------------------------------------------------------------------------- !! this = (T, io) layer_type !! param = (R, in) learnable parameters @@ -302,6 +315,7 @@ end subroutine set_params !!-------------------------------------------------------------------------- !! get parameter gradients of layer + !! procedure modified from neural-fortran library !!-------------------------------------------------------------------------- !! this = (T, in) layer_type !! clip_method = (T, in) clip method @@ -315,6 +329,7 @@ end function get_gradients !!-------------------------------------------------------------------------- !! set learnable parameters of layer + !! procedure modified from neural-fortran library !!-------------------------------------------------------------------------- !! this = (T, io) layer_type !! gradients = (R, in) learnable parameters diff --git a/src/lib/mod_base_layer_sub.f90 b/src/lib/mod_base_layer_sub.f90 index 3de3be5..79b73b6 100644 --- a/src/lib/mod_base_layer_sub.f90 +++ b/src/lib/mod_base_layer_sub.f90 @@ -35,6 +35,16 @@ !!! get_gradients - get the gradients of the layer !!! set_gradients - set the gradients of the layer !!!############################################################################# +!!! Attribution statement: +!!! The following procedures are based on code from the neural-fortran library +!!! https://github.com/modern-fortran/neural-fortran/blob/main/src/nf/nf_layer.f90 +!!! procedures: +!!! - get_num_params* +!!! - get_params* +!!! - set_params* +!!! - get_gradients* +!!! - set_gradients* +!!!############################################################################# submodule(base_layer) base_layer_submodule implicit none @@ -85,6 +95,7 @@ end subroutine set_shape_base !!!############################################################################# !!! get number of parameters in layer +!!! procedure modified from neural-fortran library !!!############################################################################# !!! this = (T, in) layer_type !!! num_params = (I, out) number of parameters @@ -125,6 +136,7 @@ end function get_num_params_batch !!!############################################################################# !!! get learnable parameters of layer +!!! procedure modified from neural-fortran library !!!############################################################################# pure module function get_params_batch(this) result(params) implicit none @@ -139,6 +151,7 @@ end function get_params_batch !!!############################################################################# !!! set learnable parameters of layer +!!! procedure modified from neural-fortran library !!!############################################################################# module subroutine set_params_batch(this, params) implicit none @@ -154,6 +167,7 @@ end subroutine set_params_batch !!!############################################################################# !!! get gradients of layer +!!! procedure modified from neural-fortran library !!!############################################################################# pure module function get_gradients_batch(this, clip_method) result(gradients) use clipper, only: clip_type @@ -172,6 +186,7 @@ end function get_gradients_batch !!!############################################################################# !!! set gradients of layer +!!! procedure modified from neural-fortran library !!!############################################################################# module subroutine set_gradients_batch(this, gradients) implicit none diff --git a/src/lib/mod_container_layer_sub.f90 b/src/lib/mod_container_layer_sub.f90 index d25c761..4015ff7 100644 --- a/src/lib/mod_container_layer_sub.f90 +++ b/src/lib/mod_container_layer_sub.f90 @@ -8,18 +8,22 @@ submodule(container_layer) container_layer_submodule use base_layer, only: learnable_layer_type, flatten_layer_type use input1d_layer, only: input1d_layer_type + use input2d_layer, only: input2d_layer_type use input3d_layer, only: input3d_layer_type use input4d_layer, only: input4d_layer_type use batchnorm1d_layer, only: batchnorm1d_layer_type use batchnorm2d_layer, only: batchnorm2d_layer_type use batchnorm3d_layer, only: batchnorm3d_layer_type + use conv1d_layer, only: conv1d_layer_type use conv2d_layer, only: conv2d_layer_type use conv3d_layer, only: conv3d_layer_type use dropout_layer, only: dropout_layer_type use dropblock2d_layer, only: dropblock2d_layer_type use dropblock3d_layer, only: dropblock3d_layer_type + use avgpool1d_layer, only: avgpool1d_layer_type use avgpool2d_layer, only: avgpool2d_layer_type use avgpool3d_layer, only: avgpool3d_layer_type + use maxpool1d_layer, only: maxpool1d_layer_type use maxpool2d_layer, only: maxpool2d_layer_type use maxpool3d_layer, only: maxpool3d_layer_type use full_layer, only: full_layer_type @@ -34,6 +38,8 @@ pure module subroutine forward(this, input) select type(previous => input%layer) type is(input1d_layer_type) call this%layer%forward(previous%output) + type is(input2d_layer_type) + call this%layer%forward(previous%output) type is(input3d_layer_type) call this%layer%forward(previous%output) type is(input4d_layer_type) @@ -46,6 +52,8 @@ pure module subroutine forward(this, input) type is(batchnorm3d_layer_type) call this%layer%forward(previous%output) + type is(conv1d_layer_type) + call this%layer%forward(previous%output) type is(conv2d_layer_type) call this%layer%forward(previous%output) type is(conv3d_layer_type) @@ -58,10 +66,14 @@ pure module subroutine forward(this, input) type is(dropblock3d_layer_type) call this%layer%forward(previous%output) + type is(avgpool1d_layer_type) + call this%layer%forward(previous%output) type is(avgpool2d_layer_type) call this%layer%forward(previous%output) type is(avgpool3d_layer_type) call this%layer%forward(previous%output) + type is(maxpool1d_layer_type) + call this%layer%forward(previous%output) type is(maxpool2d_layer_type) call this%layer%forward(previous%output) type is(maxpool3d_layer_type) @@ -86,6 +98,8 @@ pure module subroutine backward(this, input, gradient) select type(previous => input%layer) type is(input1d_layer_type) call this%layer%backward(previous%output, gradient) + type is(input2d_layer_type) + call this%layer%backward(previous%output, gradient) type is(input3d_layer_type) call this%layer%backward(previous%output, gradient) type is(input4d_layer_type) @@ -98,6 +112,8 @@ pure module subroutine backward(this, input, gradient) type is(batchnorm3d_layer_type) call this%layer%backward(previous%output, gradient) + type is(conv1d_layer_type) + call this%layer%backward(previous%output, gradient) type is(conv2d_layer_type) call this%layer%backward(previous%output, gradient) type is(conv3d_layer_type) @@ -110,10 +126,14 @@ pure module subroutine backward(this, input, gradient) type is(dropblock3d_layer_type) call this%layer%backward(previous%output, gradient) + type is(avgpool1d_layer_type) + call this%layer%backward(previous%output, gradient) type is(avgpool2d_layer_type) call this%layer%backward(previous%output, gradient) type is(avgpool3d_layer_type) call this%layer%backward(previous%output, gradient) + type is(maxpool1d_layer_type) + call this%layer%backward(previous%output, gradient) type is(maxpool2d_layer_type) call this%layer%backward(previous%output, gradient) type is(maxpool3d_layer_type) diff --git a/src/lib/mod_conv1d_layer.f90 b/src/lib/mod_conv1d_layer.f90 new file mode 100644 index 0000000..842d6b6 --- /dev/null +++ b/src/lib/mod_conv1d_layer.f90 @@ -0,0 +1,960 @@ +!!!############################################################################# +!!! Code written by Ned Thaddeus Taylor +!!! Code part of the ATHENA library - a feedforward neural network library +!!!############################################################################# +!!! module contains implementation of a 1D convolutional layer +!!!############################################################################# +!!! Attribution statement: +!!! The following procedures are based on code from the neural-fortran library +!!! https://github.com/modern-fortran/neural-fortran/blob/main/src/nf/nf_layer.f90 +!!! procedures: +!!! - get_params* +!!! - set_params* +!!! - get_gradients* +!!! - set_gradients* +!!!############################################################################# +module conv1d_layer + use constants, only: real12 + use base_layer, only: learnable_layer_type, conv_layer_type + use custom_types, only: initialiser_type + implicit none + + + type, extends(conv_layer_type) :: conv1d_layer_type + real(real12), allocatable, dimension(:,:,:) :: weight + real(real12), allocatable, dimension(:,:,:,:) :: dw ! weight gradient + real(real12), allocatable, dimension(:,:,:) :: output, z + real(real12), allocatable, dimension(:,:,:) :: di ! input gradient + contains + procedure, pass(this) :: get_params => get_params_conv1d + procedure, pass(this) :: set_params => set_params_conv1d + procedure, pass(this) :: get_gradients => get_gradients_conv1d + procedure, pass(this) :: set_gradients => set_gradients_conv1d + procedure, pass(this) :: get_output => get_output_conv1d + + procedure, pass(this) :: init => init_conv1d + procedure, pass(this) :: set_batch_size => set_batch_size_conv1d + procedure, pass(this) :: print => print_conv1d + + procedure, pass(this) :: forward => forward_rank + procedure, pass(this) :: backward => backward_rank + procedure, private, pass(this) :: forward_3d + procedure, private, pass(this) :: backward_3d + + procedure, pass(this) :: reduce => layer_reduction + procedure, pass(this) :: merge => layer_merge + procedure :: add_t_t => layer_add !t = type, r = real, i = int + generic :: operator(+) => add_t_t !, public + end type conv1d_layer_type + + +!!!----------------------------------------------------------------------------- +!!! interface for layer set up +!!!----------------------------------------------------------------------------- + interface conv1d_layer_type + module function layer_setup( & + input_shape, batch_size, & + num_filters, kernel_size, stride, padding, & + activation_function, activation_scale, & + kernel_initialiser, bias_initialiser, & + calc_input_gradients) result(layer) + integer, dimension(:), optional, intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, optional, intent(in) :: num_filters + integer, dimension(..), optional, intent(in) :: kernel_size + integer, dimension(..), optional, intent(in) :: stride + real(real12), optional, intent(in) :: activation_scale + character(*), optional, intent(in) :: activation_function, & + kernel_initialiser, bias_initialiser, padding + logical, optional, intent(in) :: calc_input_gradients + type(conv1d_layer_type) :: layer + end function layer_setup + end interface conv1d_layer_type + + + private + public :: conv1d_layer_type + public :: read_conv1d_layer + + +contains + +!!!############################################################################# +!!! layer reduction +!!!############################################################################# + subroutine layer_reduction(this, rhs) + implicit none + class(conv1d_layer_type), intent(inout) :: this + class(learnable_layer_type), intent(in) :: rhs + + select type(rhs) + class is(conv1d_layer_type) + this%db = this%db + rhs%db + this%dw = this%dw + rhs%dw + end select + + end subroutine layer_reduction +!!!############################################################################# + + +!!!############################################################################# +!!! layer addition +!!!############################################################################# + function layer_add(a, b) result(output) + implicit none + class(conv1d_layer_type), intent(in) :: a, b + type(conv1d_layer_type) :: output + + output = a + output%dw = output%dw + b%dw + output%db = output%db + b%db + + end function layer_add +!!!############################################################################# + + +!!!############################################################################# +!!! layer merge +!!!############################################################################# + subroutine layer_merge(this, input) + implicit none + class(conv1d_layer_type), intent(inout) :: this + class(learnable_layer_type), intent(in) :: input + + select type(input) + class is(conv1d_layer_type) + this%dw = this%dw + input%dw + this%db = this%db + input%db + end select + + end subroutine layer_merge +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! get learnable parameters +!!! procedure modified from neural-fortran library +!!!############################################################################# + pure function get_params_conv1d(this) result(params) + implicit none + class(conv1d_layer_type), intent(in) :: this + real(real12), allocatable, dimension(:) :: params + + params = [ reshape( & + this%weight, & + [ this%num_filters * this%num_channels * product(this%knl) ]), & + this%bias ] + + end function get_params_conv1d +!!!############################################################################# + + +!!!############################################################################# +!!! set learnable parameters +!!! procedure modified from neural-fortran library +!!!############################################################################# + subroutine set_params_conv1d(this, params) + implicit none + class(conv1d_layer_type), intent(inout) :: this + real(real12), dimension(:), intent(in) :: params + + this%weight = reshape( & + params(1:this%num_filters * this%num_channels * product(this%knl)), & + shape(this%weight)) + this%bias = params(& + this%num_filters * this%num_channels * product(this%knl) + 1 : ) + + end subroutine set_params_conv1d +!!!############################################################################# + + +!!!############################################################################# +!!! get sample-average gradients +!!! sum over batch dimension and divide by batch size +!!! procedure modified from neural-fortran library +!!!############################################################################# + pure function get_gradients_conv1d(this, clip_method) result(gradients) + use clipper, only: clip_type + implicit none + class(conv1d_layer_type), intent(in) :: this + type(clip_type), optional, intent(in) :: clip_method + real(real12), allocatable, dimension(:) :: gradients + + gradients = [ reshape( & + sum(this%dw,dim=4)/this%batch_size, & + [ this%num_filters * this%num_channels * product(this%knl) ]), & + sum(this%db,dim=2)/this%batch_size ] + + if(present(clip_method)) call clip_method%apply(size(gradients),gradients) + + end function get_gradients_conv1d +!!!############################################################################# + + +!!!############################################################################# +!!! set gradients +!!! procedure modified from neural-fortran library +!!!############################################################################# + subroutine set_gradients_conv1d(this, gradients) + implicit none + class(conv1d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: gradients + + integer :: s + + select rank(gradients) + rank(0) + this%dw = gradients + this%db = gradients + rank(1) + do s=1,this%batch_size + this%dw(:,:,:,s) = reshape(gradients(:& + this%num_filters * this%num_channels * product(this%knl)), & + shape(this%dw(:,:,:,s))) + this%db(:,s) = gradients(& + this%num_filters * this%num_channels * product(this%knl)+1:) + end do + end select + + end subroutine set_gradients_conv1d +!!!############################################################################# + + +!!!############################################################################# +!!! get layer outputs +!!!############################################################################# + pure subroutine get_output_conv1d(this, output) + implicit none + class(conv1d_layer_type), intent(in) :: this + real(real12), allocatable, dimension(..), intent(out) :: output + + select rank(output) + rank(1) + output = reshape(this%output, [size(this%output)]) + rank(2) + output = & + reshape(this%output, [product(this%output_shape),this%batch_size]) + rank(3) + output = this%output + end select + + end subroutine get_output_conv1d +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! forward propagation assumed rank handler +!!!############################################################################# + pure subroutine forward_rank(this, input) + implicit none + class(conv1d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: input + + select rank(input); rank(3) + call forward_3d(this, input) + end select + end subroutine forward_rank +!!!############################################################################# + + +!!!############################################################################# +!!! backward propagation assumed rank handler +!!!############################################################################# + pure subroutine backward_rank(this, input, gradient) + implicit none + class(conv1d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: input + real(real12), dimension(..), intent(in) :: gradient + + select rank(input); rank(3) + select rank(gradient) + rank(1) + call backward_3d(this, input, gradient) + rank(2) + call backward_3d(this, input, gradient) + rank(3) + call backward_3d(this, input, gradient) + end select + end select + end subroutine backward_rank +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! set up layer +!!!############################################################################# + module function layer_setup( & + input_shape, batch_size, & + num_filters, kernel_size, stride, padding, & + activation_function, activation_scale, & + kernel_initialiser, bias_initialiser, & + calc_input_gradients) result(layer) + !! add in dilation + use activation, only: activation_setup + use initialiser, only: get_default_initialiser + use misc_ml, only: set_padding + implicit none + integer, dimension(:), optional, intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, optional, intent(in) :: num_filters + integer, dimension(..), optional, intent(in) :: kernel_size + integer, dimension(..), optional, intent(in) :: stride + real(real12), optional, intent(in) :: activation_scale + character(*), optional, intent(in) :: activation_function, & + kernel_initialiser, bias_initialiser, padding + logical, optional, intent(in) :: calc_input_gradients + + type(conv1d_layer_type) :: layer + + integer :: i + real(real12) :: scale + character(len=10) :: activation_function_ + character(len=20) :: padding_ + + + layer%name = "conv1d" + layer%input_rank = 2 + allocate( & + layer%knl(layer%input_rank-1), & + layer%stp(layer%input_rank-1), & + layer%hlf(layer%input_rank-1), & + layer%pad(layer%input_rank-1), & + layer%cen(layer%input_rank-1) ) + !!-------------------------------------------------------------------------- + !! initialise batch size + !!-------------------------------------------------------------------------- + if(present(batch_size)) layer%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! determine whether to calculate input gradients + !!-------------------------------------------------------------------------- + if(present(calc_input_gradients))then + layer%calc_input_gradients = calc_input_gradients + write(*,*) "CONV1D input gradients turned off" + else + layer%calc_input_gradients = .true. + end if + + + !!-------------------------------------------------------------------------- + !! set up number of filters + !!-------------------------------------------------------------------------- + if(present(num_filters))then + layer%num_filters = num_filters + else + layer%num_filters = 32 + end if + + + !!-------------------------------------------------------------------------- + !! set up kernel size + !!-------------------------------------------------------------------------- + if(present(kernel_size))then + select rank(kernel_size) + rank(0) + layer%knl = kernel_size + rank(1) + layer%knl(1) = kernel_size(1) + end select + else + layer%knl = 3 + end if + !! odd or even kernel/filter size + !!-------------------------------------------------------------------------- + layer%cen = 2 - mod(layer%knl, 2) + layer%hlf = (layer%knl-1)/2 + + if(present(padding))then + padding_ = padding + else + padding_ = "valid" + end if + call set_padding(layer%pad(1), layer%knl(1), padding_) + + + !!-------------------------------------------------------------------------- + !! set up stride + !!-------------------------------------------------------------------------- + if(present(stride))then + select rank(stride) + rank(0) + layer%stp = stride + rank(1) + layer%stp(1) = stride(1) + end select + else + layer%stp = 1 + end if + + + !!-------------------------------------------------------------------------- + !! set activation and derivative functions based on input name + !!-------------------------------------------------------------------------- + if(present(activation_function))then + activation_function_ = activation_function + else + activation_function_ = "none" + end if + if(present(activation_scale))then + scale = activation_scale + else + scale = 1._real12 + end if + write(*,'("CONV1D activation function: ",A)') trim(activation_function_) + allocate(layer%transfer, & + source=activation_setup(activation_function_, scale)) + + + !!-------------------------------------------------------------------------- + !! define weights (kernels) and biases initialisers + !!-------------------------------------------------------------------------- + if(present(kernel_initialiser)) layer%kernel_initialiser =kernel_initialiser + if(trim(layer%kernel_initialiser).eq.'') & + layer%kernel_initialiser=get_default_initialiser(activation_function_) + write(*,'("CONV1D kernel initialiser: ",A)') trim(layer%kernel_initialiser) + if(present(bias_initialiser)) layer%bias_initialiser = bias_initialiser + if(trim(layer%bias_initialiser).eq.'') & + layer%bias_initialiser = get_default_initialiser(& + activation_function_, is_bias=.true.) + write(*,'("CONV1D bias initialiser: ",A)') trim(layer%bias_initialiser) + + + !!-------------------------------------------------------------------------- + !! initialise layer shape + !!-------------------------------------------------------------------------- + if(present(input_shape)) call layer%init(input_shape=input_shape) + + end function layer_setup +!!!############################################################################# + + +!!!############################################################################# +!!! initialise layer +!!!############################################################################# + subroutine init_conv1d(this, input_shape, batch_size, verbose) + use initialiser, only: initialiser_setup + implicit none + class(conv1d_layer_type), intent(inout) :: this + integer, dimension(:), intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, optional, intent(in) :: verbose + + integer :: verbose_ = 0 + integer :: end_idx + class(initialiser_type), allocatable :: initialiser_ + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + if(present(batch_size)) this%batch_size = batch_size + + + !!------------------------------------------------------------------------- + !! initialise input shape + !!-------------------------------------------------------------------------- + if(.not.allocated(this%input_shape)) call this%set_shape(input_shape) + + + !!-------------------------------------------------------------------------- + !! allocate output, activation, bias, and weight shapes + !!-------------------------------------------------------------------------- + !! NOTE: INPUT SHAPE DOES NOT INCLUDE PADDING WIDTH + !! THIS IS HANDLED AUTOMATICALLY BY THE CODE + !! ... provide the initial input data shape and let us deal with the padding + this%num_channels = this%input_shape(2) + allocate(this%output_shape(2)) + this%output_shape(2) = this%num_filters + this%output_shape(1) = floor(& + (this%input_shape(1) + & + 2.0 * this%pad(1) - this%knl(1))/real(this%stp(1)) ) + 1 + + allocate(this%bias(this%num_filters), source=0._real12) + + end_idx = this%hlf(1) + (this%cen(1) - 1) + allocate(this%weight( & + -this%hlf(1):end_idx, & + this%num_channels,this%num_filters), source=0._real12) + + + !!-------------------------------------------------------------------------- + !! initialise weights (kernels) + !!-------------------------------------------------------------------------- + allocate(initialiser_, source=initialiser_setup(this%kernel_initialiser)) + call initialiser_%initialise(this%weight, & + fan_in=product(this%knl)+1, fan_out=1) + deallocate(initialiser_) + + !! initialise biases + !!-------------------------------------------------------------------------- + allocate(initialiser_, source=initialiser_setup(this%bias_initialiser)) + call initialiser_%initialise(this%bias, & + fan_in=product(this%knl)+1, fan_out=1) + deallocate(initialiser_) + + + !!-------------------------------------------------------------------------- + !! initialise batch size-dependent arrays + !!-------------------------------------------------------------------------- + if(this%batch_size.gt.0) call this%set_batch_size(this%batch_size) + + end subroutine init_conv1d +!!!############################################################################# + + +!!!############################################################################# +!!! set batch size +!!!############################################################################# + subroutine set_batch_size_conv1d(this, batch_size, verbose) + implicit none + class(conv1d_layer_type), intent(inout) :: this + integer, intent(in) :: batch_size + integer, optional, intent(in) :: verbose + + integer :: verbose_ = 0 + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + this%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! allocate arrays + !!-------------------------------------------------------------------------- + if(allocated(this%input_shape))then + if(allocated(this%output)) deallocate(this%output) + allocate(this%output( & + this%output_shape(1), & + this%num_filters, & + this%batch_size), source=0._real12) + if(allocated(this%z)) deallocate(this%z) + allocate(this%z, source=this%output) + if(allocated(this%di)) deallocate(this%di) + allocate(this%di( & + this%input_shape(1), & + this%input_shape(2), & + this%batch_size), source=0._real12) + if(allocated(this%dw)) deallocate(this%dw) + allocate(this%dw( & + lbound(this%weight,1):ubound(this%weight,1), & + this%num_channels, this%num_filters, & + this%batch_size), source=0._real12) + if(allocated(this%db)) deallocate(this%db) + allocate(this%db(this%num_filters, this%batch_size), source=0._real12) + end if + + end subroutine set_batch_size_conv1d +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! print layer to file +!!!############################################################################# + subroutine print_conv1d(this, file) + implicit none + class(conv1d_layer_type), intent(in) :: this + character(*), intent(in) :: file + + integer :: l, i + integer :: unit + character(:), allocatable :: padding_type + + + !! determine padding method + !!-------------------------------------------------------------------------- + padding_type = "" + if(this%pad(1).eq.this%knl(1)-1)then + padding_type = "full" + elseif(this%pad(1).eq.0)then + padding_type = "valid" + else + padding_type = "same" + end if + + !! open file with new unit + !!-------------------------------------------------------------------------- + open(newunit=unit, file=trim(file), access='append') + + !! write convolution initial parameters + !!-------------------------------------------------------------------------- + write(unit,'("CONV1D")') + write(unit,'(3X,"INPUT_SHAPE = ",3(1X,I0))') this%input_shape + write(unit,'(3X,"NUM_FILTERS = ",I0)') this%num_filters + write(unit,'(3X,"KERNEL_SIZE =",1X,I0)') this%knl(1) + write(unit,'(3X,"STRIDE =",1X,I0)') this%stp(1) + write(unit,'(3X,"PADDING = ",A)') padding_type + + write(unit,'(3X,"ACTIVATION = ",A)') trim(this%transfer%name) + write(unit,'(3X,"ACTIVATION_SCALE = ",F0.9)') this%transfer%scale + + !! write convolution weights and biases + !!-------------------------------------------------------------------------- + write(unit,'("WEIGHTS")') + do l=1,this%num_filters + write(unit,'(5(E16.8E2))', advance="no") this%weight(:,:,l) + if(mod(size(this%weight(:,:,l)),5).eq.0) write(unit,*) + write(unit,'(E16.8E2)') this%bias(l) + end do + write(unit,'("END WEIGHTS")') + write(unit,'("END CONV1D")') + + !! close unit + !!-------------------------------------------------------------------------- + close(unit) + + end subroutine print_conv1d +!!!############################################################################# + + +!!!############################################################################# +!!! read layer from file +!!!############################################################################# + function read_conv1d_layer(unit, verbose) result(layer) + use infile_tools, only: assign_val, assign_vec + use misc, only: to_lower, icount + implicit none + integer, intent(in) :: unit + integer, optional, intent(in) :: verbose + + class(conv1d_layer_type), allocatable :: layer + + integer :: stat, verbose_ = 0 + integer :: j, k, l, c, itmp1 + integer :: num_filters, num_inputs + real(real12) :: activation_scale + logical :: found_weights = .false. + character(14) :: kernel_initialiser='', bias_initialiser='' + character(20) :: padding, activation_function + character(256) :: buffer, tag + + integer, dimension(1) :: kernel_size, stride + integer, dimension(2) :: input_shape + real(real12), allocatable, dimension(:) :: data_list + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + + + !!-------------------------------------------------------------------------- + !! loop over tags in layer card + !!-------------------------------------------------------------------------- + tag_loop: do + + !! check for end of file + !!----------------------------------------------------------------------- + read(unit,'(A)',iostat=stat) buffer + if(stat.ne.0)then + write(0,*) "ERROR: file encountered error (EoF?) before END CONV1D" + stop "Exiting..." + end if + if(trim(adjustl(buffer)).eq."") cycle tag_loop + + !! check for end of layer card + !!----------------------------------------------------------------------- + if(trim(adjustl(buffer)).eq."END CONV1D")then + backspace(unit) + exit tag_loop + end if + + tag=trim(adjustl(buffer)) + if(scan(buffer,"=").ne.0) tag=trim(tag(:scan(tag,"=")-1)) + + !! read parameters from save file + !!----------------------------------------------------------------------- + select case(trim(tag)) + case("INPUT_SHAPE") + call assign_vec(buffer, input_shape, itmp1) + case("NUM_FILTERS") + call assign_val(buffer, num_filters, itmp1) + case("KERNEL_SIZE") + call assign_vec(buffer, kernel_size, itmp1) + case("STRIDE") + call assign_vec(buffer, stride, itmp1) + case("PADDING") + call assign_val(buffer, padding, itmp1) + padding = to_lower(padding) + case("ACTIVATION") + call assign_val(buffer, activation_function, itmp1) + case("ACTIVATION_SCALE") + call assign_val(buffer, activation_scale, itmp1) + case("KERNEL_INITIALISER") + call assign_val(buffer, kernel_initialiser, itmp1) + case("BIAS_INITIALISER") + call assign_val(buffer, bias_initialiser, itmp1) + case("WEIGHTS") + found_weights = .true. + kernel_initialiser = 'zeros' + bias_initialiser = 'zeros' + exit tag_loop + case default + !! don't look for "e" due to scientific notation of numbers + !! ... i.e. exponent (E+00) + if(scan(to_lower(trim(adjustl(buffer))),& + 'abcdfghijklmnopqrstuvwxyz').eq.0)then + cycle tag_loop + elseif(tag(:3).eq.'END')then + cycle tag_loop + end if + stop "Unrecognised line in input file: "//trim(adjustl(buffer)) + end select + end do tag_loop + + + !!-------------------------------------------------------------------------- + !! allocate layer + !!-------------------------------------------------------------------------- + layer = conv1d_layer_type( & + input_shape = input_shape, & + num_filters = num_filters, & + kernel_size = kernel_size, stride = stride, & + padding = padding, & + activation_function = activation_function, & + activation_scale = activation_scale, & + kernel_initialiser = kernel_initialiser, & + bias_initialiser = bias_initialiser) + + + !!-------------------------------------------------------------------------- + !! check if WEIGHTS card was found + !!-------------------------------------------------------------------------- + if(.not.found_weights)then + write(0,*) "WARNING: WEIGHTS card in CONV1D not found" + else + do l=1,num_filters + num_inputs = product(layer%knl) + 1 !+1 for bias + allocate(data_list(num_inputs), source=0._real12) + c = 1 + k = 1 + data_concat_loop: do while(c.le.num_inputs) + read(unit,'(A)',iostat=stat) buffer + if(stat.ne.0) exit data_concat_loop + k = icount(buffer) + read(buffer,*,iostat=stat) (data_list(j),j=c,c+k-1) + c = c + k + end do data_concat_loop + layer%weight(:,:,l) = & + reshape(& + data_list(1:num_inputs-1),& + shape(layer%weight(:,:,l))) + layer%bias(l) = data_list(num_inputs) + deallocate(data_list) + end do + + !! check for end of weights card + !!----------------------------------------------------------------------- + read(unit,'(A)') buffer + if(trim(adjustl(buffer)).ne."END WEIGHTS")then + write(*,*) trim(adjustl(buffer)) + stop "ERROR: END WEIGHTS not where expected" + end if + end if + + + !!-------------------------------------------------------------------------- + !! check for end of layer card + !!-------------------------------------------------------------------------- + read(unit,'(A)') buffer + if(trim(adjustl(buffer)).ne."END CONV1D")then + write(*,*) trim(adjustl(buffer)) + stop "ERROR: END CONV1D not where expected" + end if + + end function read_conv1d_layer +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! forward propagation +!!!############################################################################# + pure subroutine forward_3d(this, input) + implicit none + class(conv1d_layer_type), intent(inout) :: this + real(real12), & + dimension( & + -this%pad(1)+1:this%input_shape(1)+this%pad(1), & + this%num_channels,this%batch_size), & + intent(in) :: input + + integer :: i, l, s + integer :: stp_idx, start_idx, end_idx + + + !! perform the convolution operation + !!-------------------------------------------------------------------------- + do concurrent(i=1:this%output_shape(1):1) + stp_idx = (i-1)*this%stp(1) + 1 + (this%hlf(1) - this%pad(1)) + start_idx = stp_idx - this%hlf(1) + end_idx = start_idx + this%knl(1) - 1 + + do concurrent(s=1:this%batch_size) + this%z(i,:,s) = this%bias(:) + end do + + do concurrent(l=1:this%num_filters, s=1:this%batch_size) + this%z(i,l,s) = this%z(i,l,s) + & + sum( & + input( & + start_idx:end_idx,:,s) * & + this%weight(:,:,l) & + ) + end do + end do + + + !! apply activation function to activation values (z) + !!-------------------------------------------------------------------------- + this%output = this%transfer%activate(this%z) + + end subroutine forward_3d +!!!############################################################################# + + +!!!############################################################################# +!!! backward propagation +!!! method : gradient descent +!!!############################################################################# + pure subroutine backward_3d(this, input, gradient) + implicit none + class(conv1d_layer_type), intent(inout) :: this + real(real12), & + dimension( & + -this%pad(1)+1:this%input_shape(1)+this%pad(1), & + this%num_channels,this%batch_size), & + intent(in) :: input + real(real12), & + dimension( & + this%output_shape(1), & + this%num_filters,this%batch_size), & + intent(in) :: gradient + + integer :: l, m, i, x, s + integer :: stp_idx, offset, adjust, end_idx, n_stp + integer, dimension(2) :: lim, lim_w, lim_g + real(real12), & + dimension( & + this%output_shape(1),this%num_filters, & + this%batch_size) :: grad_dz + + + real(real12), dimension(1) :: bias_diff + bias_diff = this%transfer%differentiate([1._real12]) + + + !! get size of the input and output feature maps + !!-------------------------------------------------------------------------- + end_idx = this%hlf(1) + (this%cen(1) - 1) + offset = 1 + this%hlf(1) - this%pad(1) + adjust = 2 * max(this%pad(1), this%hlf(1)) + + + !! get gradient multiplied by differential of Z + !!-------------------------------------------------------------------------- + grad_dz = 0._real12 + grad_dz(& + 1:this%output_shape(1),:,:) = gradient * & + this%transfer%differentiate(this%z) + do concurrent(l=1:this%num_filters, s=1:this%batch_size) + this%db(l,s) = this%db(l,s) + sum(grad_dz(:,l,s)) * bias_diff(1) + end do + + !! apply convolution to compute weight gradients + !! offset applied as centre of kernel is 0 ... + !! ... whilst the starting index for input is 1 + !!-------------------------------------------------------------------------- + do concurrent( & + s=1:this%batch_size, & + l=1:this%num_filters, & + m=1:this%num_channels, & + x=-this%hlf(1):end_idx:1 & + ) + this%dw(x,m,l,s) = this%dw(x,m,l,s) + & + sum(grad_dz(:,l,s) * & + input( & + x+offset:x+offset-1+size(input,1)-adjust:this%stp(1),m,s)) + end do + + + !! apply strided convolution to obtain input gradients + !!-------------------------------------------------------------------------- + if(this%calc_input_gradients)then + lim(1) = this%knl(1) - 1 + lim(2) = (this%output_shape(1) - 1) * this%stp(1) + 1 + end_idx + n_stp = this%output_shape(1) * this%stp(1) + this%di = 0._real12 + !! all elements of the output are separated by stride_x (stride_y) + do concurrent( & + s=1:this%batch_size, & + l=1:this%num_filters, & + m=1:this%num_channels, & + i=1:size(this%di,dim=1):1 & + ) + + !! set weight bounds + stp_idx = ( i - offset )/this%stp(1) + 1 + !! max( ... + !! ... 1. offset of 1st o/p idx from centre of knl (lim) + !! ... 2. lwst o/p idx overlap with <<- knl idx (rpt. pattern) + !! ...) + lim_w(2) = max(lim(1)-i, -this%hlf(1) + & + mod(n_stp+this%knl(1)-i,this%stp(1))) + !! min( ... + !! ... 1. offset of last o/p idx from centre of knl (lim) + !! ... 2. hghst o/p idx overlap with ->> knl idx (rpt. pattern) + !! ...) + lim_w(1) = min(lim(2)-i, end_idx - mod(n_stp-1+i,this%stp(1))) + if(lim_w(2).gt.lim_w(1)) cycle + + !! set gradient bounds + lim_g(1) = max(1, i - offset) + lim_g(2) = min(this%output_shape(1), i - offset + this%knl(1) - 1) + + !! apply full convolution to compute input gradients + this%di(i,m,s) = & + this%di(i,m,s) + & + sum( & + grad_dz( & + lim_g(1):lim_g(2),l,s) * & + this%weight(& + lim_w(1):lim_w(2):-this%stp(1),m,l) ) + + end do + end if + + end subroutine backward_3d +!!!############################################################################# + +end module conv1d_layer +!!!############################################################################# diff --git a/src/lib/mod_conv2d_layer.f90 b/src/lib/mod_conv2d_layer.f90 index f4480f4..8143879 100644 --- a/src/lib/mod_conv2d_layer.f90 +++ b/src/lib/mod_conv2d_layer.f90 @@ -4,6 +4,15 @@ !!!############################################################################# !!! module contains implementation of a 2D convolutional layer !!!############################################################################# +!!! Attribution statement: +!!! The following procedures are based on code from the neural-fortran library +!!! https://github.com/modern-fortran/neural-fortran/blob/main/src/nf/nf_layer.f90 +!!! procedures: +!!! - get_params* +!!! - set_params* +!!! - get_gradients* +!!! - set_gradients* +!!!############################################################################# module conv2d_layer use constants, only: real12 use base_layer, only: learnable_layer_type, conv_layer_type @@ -129,6 +138,7 @@ end subroutine layer_merge !!!############################################################################# !!! get learnable parameters +!!! procedure modified from neural-fortran library !!!############################################################################# pure function get_params_conv2d(this) result(params) implicit none @@ -146,6 +156,7 @@ end function get_params_conv2d !!!############################################################################# !!! set learnable parameters +!!! procedure modified from neural-fortran library !!!############################################################################# subroutine set_params_conv2d(this, params) implicit none @@ -165,6 +176,7 @@ end subroutine set_params_conv2d !!!############################################################################# !!! get sample-average gradients !!! sum over batch dimension and divide by batch size +!!! procedure modified from neural-fortran library !!!############################################################################# pure function get_gradients_conv2d(this, clip_method) result(gradients) use clipper, only: clip_type @@ -186,6 +198,7 @@ end function get_gradients_conv2d !!!############################################################################# !!! set gradients +!!! procedure modified from neural-fortran library !!!############################################################################# subroutine set_gradients_conv2d(this, gradients) implicit none @@ -970,19 +983,37 @@ pure subroutine backward_4d(this, input, gradient) !! ... 1. offset of 1st o/p idx from centre of knl (lim) !! ... 2. lwst o/p idx overlap with <<- knl idx (rpt. pattern) !! ...) - lim_w(2,:) = max(lim(1,:)-[i,j], -this%hlf + & - mod(n_stp+this%knl-[i,j],this%stp)) !! min( ... !! ... 1. offset of last o/p idx from centre of knl (lim) !! ... 2. hghst o/p idx overlap with ->> knl idx (rpt. pattern) !! ...) +#if defined(GFORTRAN) + lim_w(2,:) = max(lim(1,:)-[i,j], -this%hlf + & + mod(n_stp+this%knl-[i,j],this%stp)) lim_w(1,:) = min(lim(2,:)-[i,j], end_idx - & mod(n_stp-1+[i,j],this%stp)) +#else + lim_w(2,1) = max(lim(1,1) - i, -this%hlf(1) + & + mod(n_stp(1)+this%knl(1)-i,this%stp(1))) + lim_w(2,2) = max(lim(1,2) - j, -this%hlf(2) + & + mod(n_stp(2)+this%knl(2)-j,this%stp(2))) + lim_w(1,1) = min(lim(2,1)-i, end_idx(1) - & + mod(n_stp(1)-1+i,this%stp(1))) + lim_w(1,2) = min(lim(2,2)-j, end_idx(2) - & + mod(n_stp(2)-1+j,this%stp(2))) +#endif if(any(lim_w(2,:).gt.lim_w(1,:))) cycle !! set gradient bounds +#if defined(GFORTRAN) lim_g(1,:) = max(1, [i,j] - offset) lim_g(2,:) = min(this%output_shape(:2), [i,j] - offset + this%knl - 1) +#else + lim_g(1,1) = max(1, i - offset(1)) + lim_g(1,2) = max(1, j - offset(2)) + lim_g(2,1) = min(this%output_shape(1), i - offset(1) + this%knl(1) -1) + lim_g(2,2) = min(this%output_shape(2), j - offset(2) + this%knl(2) -1) +#endif !! apply full convolution to compute input gradients this%di(i,j,m,s) = & diff --git a/src/lib/mod_conv3d_layer.f90 b/src/lib/mod_conv3d_layer.f90 index a5839aa..eb21137 100644 --- a/src/lib/mod_conv3d_layer.f90 +++ b/src/lib/mod_conv3d_layer.f90 @@ -4,6 +4,15 @@ !!!############################################################################# !!! module contains implementation of a 3D convolutional layer !!!############################################################################# +!!! Attribution statement: +!!! The following procedures are based on code from the neural-fortran library +!!! https://github.com/modern-fortran/neural-fortran/blob/main/src/nf/nf_layer.f90 +!!! procedures: +!!! - get_params* +!!! - set_params* +!!! - get_gradients* +!!! - set_gradients* +!!!############################################################################# module conv3d_layer use constants, only: real12 use base_layer, only: learnable_layer_type, conv_layer_type @@ -129,6 +138,7 @@ end subroutine layer_merge !!!############################################################################# !!! get learnable parameters +!!! procedure modified from neural-fortran library !!!############################################################################# pure function get_params_conv3d(this) result(params) implicit none @@ -146,6 +156,7 @@ end function get_params_conv3d !!!############################################################################# !!! set learnable parameters +!!! procedure modified from neural-fortran library !!!############################################################################# subroutine set_params_conv3d(this, params) implicit none @@ -165,6 +176,7 @@ end subroutine set_params_conv3d !!!############################################################################# !!! get sample-average gradients !!! sum over batch dimension and divide by batch size +!!! procedure modified from neural-fortran library !!!############################################################################# pure function get_gradients_conv3d(this, clip_method) result(gradients) use clipper, only: clip_type @@ -186,6 +198,7 @@ end function get_gradients_conv3d !!!############################################################################# !!! set gradients +!!! procedure modified from neural-fortran library !!!############################################################################# subroutine set_gradients_conv3d(this, gradients) implicit none @@ -994,19 +1007,43 @@ pure subroutine backward_5d(this, input, gradient) !! ... 1. offset of 1st o/p idx from centre of knl (lim) !! ... 2. lwst o/p idx overlap with <<- knl idx (rpt. pattern) !! ...) - lim_w(2,:) = max(lim(1,:)-[i,j,k], -this%hlf + & - mod(n_stp+this%knl-[i,j,k],this%stp)) !! min( ... !! ... 1. offset of last o/p idx from centre of knl (lim) !! ... 2. hghst o/p idx overlap with ->> knl idx (rpt. pattern) !! ...) +#if defined(GFORTRAN) + lim_w(2,:) = max(lim(1,:)-[i,j,k], -this%hlf + & + mod(n_stp+this%knl-[i,j,k],this%stp)) lim_w(1,:) = min(lim(2,:)-[i,j,k], end_idx - & mod(n_stp-1+[i,j,k],this%stp)) +#else + lim_w(2,1) = max(lim(1,1) - i, -this%hlf(1) + & + mod(n_stp(1)+this%knl(1)-i,this%stp(1))) + lim_w(2,2) = max(lim(1,2) - j, -this%hlf(2) + & + mod(n_stp(2)+this%knl(2)-j,this%stp(2))) + lim_w(2,3) = max(lim(1,3) - k, -this%hlf(3) + & + mod(n_stp(3)+this%knl(3)-k,this%stp(3))) + lim_w(1,1) = min(lim(2,1)-i, end_idx(1) - & + mod(n_stp(1)-1+i,this%stp(1))) + lim_w(1,2) = min(lim(2,2)-j, end_idx(2) - & + mod(n_stp(2)-1+j,this%stp(2))) + lim_w(1,3) = min(lim(2,3)-k, end_idx(3) - & + mod(n_stp(3)-1+j,this%stp(3))) +#endif if(any(lim_w(2,:).gt.lim_w(1,:))) cycle !! set gradient bounds +#if defined(GFORTRAN) lim_g(1,:) = max(1, [i,j,k] - offset) lim_g(2,:) = min(this%output_shape(:3), [i,j,k] - offset + this%knl-1) +#else + lim_g(1,1) = max(1, i - offset(1)) + lim_g(1,2) = max(1, j - offset(2)) + lim_g(1,3) = max(1, k - offset(3)) + lim_g(2,1) = min(this%output_shape(1), i - offset(1) + this%knl(1)-1) + lim_g(2,2) = min(this%output_shape(2), j - offset(2) + this%knl(2)-1) + lim_g(2,3) = min(this%output_shape(3), k - offset(3) + this%knl(3)-1) +#endif !! apply full convolution to compute input gradients this%di(i,j,k,m,s) = & diff --git a/src/lib/mod_full_layer.f90 b/src/lib/mod_full_layer.f90 index 159ed9b..3fec3e6 100644 --- a/src/lib/mod_full_layer.f90 +++ b/src/lib/mod_full_layer.f90 @@ -4,6 +4,16 @@ !!!############################################################################# !!! module contains implementation of a fully connected (dense) layer !!!############################################################################# +!!! Attribution statement: +!!! The following procedures are based on code from the neural-fortran library +!!! https://github.com/modern-fortran/neural-fortran/blob/main/src/nf/nf_layer.f90 +!!! procedures: +!!! - get_num_params* +!!! - get_params* +!!! - set_params* +!!! - get_gradients* +!!! - set_gradients* +!!!############################################################################# module full_layer use constants, only: real12 use base_layer, only: learnable_layer_type @@ -128,6 +138,7 @@ end subroutine layer_merge !!!############################################################################# !!! get number of parameters +!!! procedure modified from neural-fortran library !!!############################################################################# pure function get_num_params_full(this) result(num_params) implicit none @@ -142,6 +153,7 @@ end function get_num_params_full !!!############################################################################# !!! get number of parameters +!!! procedure modified from neural-fortran library !!!############################################################################# pure function get_params_full(this) result(params) implicit none @@ -156,6 +168,7 @@ end function get_params_full !!!############################################################################# !!! get number of parameters +!!! procedure modified from neural-fortran library !!!############################################################################# subroutine set_params_full(this, params) implicit none @@ -170,6 +183,7 @@ end subroutine set_params_full !!!############################################################################# !!! get number of parameters +!!! procedure modified from neural-fortran library !!!############################################################################# pure function get_gradients_full(this, clip_method) result(gradients) use clipper, only: clip_type @@ -189,6 +203,7 @@ end function get_gradients_full !!!############################################################################# !!! set gradients +!!! procedure modified from neural-fortran library !!!############################################################################# subroutine set_gradients_full(this, gradients) implicit none diff --git a/src/lib/mod_input2d_layer.f90 b/src/lib/mod_input2d_layer.f90 new file mode 100644 index 0000000..6a77dce --- /dev/null +++ b/src/lib/mod_input2d_layer.f90 @@ -0,0 +1,256 @@ +!!!############################################################################# +!!! Code written by Ned Thaddeus Taylor +!!! Code part of the ATHENA library - a feedforward neural network library +!!!############################################################################# +!!! module contains implementation of a 2D input layer +!!!############################################################################# +module input2d_layer + use constants, only: real12 + use base_layer, only: input_layer_type + implicit none + + + type, extends(input_layer_type) :: input2d_layer_type + real(real12), allocatable, dimension(:,:,:) :: output + contains + procedure, pass(this) :: get_output => get_output_input2d + procedure, pass(this) :: init => init_input2d + procedure, pass(this) :: set_batch_size => set_batch_size_input2d + procedure, pass(this) :: forward => forward_rank + procedure, pass(this) :: backward => backward_rank + procedure, pass(this) :: set => set_input2d + end type input2d_layer_type + + interface input2d_layer_type + module function layer_setup(input_shape, batch_size) result(layer) + integer, dimension(:), optional, intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + type(input2d_layer_type) :: layer + end function layer_setup + end interface input2d_layer_type + + + private + public :: input2d_layer_type + + +contains + +!!!############################################################################# +!!! get layer outputs +!!!############################################################################# + pure subroutine get_output_input2d(this, output) + implicit none + class(input2d_layer_type), intent(in) :: this + real(real12), allocatable, dimension(..), intent(out) :: output + + select rank(output) + rank(1) + output = reshape(this%output, [size(this%output)]) + rank(2) + output = & + reshape(this%output, [product(this%output_shape),this%batch_size]) + rank(3) + output = this%output + end select + + end subroutine get_output_input2d +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! forward propagation assumed rank handler +!!! placeholder to satisfy deferred +!!!############################################################################# + pure subroutine forward_rank(this, input) + implicit none + class(input2d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: input + + select rank(input) + rank(1) + this%output = reshape(input, shape=shape(this%output)) + rank(2) + this%output = reshape(input, shape=shape(this%output)) + rank(3) + this%output = reshape(input, shape=shape(this%output)) + rank(4) + this%output = reshape(input, shape=shape(this%output)) + rank(5) + this%output = reshape(input, shape=shape(this%output)) + rank(6) + this%output = reshape(input, shape=shape(this%output)) + end select + end subroutine forward_rank +!!!############################################################################# + + +!!!############################################################################# +!!! backward propagation assumed rank handler +!!! placeholder to satisfy deferred +!!!############################################################################# + pure subroutine backward_rank(this, input, gradient) + implicit none + class(input2d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: input + real(real12), dimension(..), intent(in) :: gradient + return + end subroutine backward_rank +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! set up layer +!!!############################################################################# +#if defined(GFORTRAN) + module function layer_setup(input_shape, batch_size) result(layer) + implicit none + integer, dimension(:), optional, intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + + type(input2d_layer_type) :: layer +#else + module procedure layer_setup + implicit none +#endif + + + layer%name = "input2d" + layer%input_rank = 2 + !!-------------------------------------------------------------------------- + !! initialise batch size + !!-------------------------------------------------------------------------- + if(present(batch_size)) layer%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! initialise layer shape + !!-------------------------------------------------------------------------- + if(present(input_shape)) call layer%init(input_shape=input_shape) + +#if defined(GFORTRAN) + end function layer_setup +#else + end procedure layer_setup +#endif +!!!############################################################################# + + +!!!############################################################################# +!!! initialise layer +!!!############################################################################# + subroutine init_input2d(this, input_shape, batch_size, verbose) + implicit none + class(input2d_layer_type), intent(inout) :: this + integer, dimension(:), intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, optional, intent(in) :: verbose + + integer :: verbose_ = 0 + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + if(present(batch_size)) this%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! initialise input shape + !!-------------------------------------------------------------------------- + if(.not.allocated(this%input_shape)) call this%set_shape(input_shape) + + this%output_shape = this%input_shape + this%num_outputs = product(this%input_shape) + + + !!-------------------------------------------------------------------------- + !! initialise batch size-dependent arrays + !!-------------------------------------------------------------------------- + if(this%batch_size.gt.0) call this%set_batch_size(this%batch_size) + + end subroutine init_input2d +!!!############################################################################# + + +!!!############################################################################# +!!! set batch size +!!!############################################################################# + subroutine set_batch_size_input2d(this, batch_size, verbose) + implicit none + class(input2d_layer_type), intent(inout) :: this + integer, intent(in) :: batch_size + integer, optional, intent(in) :: verbose + + integer :: verbose_ = 0 + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + this%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! allocate arrays + !!-------------------------------------------------------------------------- + if(allocated(this%input_shape))then + if(allocated(this%output)) deallocate(this%output) + allocate(this%output( & + this%input_shape(1), & + this%input_shape(2), & + this%batch_size)) + end if + + end subroutine set_batch_size_input2d +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! set input layer values +!!!############################################################################# + pure subroutine set_input2d(this, input) + implicit none + class(input2d_layer_type), intent(inout) :: this + real(real12), & + dimension(..), intent(in) :: input + !dimension(this%batch_size * this%num_outputs), intent(in) :: input + + select rank(input) + rank(1) + this%output = reshape(input, shape=shape(this%output)) + rank(2) + this%output = reshape(input, shape=shape(this%output)) + rank(3) + this%output = reshape(input, shape=shape(this%output)) + rank(4) + this%output = reshape(input, shape=shape(this%output)) + rank(5) + this%output = reshape(input, shape=shape(this%output)) + rank(6) + this%output = reshape(input, shape=shape(this%output)) + end select + + end subroutine set_input2d +!!!############################################################################# + + +end module input2d_layer +!!!############################################################################# \ No newline at end of file diff --git a/src/lib/mod_loss.f90 b/src/lib/mod_loss.f90 index 23a5a49..b17a6f3 100644 --- a/src/lib/mod_loss.f90 +++ b/src/lib/mod_loss.f90 @@ -185,6 +185,7 @@ pure function compute_loss_nll(predicted, expected) result(output) real(real12), dimension(size(predicted,1),size(predicted,2)) :: output real(real12) :: epsilon + epsilon = 1.E-10_real12 output = - log(expected - predicted + epsilon) end function compute_loss_nll diff --git a/src/lib/mod_maxpool1d_layer.f90 b/src/lib/mod_maxpool1d_layer.f90 new file mode 100644 index 0000000..78166e2 --- /dev/null +++ b/src/lib/mod_maxpool1d_layer.f90 @@ -0,0 +1,465 @@ +!!!############################################################################# +!!! Code written by Ned Thaddeus Taylor +!!! Code part of the ATHENA library - a feedforward neural network library +!!!############################################################################# +!!! module contains implementation of a 1D maxpooling layer +!!!############################################################################# +module maxpool1d_layer + use constants, only: real12 + use base_layer, only: pool_layer_type + implicit none + + + type, extends(pool_layer_type) :: maxpool1d_layer_type + real(real12), allocatable, dimension(:,:,:) :: output + real(real12), allocatable, dimension(:,:,:) :: di ! gradient of input (i.e. delta) + contains + procedure, pass(this) :: get_output => get_output_maxpool1d + procedure, pass(this) :: init => init_maxpool1d + procedure, pass(this) :: set_batch_size => set_batch_size_maxpool1d + procedure, pass(this) :: print => print_maxpool1d + procedure, pass(this) :: forward => forward_rank + procedure, pass(this) :: backward => backward_rank + procedure, private, pass(this) :: forward_3d + procedure, private, pass(this) :: backward_3d + end type maxpool1d_layer_type + + + interface maxpool1d_layer_type + module function layer_setup( & + input_shape, batch_size, & + pool_size, stride) result(layer) + integer, dimension(:), optional, intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, dimension(..), optional, intent(in) :: pool_size + integer, dimension(..), optional, intent(in) :: stride + type(maxpool1d_layer_type) :: layer + end function layer_setup + end interface maxpool1d_layer_type + + + private + public :: maxpool1d_layer_type + public :: read_maxpool1d_layer + + +contains + +!!!############################################################################# +!!! get layer outputs +!!!############################################################################# + pure subroutine get_output_maxpool1d(this, output) + implicit none + class(maxpool1d_layer_type), intent(in) :: this + real(real12), allocatable, dimension(..), intent(out) :: output + + select rank(output) + rank(1) + output = reshape(this%output, [size(this%output)]) + rank(2) + output = & + reshape(this%output, [product(this%output_shape),this%batch_size]) + rank(3) + output = this%output + end select + + end subroutine get_output_maxpool1d +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! forward propagation assumed rank handler +!!!############################################################################# + pure subroutine forward_rank(this, input) + implicit none + class(maxpool1d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: input + + select rank(input); rank(3) + call forward_3d(this, input) + end select + end subroutine forward_rank +!!!############################################################################# + + +!!!############################################################################# +!!! backward propagation assumed rank handler +!!!############################################################################# + pure subroutine backward_rank(this, input, gradient) + implicit none + class(maxpool1d_layer_type), intent(inout) :: this + real(real12), dimension(..), intent(in) :: input + real(real12), dimension(..), intent(in) :: gradient + + select rank(input); rank(3) + select rank(gradient); rank(3) + call backward_3d(this, input, gradient) + end select + end select + end subroutine backward_rank +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! set up layer +!!!############################################################################# +#if defined(GFORTRAN) + module function layer_setup( & + input_shape, batch_size, & + pool_size, stride) result(layer) + implicit none + integer, dimension(:), optional, intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, dimension(..), optional, intent(in) :: pool_size + integer, dimension(..), optional, intent(in) :: stride + + type(maxpool1d_layer_type) :: layer +#else + module procedure layer_setup + implicit none +#endif + + + layer%name = "maxpool1d" + layer%input_rank = 2 + allocate( & + layer%pool(layer%input_rank-1), & + layer%strd(layer%input_rank-1) ) + !!----------------------------------------------------------------------- + !! set up pool size + !!----------------------------------------------------------------------- + if(present(pool_size))then + select rank(pool_size) + rank(0) + layer%pool = pool_size + rank(1) + layer%pool = pool_size + end select + else + layer%pool = 2 + end if + + + !!----------------------------------------------------------------------- + !! set up stride + !!----------------------------------------------------------------------- + if(present(stride))then + select rank(stride) + rank(0) + layer%strd = stride + rank(1) + layer%strd = stride + end select + else + layer%strd = 2 + end if + + + !!-------------------------------------------------------------------------- + !! initialise batch size + !!-------------------------------------------------------------------------- + if(present(batch_size)) layer%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! initialise layer shape + !!-------------------------------------------------------------------------- + if(present(input_shape)) call layer%init(input_shape=input_shape) + +#if defined(GFORTRAN) + end function layer_setup +#else + end procedure layer_setup +#endif +!!!############################################################################# + + +!!!############################################################################# +!!! initialise layer +!!!############################################################################# + subroutine init_maxpool1d(this, input_shape, batch_size, verbose) + implicit none + class(maxpool1d_layer_type), intent(inout) :: this + integer, dimension(:), intent(in) :: input_shape + integer, optional, intent(in) :: batch_size + integer, optional, intent(in) :: verbose + + integer :: verbose_ = 0 + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + if(present(batch_size)) this%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! initialise input shape + !!-------------------------------------------------------------------------- + if(.not.allocated(this%input_shape)) call this%set_shape(input_shape) + + + !!----------------------------------------------------------------------- + !! set up number of channels, width, height + !!----------------------------------------------------------------------- + this%num_channels = this%input_shape(2) + allocate(this%output_shape(2)) + this%output_shape(2) = this%input_shape(2) + this%output_shape(1) = & + floor( (this%input_shape(1) - this%pool(1))/real(this%strd(1))) + 1 + + + !!-------------------------------------------------------------------------- + !! initialise batch size-dependent arrays + !!-------------------------------------------------------------------------- + if(this%batch_size.gt.0) call this%set_batch_size(this%batch_size) + + end subroutine init_maxpool1d +!!!############################################################################# + + +!!!############################################################################# +!!! set batch size +!!!############################################################################# + subroutine set_batch_size_maxpool1d(this, batch_size, verbose) + implicit none + class(maxpool1d_layer_type), intent(inout) :: this + integer, intent(in) :: batch_size + integer, optional, intent(in) :: verbose + + integer :: verbose_ = 0 + + + !!-------------------------------------------------------------------------- + !! initialise optional arguments + !!-------------------------------------------------------------------------- + if(present(verbose)) verbose_ = verbose + this%batch_size = batch_size + + + !!-------------------------------------------------------------------------- + !! allocate arrays + !!-------------------------------------------------------------------------- + if(allocated(this%input_shape))then + if(allocated(this%output)) deallocate(this%output) + allocate(this%output( & + this%output_shape(1), this%num_channels, & + this%batch_size), & + source=0._real12) + if(allocated(this%di)) deallocate(this%di) + allocate(this%di( & + this%input_shape(1), & + this%input_shape(2), & + this%batch_size), & + source=0._real12) + end if + + end subroutine set_batch_size_maxpool1d +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! print layer to file +!!!############################################################################# + subroutine print_maxpool1d(this, file) + implicit none + class(maxpool1d_layer_type), intent(in) :: this + character(*), intent(in) :: file + + integer :: unit + + !! open file with new unit + !!-------------------------------------------------------------------------- + open(newunit=unit, file=trim(file), access='append') + + !! write convolution initial parameters + !!-------------------------------------------------------------------------- + write(unit,'("MAXPOOL1D")') + write(unit,'(3X,"INPUT_SHAPE = ",3(1X,I0))') this%input_shape + write(unit,'(3X,"POOL_SIZE =",1X,I0)') this%pool(1) + write(unit,'(3X,"STRIDE =",1X,I0)') this%strd(1) + write(unit,'("END MAXPOOL1D")') + + !! close unit + !!-------------------------------------------------------------------------- + close(unit) + + end subroutine print_maxpool1d +!!!############################################################################# + + +!!!############################################################################# +!!! read layer from file +!!!############################################################################# + function read_maxpool1d_layer(unit) result(layer) + use infile_tools, only: assign_val, assign_vec + use misc, only: to_lower, icount + implicit none + integer, intent(in) :: unit + + class(maxpool1d_layer_type), allocatable :: layer + + integer :: stat + integer :: itmp1 + integer, dimension(1) :: pool_size, stride + integer, dimension(2) :: input_shape + character(256) :: buffer, tag + + + !! loop over tags in layer card + tag_loop: do + + !! check for end of file + read(unit,'(A)',iostat=stat) buffer + if(stat.ne.0)then + write(0,*) "ERROR: file encountered error (EoF?) before END MAXPOOL1D" + stop "Exiting..." + end if + if(trim(adjustl(buffer)).eq."") cycle tag_loop + + !! check for end of convolution card + if(trim(adjustl(buffer)).eq."END MAXPOOL1D")then + backspace(unit) + exit tag_loop + end if + + tag=trim(adjustl(buffer)) + if(scan(buffer,"=").ne.0) tag=trim(tag(:scan(tag,"=")-1)) + + !! read parameters from save file + select case(trim(tag)) + case("INPUT_SHAPE") + call assign_vec(buffer, input_shape, itmp1) + case("POOL_SIZE") + call assign_vec(buffer, pool_size, itmp1) + case("STRIDE") + call assign_vec(buffer, stride, itmp1) + case default + !! don't look for "e" due to scientific notation of numbers + !! ... i.e. exponent (E+00) + if(scan(to_lower(trim(adjustl(buffer))),& + 'abcdfghijklmnopqrstuvwxyz').eq.0)then + cycle tag_loop + elseif(tag(:3).eq.'END')then + cycle tag_loop + end if + stop "Unrecognised line in input file: "//trim(adjustl(buffer)) + end select + end do tag_loop + + !! set transfer activation function + + layer = maxpool1d_layer_type( input_shape=input_shape, & + pool_size = pool_size, stride = stride & + ) + + !! check for end of layer card + read(unit,'(A)') buffer + if(trim(adjustl(buffer)).ne."END MAXPOOL1D")then + write(*,*) trim(adjustl(buffer)) + stop "ERROR: END MAXPOOL1D not where expected" + end if + + end function read_maxpool1d_layer +!!!############################################################################# + + +!!!##########################################################################!!! +!!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * !!! +!!!##########################################################################!!! + + +!!!############################################################################# +!!! forward propagation +!!!############################################################################# + pure subroutine forward_3d(this, input) + implicit none + class(maxpool1d_layer_type), intent(inout) :: this + real(real12), dimension( & + this%input_shape(1), & + this%num_channels, & + this%batch_size), & + intent(in) :: input + + integer :: i, m, s + integer :: stride_idx + + + this%output = 0._real12 + !! perform the pooling operation + do concurrent(& + s = 1:this%batch_size, & + m = 1:this%num_channels, & + i = 1:this%output_shape(1)) + stride_idx = (i-1) * this%strd(1) + 1 + this%output(i, m, s) = maxval(& + input( & + stride_idx:stride_idx+this%pool(1)-1, m, s)) + end do + + end subroutine forward_3d +!!!############################################################################# + + +!!!############################################################################# +!!! backward propagation +!!!############################################################################# + pure subroutine backward_3d(this, input, gradient) + implicit none + class(maxpool1d_layer_type), intent(inout) :: this + real(real12), dimension( & + this%input_shape(1), & + this%num_channels, & + this%batch_size), & + intent(in) :: input + real(real12), & + dimension(& + this%output_shape(1), & + this%num_channels, & + this%batch_size), & + intent(in) :: gradient + + integer :: i, m, s + integer :: stride_idx, max_idx + + + this%di = 0._real12 + !! compute gradients for input feature map + do concurrent( & + s = 1:this%batch_size, & + m = 1:this%num_channels, & + i = 1:this%output_shape(1)) + stride_idx = (i-1) * this%strd(1) + !! find the index of the maximum value in the corresponding pooling window + max_idx = maxloc(input( & + stride_idx+1:stride_idx+this%pool(1), m, s), dim=1) + + !! compute gradients for input feature map + this%di( & + stride_idx+max_idx, m, s) = & + this%di( & + stride_idx+max_idx, m, s) + gradient(i, m, s) + + end do + + end subroutine backward_3d +!!!############################################################################# + +end module maxpool1d_layer +!!!############################################################################# diff --git a/src/lib/mod_network_sub.f90 b/src/lib/mod_network_sub.f90 index ca4f8d9..cfd917d 100644 --- a/src/lib/mod_network_sub.f90 +++ b/src/lib/mod_network_sub.f90 @@ -25,6 +25,7 @@ !! input layer types use input1d_layer, only: input1d_layer_type + use input2d_layer, only: input2d_layer_type use input3d_layer, only: input3d_layer_type use input4d_layer, only: input4d_layer_type @@ -34,6 +35,7 @@ use batchnorm3d_layer, only: batchnorm3d_layer_type, read_batchnorm3d_layer !! convolution layer types + use conv1d_layer, only: conv1d_layer_type, read_conv1d_layer use conv2d_layer, only: conv2d_layer_type, read_conv2d_layer use conv3d_layer, only: conv3d_layer_type, read_conv3d_layer @@ -43,8 +45,10 @@ use dropblock3d_layer, only: dropblock3d_layer_type, read_dropblock3d_layer !! pooling layer types + use avgpool1d_layer, only: avgpool1d_layer_type, read_avgpool1d_layer use avgpool2d_layer, only: avgpool2d_layer_type, read_avgpool2d_layer use avgpool3d_layer, only: avgpool3d_layer_type, read_avgpool3d_layer + use maxpool1d_layer, only: maxpool1d_layer_type, read_maxpool1d_layer use maxpool2d_layer, only: maxpool2d_layer_type, read_maxpool2d_layer use maxpool3d_layer, only: maxpool3d_layer_type, read_maxpool3d_layer @@ -171,6 +175,8 @@ module subroutine read(this, file) call this%add(read_batchnorm2d_layer(unit)) case("BATCHNORM3D") call this%add(read_batchnorm3d_layer(unit)) + case("CONV1D") + call this%add(read_conv1d_layer(unit)) case("CONV2D") call this%add(read_conv2d_layer(unit)) case("CONV3D") @@ -181,10 +187,14 @@ module subroutine read(this, file) call this%add(read_dropblock2d_layer(unit)) case("DROPBLOCK3D") call this%add(read_dropblock3d_layer(unit)) + case("AVGPOOL1D") + call this%add(read_avgpool1d_layer(unit)) case("AVGPOOL2D") call this%add(read_avgpool2d_layer(unit)) case("AVGPOOL3D") call this%add(read_avgpool3d_layer(unit)) + case("MAXPOOL1D") + call this%add(read_maxpool1d_layer(unit)) case("MAXPOOL2D") call this%add(read_maxpool2d_layer(unit)) case("MAXPOOL3D") @@ -490,6 +500,18 @@ module subroutine compile(this, optimiser, loss_method, metrics, batch_size, & case(1) t_input_layer = input1d_layer_type(input_shape = next%input_shape) allocate(this%model(1)%layer, source = t_input_layer) + case(2) + select type(next) + type is(conv1d_layer_type) + t_input_layer = input2d_layer_type(& + input_shape = next%input_shape + & + [2*next%pad,0]) + allocate(this%model(1)%layer, source = t_input_layer) + class default + t_input_layer = input2d_layer_type(& + input_shape = next%input_shape) + allocate(this%model(1)%layer, source = t_input_layer) + end select case(3) select type(next) type is(conv2d_layer_type) @@ -524,9 +546,7 @@ module subroutine compile(this, optimiser, loss_method, metrics, batch_size, & !!! ignore calcuation of input gradients for 1st non-input layer !!!----------------------------------------------------------------------------- select type(second => this%model(2)%layer) - type is(conv2d_layer_type) - second%calc_input_gradients = .false. - type is(conv3d_layer_type) + class is(conv_layer_type) second%calc_input_gradients = .false. end select @@ -941,6 +961,8 @@ pure module subroutine backward_1d(this, output) type is(batchnorm3d_layer_type) call this%model(i)%backward(this%model(i-1),next%di) + type is(conv1d_layer_type) + call this%model(i)%backward(this%model(i-1),next%di) type is(conv2d_layer_type) call this%model(i)%backward(this%model(i-1),next%di) type is(conv3d_layer_type) @@ -953,10 +975,14 @@ pure module subroutine backward_1d(this, output) type is(dropblock3d_layer_type) call this%model(i)%backward(this%model(i-1),next%di) + type is(avgpool1d_layer_type) + call this%model(i)%backward(this%model(i-1),next%di) type is(avgpool2d_layer_type) call this%model(i)%backward(this%model(i-1),next%di) type is(avgpool3d_layer_type) call this%model(i)%backward(this%model(i-1),next%di) + type is(maxpool1d_layer_type) + call this%model(i)%backward(this%model(i-1),next%di) type is(maxpool2d_layer_type) call this%model(i)%backward(this%model(i-1),next%di) type is(maxpool3d_layer_type) diff --git a/src/lib/mod_optimiser.f90 b/src/lib/mod_optimiser.f90 index e01f87f..5d71ffc 100644 --- a/src/lib/mod_optimiser.f90 +++ b/src/lib/mod_optimiser.f90 @@ -15,6 +15,14 @@ !!! - minimise - minimise the loss function by applying gradients to ... !!! ... the parameters !!!############################################################################# +!!! Attribution statement: +!!! The following module is based on code from the neural-fortran library +!!! https://github.com/modern-fortran/neural-fortran/blob/main/src/nf/nf_optimizers.f90 +!!! The implementation of optimiser_base_type is based on the ... +!!! ... optimizer_base_type from the neural-fortran library +!!! The same applies to the implementation of the sgd_optimiser_type, ... +!!! ... rmsprop_optimiser_type, adagrad_optimiser_type, and adam_optimiser_type +!!!############################################################################# module optimiser use constants, only: real12 use clipper, only: clip_type @@ -218,6 +226,8 @@ module function optimiser_setup_base( & if(present(lr_decay)) then if(allocated(optimiser%lr_decay)) deallocate(optimiser%lr_decay) allocate(optimiser%lr_decay, source = lr_decay) + else + allocate(optimiser%lr_decay, source = base_lr_decay_type()) end if !! initialise gradients @@ -342,6 +352,8 @@ module function optimiser_setup_sgd( & if(present(lr_decay)) then if(allocated(optimiser%lr_decay)) deallocate(optimiser%lr_decay) allocate(optimiser%lr_decay, source = lr_decay) + else + allocate(optimiser%lr_decay, source = base_lr_decay_type()) end if !! initialise nesterov boolean @@ -459,6 +471,8 @@ module function optimiser_setup_rmsprop( & if(present(lr_decay)) then if(allocated(optimiser%lr_decay)) deallocate(optimiser%lr_decay) allocate(optimiser%lr_decay, source = lr_decay) + else + allocate(optimiser%lr_decay, source = base_lr_decay_type()) end if !! initialise adam parameters @@ -569,6 +583,8 @@ module function optimiser_setup_adagrad( & if(present(lr_decay)) then if(allocated(optimiser%lr_decay)) deallocate(optimiser%lr_decay) allocate(optimiser%lr_decay, source = lr_decay) + else + allocate(optimiser%lr_decay, source = base_lr_decay_type()) end if !! initialise adam parameters @@ -677,6 +693,8 @@ module function optimiser_setup_adam( & if(present(lr_decay)) then if(allocated(optimiser%lr_decay)) deallocate(optimiser%lr_decay) allocate(optimiser%lr_decay, source = lr_decay) + else + allocate(optimiser%lr_decay, source = base_lr_decay_type()) end if !! initialise adam parameters diff --git a/src/lib/mod_random.f90 b/src/lib/mod_random.f90 index c63a0a6..76d0438 100644 --- a/src/lib/mod_random.f90 +++ b/src/lib/mod_random.f90 @@ -23,7 +23,7 @@ module random subroutine random_setup(seed, num_seed, restart, already_initialised) implicit none integer, dimension(..), optional, intent(in) :: seed !dimension(..1) - integer, optional, intent(in) :: num_seed + integer, optional, intent(out) :: num_seed logical, optional, intent(in) :: restart logical, optional, intent(out) :: already_initialised @@ -41,30 +41,6 @@ subroutine random_setup(seed, num_seed, restart, already_initialised) end if if(present(already_initialised)) already_initialised = .false. - !! define number of seeds - if(present(num_seed))then - if(present(seed))then - select rank(seed) - rank(0) - num_seed_ = num_seed - rank(1) - if(size(seed,dim=1).ne.1.and.size(seed,dim=1).ne.num_seed)then - write(0,*) "ERROR: seed and num_seed provided to random_setup" - write(0,*) " Cannot decide which to listen to" - stop "Exiting..." - end if - end select - else - num_seed_ = num_seed - end if - else - if(present(seed))then - num_seed_ = size(seed,dim=1) - else - num_seed_ = 1 - end if - end if - !! check if already initialised if(l_random_initialised.and..not.restart_)then if(present(already_initialised)) already_initialised = .true. @@ -77,8 +53,16 @@ subroutine random_setup(seed, num_seed, restart, already_initialised) rank(0) seed_arr = seed rank(1) - if(size(seed,dim=1).gt.1)then - seed_arr = seed + if(size(seed,dim=1).ne.1)then + if(size(seed,dim=1).eq.num_seed_)then + seed_arr = seed + else + write(0,*) "ERROR: seed size not consistent with & + &seed size returned by implementation" + write(0,*) "Cannot resolve" + write(0,*) "Exiting..." + stop 1 + end if else seed_arr = seed(1) end if @@ -90,6 +74,8 @@ subroutine random_setup(seed, num_seed, restart, already_initialised) call random_seed(put=seed_arr) l_random_initialised = .true. end if + + if(present(num_seed)) num_seed = num_seed_ end subroutine random_setup !!!############################################################################# diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5bcb727..d2d264a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,14 +1,18 @@ message(STATUS "Building tests") foreach(execid input_layer + conv1d_layer + conv1d_network conv2d_layer conv2d_network conv3d_layer conv3d_network full_layer full_network + maxpool1d_layer maxpool2d_layer maxpool3d_layer + avgpool1d_layer avgpool2d_layer avgpool3d_layer batchnorm1d_layer diff --git a/test/test_activations.f90 b/test/test_activations.f90 index 4fa1ff1..3a49542 100644 --- a/test/test_activations.f90 +++ b/test/test_activations.f90 @@ -14,7 +14,7 @@ program test_activations implicit none class(base_layer_type), allocatable :: full_layer, conv2d_layer, conv3d_layer - class(activation_type), allocatable :: activation + class(activation_type), allocatable :: activation_var logical :: success = .true. integer :: i @@ -37,6 +37,9 @@ program test_activations real, dimension(9) :: activate, differentiate real, dimension(1) :: value_1d, rtmp1_1d real, dimension(1,1,1) :: value_3d, rtmp1_3d + real, dimension(:,:), allocatable :: output_2d + real, dimension(:,:,:,:), allocatable :: output_4d + real, dimension(:,:,:,:,:), allocatable :: output_5d !!!----------------------------------------------------------------------------- @@ -83,12 +86,12 @@ program test_activations !!!----------------------------------------------------------------------------- !!! check gaussian setup !!!----------------------------------------------------------------------------- - activation = gaussian_setup(threshold = 2.E0, sigma = 2.E0) - if(.not. activation%name .eq. 'gaussian')then + activation_var = gaussian_setup(threshold = 2.E0, sigma = 2.E0) + if(.not. activation_var%name .eq. 'gaussian')then success = .false. write(0,*) 'activation has wrong name for gaussian' else - if (activation%threshold .ne. 2.E0) then + if (activation_var%threshold .ne. 2.E0) then success = .false. write(0,*) 'activation has wrong threshold for gaussian' end if @@ -98,8 +101,8 @@ program test_activations !!!----------------------------------------------------------------------------- !!! check piecewise setup !!!----------------------------------------------------------------------------- - activation = piecewise_setup(intercept = 2.E0) - if(.not. activation%name .eq. 'piecewise')then + activation_var = piecewise_setup(intercept = 2.E0) + if(.not. activation_var%name .eq. 'piecewise')then success = .false. write(0,*) 'activation has wrong name for piecewise' end if @@ -108,12 +111,12 @@ program test_activations !!!----------------------------------------------------------------------------- !!! check sigmoid setup !!!----------------------------------------------------------------------------- - activation = sigmoid_setup(threshold = 2.E0) - if(.not. activation%name .eq. 'sigmoid')then + activation_var = sigmoid_setup(threshold = 2.E0) + if(.not. activation_var%name .eq. 'sigmoid')then success = .false. write(0,*) 'activation has wrong name for sigmoid' else - if (activation%threshold .ne. 2.E0) then + if (activation_var%threshold .ne. 2.E0) then success = .false. write(0,*) 'activation has wrong threshold for sigmoid' end if @@ -123,12 +126,12 @@ program test_activations !!!----------------------------------------------------------------------------- !!! check softmax setup !!!----------------------------------------------------------------------------- - activation = softmax_setup(threshold = 2.E0) - if(.not. activation%name .eq. 'softmax')then + activation_var = softmax_setup(threshold = 2.E0) + if(.not. activation_var%name .eq. 'softmax')then success = .false. write(0,*) 'activation has wrong name for softmax' else - if (activation%threshold .ne. 2.E0) then + if (activation_var%threshold .ne. 2.E0) then success = .false. write(0,*) 'activation has wrong threshold for softmax' end if @@ -138,12 +141,12 @@ program test_activations !!!----------------------------------------------------------------------------- !!! check tanh setup !!!----------------------------------------------------------------------------- - activation = tanh_setup(threshold = 2.E0) - if(.not. activation%name .eq. 'tanh')then + activation_var = tanh_setup(threshold = 2.E0) + if(.not. activation_var%name .eq. 'tanh')then success = .false. write(0,*) 'activation has wrong name for tanh' else - if (activation%threshold .ne. 2.E0) then + if (activation_var%threshold .ne. 2.E0) then success = .false. write(0,*) 'activation has wrong threshold for tanh' end if @@ -154,27 +157,27 @@ program test_activations !!! check for different scales, and ranks !!!----------------------------------------------------------------------------- do i = 1, size(activation_names) - deallocate(activation) - allocate(activation, & + deallocate(activation_var) + allocate(activation_var, & source=activation_setup(activation_names(i), scale = scale)) - if(.not. activation%name .eq. trim(activation_names(i)))then + if(.not. activation_var%name .eq. trim(activation_names(i)))then success = .false. write(0,*) 'activation has wrong name for ', & trim(activation_names(i)) else - if (abs(activation%scale - scale).gt.1.E-6) then + if (abs(activation_var%scale - scale).gt.1.E-6) then success = .false. write(0,*) 'activation has wrong scale for ', & trim(activation_names(i)) end if - rtmp1_1d = activation%activate_1d(value_1d) + rtmp1_1d = activation_var%activate_1d(value_1d) if (any(abs(rtmp1_1d - activate(i)).gt.1.E-6)) then success = .false. write(0,*) 'activation has wrong activation for ', & trim(activation_names(i)) write(*,*) rtmp1_1d, activate(i) end if - rtmp1_1d = activation%differentiate_1d(value_1d) + rtmp1_1d = activation_var%differentiate_1d(value_1d) if (any(abs(rtmp1_1d - differentiate(i)).gt.1.E-6)) then success = .false. write(0,*) 'activation has wrong differentiation for ', & @@ -182,13 +185,13 @@ program test_activations write(*,*) rtmp1_1d, differentiate(i) end if - rtmp1_3d = activation%activate_3d(value_3d) + rtmp1_3d = activation_var%activate_3d(value_3d) if (any(abs(rtmp1_3d - activate(i)).gt.1.E-6)) then success = .false. write(0,*) 'activation has wrong activation for ', & trim(activation_names(i)) end if - rtmp1_3d = activation%differentiate_3d(value_3d) + rtmp1_3d = activation_var%differentiate_3d(value_3d) if (any(abs(rtmp1_3d - differentiate(i)).gt.1.E-6)) then success = .false. write(0,*) 'activation has wrong differentiation for ', & @@ -220,8 +223,9 @@ program test_activations full_layer%output, input_data, activation_names(i), success) full_layer%output = 1.E0 + output_2d = full_layer%transfer%differentiate(full_layer%output) call compare_derivative( & - full_layer%transfer%differentiate(full_layer%output), & + output_2d, & full_layer%output, & activation_names(i), success) end if @@ -259,8 +263,9 @@ program test_activations input_data_conv2d, activation_names(i), success) conv2d_layer%output = 1.E0 + output_4d = conv2d_layer%transfer%differentiate(conv2d_layer%output) call compare_derivative( & - conv2d_layer%transfer%differentiate(conv2d_layer%output), & + output_4d, & conv2d_layer%output, & activation_names(i), success) end if @@ -298,8 +303,9 @@ program test_activations input_data_conv3d, activation_names(i), success) conv3d_layer%output = 1.E0 + output_5d = conv3d_layer%transfer%differentiate(conv3d_layer%output) call compare_derivative( & - conv3d_layer%transfer%differentiate(conv3d_layer%output), & + output_5d, & conv3d_layer%output, & activation_names(i), success) end if diff --git a/test/test_avgpool1d_layer.f90 b/test/test_avgpool1d_layer.f90 new file mode 100644 index 0000000..f3fbfd2 --- /dev/null +++ b/test/test_avgpool1d_layer.f90 @@ -0,0 +1,222 @@ +program test_avgpool1d_layer + use athena, only: & + avgpool1d_layer_type, & + base_layer_type, & + learnable_layer_type + implicit none + + class(base_layer_type), allocatable :: pool_layer + integer, parameter :: num_channels = 3, pool = 3, stride = 3, width = 9 + real, allocatable, dimension(:) :: output_1d + real, allocatable, dimension(:,:) :: output_2d + real, allocatable, dimension(:,:,:) :: input_data, output, gradient, & + di_compare + real, parameter :: tol = 1.E-7 + logical :: success = .true. + + integer :: i, j, output_width, max_loc + integer :: ip1, ip2 + real, parameter :: max_value = 3.0 + + + !! set up avgpool1d layer + pool_layer = avgpool1d_layer_type( & + pool_size = pool, & + stride = stride & + ) + + !! check layer name + if(.not. pool_layer%name .eq. 'avgpool1d')then + success = .false. + write(0,*) 'avgpool1d layer has wrong name' + end if + +!!!----------------------------------------------------------------------------- + + !! check layer type + select type(pool_layer) + type is(avgpool1d_layer_type) + !! check pool size + if(any(pool_layer%pool .ne. pool))then + success = .false. + write(0,*) 'avgpool1d layer has wrong pool size' + end if + + !! check stride size + if(any(pool_layer%strd .ne. stride))then + success = .false. + write(0,*) 'avgpool1d layer has wrong stride size' + end if + + !! check input shape allocated + if(allocated(pool_layer%input_shape))then + success = .false. + write(0,*) 'avgpool1d layer shape should not be allocated yet' + end if + class default + success = .false. + write(0,*) 'avgpool1d layer has wrong type' + end select + +!!!----------------------------------------------------------------------------- + + !! initialise width and output width + output_width = floor( (width - pool)/real(stride)) + 1 + max_loc = width / 2 + mod(width, 2) + + !! initialise sample input + allocate(input_data(width, num_channels, 1), source = 0.0) + input_data(max_loc, :, 1) = max_value + pool_layer = avgpool1d_layer_type( & + pool_size = pool, & + stride = stride & + ) + +!!!----------------------------------------------------------------------------- + + !! check layer input and output shape based on input layer + call pool_layer%init(shape(input_data(:,:,1)), batch_size=1) + select type(pool_layer) + type is(avgpool1d_layer_type) + if(any(pool_layer%input_shape .ne. [width,num_channels]))then + success = .false. + write(0,*) 'avgpool1d layer has wrong input_shape' + end if + if(any(pool_layer%output_shape .ne. [output_width,num_channels]))then + success = .false. + write(0,*) 'avgpool1d layer has wrong output_shape', pool_layer%output_shape + write(0,*) 'expected', [output_width,num_channels] + end if + end select + + !! run forward pass + call pool_layer%forward(input_data) + call pool_layer%get_output(output) + + !! check outputs have expected value + do i = 1, output_width + if( max_loc .ge. (i-1)*stride + 1 .and. & + max_loc .le. (i-1)*stride + pool )then + if(output(i, 1, 1) .ne. max_value/(pool))then + success = .false. + write(*,*) 'avgpool1d layer forward pass failed' + end if + else if(output(i, 1, 1) .ne. 0.0) then + success = .false. + write(*,*) 'avgpool1d layer forward pass failed' + end if + end do + + !! check 1d and 2d output are the same + call pool_layer%get_output(output_1d) + call pool_layer%get_output(output_2d) + if(any(abs(output_1d - & + reshape(output, [output_width*num_channels])) & + .gt. 1.E-6))then + success = .false. + write(*,*) 'avgpool1d layer output pass failed' + end if + if(any(abs(& + reshape(output_2d, [output_width*num_channels]) - & + reshape(output, [output_width*num_channels])) & + .gt. 1.E-6))then + success = .false. + write(*,*) 'avgpool1d layer output pass failed' + end if + +!!!----------------------------------------------------------------------------- + + !! run backward pass + allocate(gradient, source = output) + call pool_layer%backward(input_data, gradient) + allocate(di_compare(width,num_channels,1), source = 0.0) + do i = 1, output_width + ip1 = (i-1) * stride + 1 + ip2 = (i-1) * stride + pool + di_compare(ip1:ip2,1,1) = & + di_compare(ip1:ip2,1,1) + gradient(i,1,1)/pool + end do + + !! check gradient has expected value + select type(current => pool_layer) + type is(avgpool1d_layer_type) + if(any(abs(current%di(:,1,1) - di_compare(:,1,1)) .gt. tol))then + success = .false. + write(*,*) 'avgpool1d layer backward pass failed' + end if + end select + + !! check backward pass recovers input (with division by pool) + !! https://stats.stackexchange.com/questions/565032/cnn-upsampling-backprop-gradients-across-average-pooling-layer + call pool_layer%forward(di_compare) + call pool_layer%get_output(output) + + !! check outputs have expected value + if (any(abs(output(:,1,1) - (gradient(:,1,1))/real(pool)) .gt. tol)) then + success = .false. + write(*,*) 'avgpool1d layer forward pass failed' + do i = 1, width + write(*,'(18(1X,F7.3))') input_data(i,1,1) + end do + write(*,*) "----------------------------------------" + do i = 1, width + write(*,'(18(1X,F7.3))') di_compare(i,1,1) + end do + write(*,*) "----------------------------------------" + write(*,*) "----------------------------------------" + do i = 1, output_width + write(*,'(3(1X,F7.5))') gradient(i,1,1)/pool + end do + write(*,*) "----------------------------------------" + do i = 1, output_width + write(*,'(3(1X,F7.5))') output(i,1,1) + end do + + end if + + !! check expected initialisation of pool and stride + pool_layer = avgpool1d_layer_type( & + pool_size = [2], & + stride = [2] & + ) + select type(pool_layer) + type is (avgpool1d_layer_type) + if(any(pool_layer%pool .ne. [2]))then + success = .false. + write(0,*) 'avgpool1d layer has wrong pool size' + end if + if(any(pool_layer%strd .ne. [2]))then + success = .false. + write(0,*) 'avgpool1d layer has wrong stride size' + end if + end select + + !! check expected initialisation of pool and stride + pool_layer = avgpool1d_layer_type( & + pool_size = [4], & + stride = [4] & + ) + select type(pool_layer) + type is (avgpool1d_layer_type) + if(any(pool_layer%pool .ne. 4))then + success = .false. + write(0,*) 'avgpool1d layer has wrong pool size' + end if + if(any(pool_layer%strd .ne. 4))then + success = .false. + write(0,*) 'avgpool1d layer has wrong stride size' + end if + end select + +!!!----------------------------------------------------------------------------- +!!! check for any failed tests +!!!----------------------------------------------------------------------------- + write(*,*) "----------------------------------------" + if(success)then + write(*,*) 'test_avgpool1d_layer passed all tests' + else + write(0,*) 'test_avgpool1d_layer failed one or more tests' + stop 1 + end if + +end program test_avgpool1d_layer \ No newline at end of file diff --git a/test/test_batchnorm1d_layer.f90 b/test/test_batchnorm1d_layer.f90 index 5a7b130..f22c3a6 100644 --- a/test/test_batchnorm1d_layer.f90 +++ b/test/test_batchnorm1d_layer.f90 @@ -14,7 +14,7 @@ program test_batchnorm1d_layer logical :: success = .true. integer :: i, j, output_width, num_params - integer :: seed_size = 1 + integer :: seed_size real :: mean, std integer, allocatable, dimension(:) :: seed real, parameter :: max_value = 3.0 diff --git a/test/test_batchnorm2d_layer.f90 b/test/test_batchnorm2d_layer.f90 index 5382dd4..640ca5a 100644 --- a/test/test_batchnorm2d_layer.f90 +++ b/test/test_batchnorm2d_layer.f90 @@ -15,7 +15,7 @@ program test_batchnorm2d_layer logical :: success = .true. integer :: i, j, output_width, num_params - integer :: seed_size = 1 + integer :: seed_size real :: mean, std integer, allocatable, dimension(:) :: seed real, parameter :: max_value = 3.0 diff --git a/test/test_batchnorm3d_layer.f90 b/test/test_batchnorm3d_layer.f90 index c1e8174..a73f798 100644 --- a/test/test_batchnorm3d_layer.f90 +++ b/test/test_batchnorm3d_layer.f90 @@ -15,7 +15,7 @@ program test_batchnorm3d_layer logical :: success = .true. integer :: i, j, output_width, num_params - integer :: seed_size = 1 + integer :: seed_size real :: mean, std integer, allocatable, dimension(:) :: seed real, parameter :: max_value = 3.0 @@ -287,4 +287,4 @@ subroutine compare_batchnorm3d_layers(layer1, layer2, layer3, success) end subroutine compare_batchnorm3d_layers -end program test_batchnorm3d_layer \ No newline at end of file +end program test_batchnorm3d_layer diff --git a/test/test_conv1d_layer.f90 b/test/test_conv1d_layer.f90 new file mode 100644 index 0000000..bdf5b6c --- /dev/null +++ b/test/test_conv1d_layer.f90 @@ -0,0 +1,312 @@ +program test_conv1d_layer + use athena, only: & + conv1d_layer_type, & + input2d_layer_type, & + base_layer_type, & + learnable_layer_type + implicit none + + class(base_layer_type), allocatable :: conv_layer, conv_layer1, conv_layer2 + class(base_layer_type), allocatable :: input_layer + integer, parameter :: num_filters = 32, kernel_size = 3 + real, allocatable, dimension(:,:,:) :: input_data, output + real, parameter :: tol = 1.E-7 + logical :: success = .true. + + real, allocatable, dimension(:) :: params + real, allocatable, dimension(:,:) :: outputs + + +!!!----------------------------------------------------------------------------- +!!! set up layer +!!!----------------------------------------------------------------------------- + conv_layer = conv1d_layer_type( & + num_filters = num_filters, & + kernel_size = kernel_size & + ) + + !! check layer name + if(.not. conv_layer%name .eq. 'conv1d')then + success = .false. + write(0,*) 'conv1d layer has wrong name' + end if + + !! check layer type + select type(conv_layer) + type is(conv1d_layer_type) + !! check default layer transfer/activation function + if(conv_layer%transfer%name .ne. 'none')then + success = .false. + write(0,*) 'conv1d layer has wrong transfer: '//conv_layer%transfer%name + end if + + !! check number of filters + if(conv_layer%num_filters .ne. num_filters)then + success = .false. + write(0,*) 'conv1d layer has wrong num_filters' + end if + + !! check kernel size + if(any(conv_layer%knl .ne. kernel_size))then + success = .false. + write(0,*) 'conv1d layer has wrong kernel_size' + end if + + !! check input shape allocated + if(allocated(conv_layer%input_shape))then + success = .false. + write(0,*) 'conv1d layer shape should not be allocated yet' + end if + class default + success = .false. + write(0,*) 'conv1d layer has wrong type' + end select + + +!!!----------------------------------------------------------------------------- +!!! check layer input and output shape based on input layer +!!! conv1d layer: 32 x 32 pixel image, 3 channels +!!!----------------------------------------------------------------------------- + input_layer = input2d_layer_type([32,3]) + call conv_layer%init(input_layer%input_shape) + select type(conv_layer) + type is(conv1d_layer_type) + if(any(conv_layer%input_shape .ne. [32,3]))then + success = .false. + write(0,*) 'conv1d layer has wrong input_shape' + end if + if(any(conv_layer%output_shape .ne. [30,num_filters]))then + success = .false. + write(0,*) 'conv1d layer has wrong output_shape' + end if + end select + + +!!!----------------------------------------------------------------------------- +!!! test forward pass and check expected output +!!! use existing layer +!!!----------------------------------------------------------------------------- + !! initialise sample input + !! conv1d layer: 3x3 pixel image, 1 channel + allocate(input_data(3, 1, 1), source = 0.0) + input_layer = input2d_layer_type([3,1], batch_size=1) + conv_layer = conv1d_layer_type( & + num_filters = 1, & + kernel_size = 3, & + activation_function = 'sigmoid' & + ) + call conv_layer%init(input_layer%input_shape, batch_size=1) + + !! set input data in input_layer + select type (current => input_layer) + type is(input2d_layer_type) + call current%set(input_data) + end select + call input_layer%get_output(output) + + !! run forward pass + call conv_layer%forward(input_data) + call conv_layer%get_output(output) + + !! check outputs have expected value + if (any(abs(output - 0.5).gt.tol)) then + success = .false. + write(*,*) 'conv1d layer with zero input and sigmoid activation must & + &return outputs all equal to 0.5' + write(*,*) output + end if + + +!!!----------------------------------------------------------------------------- +!!! check handling of layer parameters, gradients, and outputs +!!!----------------------------------------------------------------------------- + select type(conv_layer) + class is(learnable_layer_type) + !! check layer parameter handling + params = conv_layer%get_params() + if(size(params) .eq. 0)then + success = .false. + write(0,*) 'conv1d layer has wrong number of parameters' + end if + params = 1.E0 + call conv_layer%set_params(params) + params = conv_layer%get_params() + if(any(abs(params - 1.E0).gt.tol))then + success = .false. + write(0,*) 'conv1d layer has wrong parameters' + end if + + !! check layer gradient handling + params = conv_layer%get_gradients() + if(size(params) .eq. 0)then + success = .false. + write(0,*) 'conv1d layer has wrong number of gradients' + end if + params = 1.E0 + call conv_layer%set_gradients(params) + params = conv_layer%get_gradients() + if(any(abs(params - 1.E0).gt.tol))then + success = .false. + write(0,*) 'conv1d layer has wrong gradients' + end if + call conv_layer%set_gradients(10.E0) + params = conv_layer%get_gradients() + if(any(abs(params - 10.E0).gt.tol))then + success = .false. + write(0,*) 'conv1d layer has wrong gradients' + end if + + !! check layer output handling + call conv_layer%get_output(params) + if(size(params) .ne. product(conv_layer%output_shape))then + success = .false. + write(0,*) 'conv1d layer has wrong number of outputs' + end if + call conv_layer%get_output(outputs) + if(any(shape(outputs) .ne. [product(conv_layer%output_shape), 1]))then + success = .false. + write(0,*) 'conv1d layer has wrong number of outputs' + end if + end select + + +!!!----------------------------------------------------------------------------- +!!! check expected initialisation of kernel size and stride +!!!----------------------------------------------------------------------------- + conv_layer = conv1d_layer_type( & + kernel_size = [2], & + stride = [2] & + ) + select type(conv_layer) + type is (conv1d_layer_type) + if(any(conv_layer%knl .ne. 2))then + success = .false. + write(0,*) 'conv1d layer has wrong pool size' + end if + if(any(conv_layer%stp .ne. 2))then + success = .false. + write(0,*) 'conv1d layer has wrong stride size' + end if + end select + + +!!!----------------------------------------------------------------------------- +!!! check expected initialisation of kernel size and stride +!!!----------------------------------------------------------------------------- + conv_layer = conv1d_layer_type( & + kernel_size = [4], & + stride = [4] & + ) + select type(conv_layer) + type is (conv1d_layer_type) + if(any(conv_layer%knl .ne. 4))then + success = .false. + write(0,*) 'conv1d layer has wrong pool size' + end if + if(any(conv_layer%stp .ne. 4))then + success = .false. + write(0,*) 'conv1d layer has wrong stride size' + end if + end select + + +!!!----------------------------------------------------------------------------- +!!! check layer operations +!!!----------------------------------------------------------------------------- + conv_layer1 = conv1d_layer_type( & + num_filters = 1, & + kernel_size = 3, & + activation_function = 'sigmoid' & + ) + call conv_layer1%init(input_layer%input_shape, batch_size=1) + conv_layer2 = conv1d_layer_type( & + num_filters = 1, & + kernel_size = 3, & + activation_function = 'sigmoid' & + ) + call conv_layer2%init(input_layer%input_shape, batch_size=1) + select type(conv_layer1) + type is(conv1d_layer_type) + select type(conv_layer2) + type is(conv1d_layer_type) + conv_layer = conv_layer1 + conv_layer2 + select type(conv_layer) + type is(conv1d_layer_type) + !! check layer addition + call compare_conv1d_layers(& + conv_layer, conv_layer1, success, conv_layer2) + + !! check layer reduction + conv_layer = conv_layer1 + call conv_layer%reduce(conv_layer2) + call compare_conv1d_layers(& + conv_layer, conv_layer1, success, conv_layer2) + + !! check layer merge + conv_layer = conv_layer1 + call conv_layer%merge(conv_layer2) + call compare_conv1d_layers(& + conv_layer, conv_layer1, success, conv_layer2) + class default + success = .false. + write(0,*) 'conv1d layer has wrong type' + end select + class default + success = .false. + write(0,*) 'conv1d layer has wrong type' + end select + class default + success = .false. + write(0,*) 'conv1d layer has wrong type' + end select + + +!!!----------------------------------------------------------------------------- +!!! Check for any failed tests +!!!----------------------------------------------------------------------------- + write(*,*) "----------------------------------------" + if(success)then + write(*,*) 'test_conv1d_layer passed all tests' + else + write(0,*) 'test_conv1d_layer failed one or more tests' + stop 1 + end if + + +contains + +!!!----------------------------------------------------------------------------- +!!! compare two or three layers +!!!----------------------------------------------------------------------------- + subroutine compare_conv1d_layers(layer1, layer2, success, layer3) + type(conv1d_layer_type), intent(in) :: layer1, layer2 + logical, intent(inout) :: success + type(conv1d_layer_type), optional, intent(in) :: layer3 + + if(all(layer1%knl .ne. layer2%knl))then + success = .false. + write(0,*) 'conv1d layer has wrong kernel_size' + end if + if(layer1%num_filters .ne. layer2%num_filters)then + success = .false. + write(0,*) 'conv1d layer has wrong num_filters' + end if + if(layer1%transfer%name .ne. 'sigmoid')then + success = .false. + write(0,*) 'conv1d layer has wrong transfer: '//& + layer1%transfer%name + end if + if(present(layer3))then + if(any(abs(layer1%dw-layer2%dw-layer3%dw).gt.tol))then + success = .false. + write(0,*) 'conv1d layer has wrong weights' + end if + if(any(abs(layer1%db-layer2%db-layer3%db).gt.tol))then + success = .false. + write(0,*) 'conv1d layer has wrong weights' + end if + end if + + end subroutine compare_conv1d_layers + +end program test_conv1d_layer \ No newline at end of file diff --git a/test/test_conv1d_network.f90 b/test/test_conv1d_network.f90 new file mode 100644 index 0000000..8917922 --- /dev/null +++ b/test/test_conv1d_network.f90 @@ -0,0 +1,127 @@ +program test_conv1d_network + use athena, only: & + network_type, & + conv1d_layer_type, & + base_optimiser_type + implicit none + + type(network_type) :: network + + integer :: i + integer, parameter :: num_channels = 3 + integer, parameter :: kernel_size = 3 + integer, parameter :: num_filters1 = 4 + integer, parameter :: num_filters2 = 8 + integer, parameter :: width = 7 + + real, allocatable, dimension(:,:) :: output_reshaped + real, allocatable, dimension(:,:,:) :: input_data, output, gradients_weight + real, allocatable, dimension(:) :: gradients, gradients_bias + logical :: success = .true. + + +!!!----------------------------------------------------------------------------- +!!! set up network +!!!----------------------------------------------------------------------------- + call network%add(conv1d_layer_type( & + input_shape=[width, width, num_channels], & + num_filters = num_filters1, & + kernel_size = kernel_size, & + kernel_initialiser = "ones", & + activation_function = "linear" & + )) + call network%add(conv1d_layer_type( & + num_filters = num_filters2, & + kernel_size = kernel_size, & + kernel_initialiser = "ones", & + activation_function = "linear" & + )) + + call network%compile( & + optimiser = base_optimiser_type(learning_rate=1.0), & + loss_method="mse", metrics=["loss"], verbose=1, & + batch_size=1) + + if(network%num_layers.ne.3)then + success = .false. + write(*,*) "conv1d network should have 3 layers" + end if + + call network%set_batch_size(1) + allocate(input_data(width, num_channels, 1)) + input_data = 0.0 + + call network%forward(input_data) + call network%model(3)%layer%get_output(output) + + if(any(shape(output).ne.[width-4,num_filters2,1]))then + success = .false. + write(*,*) "conv1d network output shape should be [28,28,32]" + end if + + +!!!----------------------------------------------------------------------------- +!!! check gradients +!!!----------------------------------------------------------------------------- + output = 0.E0 + output(:(width-4)/2,:,:) = 1.E0 + input_data = 0.E0 + input_data(:(width)/2,:,:) = 1.E0 + call network%forward(input_data) + output_reshaped = reshape(output, [(width-4)*num_filters2,1]) + call network%backward(output_reshaped) + select type(current => network%model(3)%layer) + type is(conv1d_layer_type) + gradients = current%get_gradients() + gradients_weight = & + reshape(& + gradients(:kernel_size*num_filters1*num_filters2), & + [kernel_size,num_filters1,num_filters2]) + gradients_bias = & + gradients(kernel_size*num_filters1*num_filters2+1:) + if(size(gradients).ne.& + (kernel_size*num_filters1 + 1) * num_filters2)then + success = .false. + write(*,*) "conv1d network gradients size should be ", & + ( kernel_size * num_filters1 + 1 ) * num_filters2 + end if + if(any(abs(gradients).lt.1.E-6))then + success = .false. + write(*,*) "conv1d network gradients should not be zero" + end if + if(any(abs(gradients_weight(1,:,:) - & + gradients_weight(1,1,1)).gt.1.E-6))then + success = .false. + write(*,*) "conv1d network gradients first column should be equivalent" + end if + if(any(abs(gradients_weight(2,:,:) - & + gradients_weight(2,1,1)).gt.1.E-6))then + success = .false. + write(*,*) "conv1d network gradients second column should be equivalent" + end if + if(any(abs(gradients_weight(3,:,:) - & + gradients_weight(3,1,1)).gt.1.E-6))then + success = .false. + write(*,*) "conv1d network gradients third column should be equivalent" + end if + if(any(abs(gradients_bias(:)-gradients_bias(1)).gt.1.E-6))then + success = .false. + write(*,*) "conv1d network gradients bias should be equivalent" + end if + class default + success = .false. + write(*,*) "conv1d network layer should be conv1d_layer_type" + end select + +!!!----------------------------------------------------------------------------- +!!! check for any failed tests +!!!----------------------------------------------------------------------------- + write(*,*) "----------------------------------------" + if(success)then + write(*,*) 'test_conv1d_network passed all tests' + else + write(0,*) 'test_conv1d_network failed one or more tests' + stop 1 + end if + +end program test_conv1d_network \ No newline at end of file diff --git a/test/test_conv2d_layer.f90 b/test/test_conv2d_layer.f90 index 1055782..6b3d53f 100644 --- a/test/test_conv2d_layer.f90 +++ b/test/test_conv2d_layer.f90 @@ -163,7 +163,7 @@ program test_conv2d_layer write(0,*) 'conv2d layer has wrong number of outputs' end if call conv_layer%get_output(outputs) - if(any(shape(outputs) .ne. conv_layer%output_shape))then + if(any(shape(outputs) .ne. [product(conv_layer%output_shape), 1]))then success = .false. write(0,*) 'conv2d layer has wrong number of outputs' end if diff --git a/test/test_conv3d_layer.f90 b/test/test_conv3d_layer.f90 index 8e653f7..72f61a2 100644 --- a/test/test_conv3d_layer.f90 +++ b/test/test_conv3d_layer.f90 @@ -163,7 +163,7 @@ program test_conv3d_layer write(0,*) 'conv3d layer has wrong number of outputs' end if call conv_layer%get_output(outputs) - if(any(shape(outputs) .ne. conv_layer%output_shape))then + if(any(shape(outputs) .ne. [product(conv_layer%output_shape), 1]))then success = .false. write(0,*) 'conv3d layer has wrong number of outputs' end if diff --git a/test/test_dropblock2d_layer.f90 b/test/test_dropblock2d_layer.f90 index 1566f65..82ba80c 100644 --- a/test/test_dropblock2d_layer.f90 +++ b/test/test_dropblock2d_layer.f90 @@ -15,7 +15,7 @@ program test_dropblock2d_layer integer :: i, j, output_width real, parameter :: max_value = 3.0 - integer :: seed_size = 1 + integer :: seed_size integer, allocatable, dimension(:) :: seed diff --git a/test/test_dropblock3d_layer.f90 b/test/test_dropblock3d_layer.f90 index 9e8b8af..fcc48b6 100644 --- a/test/test_dropblock3d_layer.f90 +++ b/test/test_dropblock3d_layer.f90 @@ -15,7 +15,7 @@ program test_dropblock3d_layer integer :: i, j, output_width real, parameter :: max_value = 3.0 - integer :: seed_size = 1 + integer :: seed_size integer, allocatable, dimension(:) :: seed @@ -144,4 +144,4 @@ program test_dropblock3d_layer stop 1 end if -end program test_dropblock3d_layer \ No newline at end of file +end program test_dropblock3d_layer diff --git a/test/test_dropout_layer.f90 b/test/test_dropout_layer.f90 index d6df2fe..abebba7 100644 --- a/test/test_dropout_layer.f90 +++ b/test/test_dropout_layer.f90 @@ -14,7 +14,7 @@ program test_dropout_layer integer :: i, j, output_width real, parameter :: max_value = 3.0 - integer :: seed_size = 1 + integer :: seed_size integer, allocatable, dimension(:) :: seed !! Initialize random number generator with a seed diff --git a/test/test_flatten1d_layer.f90 b/test/test_flatten1d_layer.f90 index d56c052..48eab35 100644 --- a/test/test_flatten1d_layer.f90 +++ b/test/test_flatten1d_layer.f90 @@ -12,7 +12,7 @@ program test_flatten1d_layer integer :: i, j, output_width - integer :: seed_size = 1 + integer :: seed_size integer, allocatable, dimension(:) :: seed diff --git a/test/test_flatten2d_layer.f90 b/test/test_flatten2d_layer.f90 index cd8f682..d4f1421 100644 --- a/test/test_flatten2d_layer.f90 +++ b/test/test_flatten2d_layer.f90 @@ -12,7 +12,7 @@ program test_flatten2d_layer integer :: i, j, output_width - integer :: seed_size = 1 + integer :: seed_size integer, allocatable, dimension(:) :: seed diff --git a/test/test_flatten3d_layer.f90 b/test/test_flatten3d_layer.f90 index 3458697..f9be24a 100644 --- a/test/test_flatten3d_layer.f90 +++ b/test/test_flatten3d_layer.f90 @@ -12,7 +12,7 @@ program test_flatten3d_layer integer :: i, j, output_width - integer :: seed_size = 1 + integer :: seed_size integer, allocatable, dimension(:) :: seed diff --git a/test/test_flatten4d_layer.f90 b/test/test_flatten4d_layer.f90 index a0fd0d8..0868c33 100644 --- a/test/test_flatten4d_layer.f90 +++ b/test/test_flatten4d_layer.f90 @@ -12,7 +12,7 @@ program test_flatten4d_layer integer :: i, j, output_width - integer :: seed_size = 1 + integer :: seed_size integer, allocatable, dimension(:) :: seed diff --git a/test/test_initialisers.f90 b/test/test_initialisers.f90 index 958d7fc..c8c85b3 100644 --- a/test/test_initialisers.f90 +++ b/test/test_initialisers.f90 @@ -8,7 +8,7 @@ program test_initialisers use initialiser, only: initialiser_setup, get_default_initialiser implicit none - class(initialiser_type), allocatable :: initialiser + class(initialiser_type), allocatable :: initialiser_var class(base_layer_type), allocatable :: full_layer, conv2d_layer, conv3d_layer logical :: success = .true. @@ -64,13 +64,13 @@ program test_initialisers !!! check initialisers work as expected for each rank !!!----------------------------------------------------------------------------- do i = 1, size(initialiser_names) - if(allocated(initialiser)) deallocate(initialiser) - allocate(initialiser, source=initialiser_setup(initialiser_names(i))) + if(allocated(initialiser_var)) deallocate(initialiser_var) + allocate(initialiser_var, source=initialiser_setup(initialiser_names(i))) if(.not.trim(initialiser_names(i)).eq."ident") & - call initialiser%initialise(input_0d, fan_in = 1, fan_out = 1) - call initialiser%initialise(input_1d, fan_in = 1, fan_out = 1) - call initialiser%initialise(input_3d, fan_in = 1, fan_out = 1) - call initialiser%initialise(input_6d, fan_in = 1, fan_out = 1) + call initialiser_var%initialise(input_0d, fan_in = 1, fan_out = 1) + call initialiser_var%initialise(input_1d, fan_in = 1, fan_out = 1) + call initialiser_var%initialise(input_3d, fan_in = 1, fan_out = 1) + call initialiser_var%initialise(input_6d, fan_in = 1, fan_out = 1) !! check for rank 2 data diff --git a/test/test_input_layer.f90 b/test/test_input_layer.f90 index b5435d7..6ffda8e 100644 --- a/test/test_input_layer.f90 +++ b/test/test_input_layer.f90 @@ -1,6 +1,7 @@ program test_input_layer use athena, only: & input1d_layer_type, & + input2d_layer_type, & input3d_layer_type, & input4d_layer_type, & base_layer_type @@ -70,14 +71,75 @@ program test_input_layer call input_layer%set(input_5d) call input_layer%forward(input_5d) if(any(abs(input_layer%output-5.E0).gt.1.E-6))then - write(0,*) 'input1d_layer forward 4d failed' + write(0,*) 'input1d_layer forward 5d failed' write(*,*) input_layer%output success = .false. end if call input_layer%set(input_6d) call input_layer%forward(input_6d) if(any(abs(input_layer%output-6.E0).gt.1.E-6))then - write(0,*) 'input1d_layer forward 4d failed' + write(0,*) 'input1d_layer forward 6d failed' + write(*,*) input_layer%output + success = .false. + end if + call input_layer%backward([1.E0], [1.E0]) + end select + + +!!!----------------------------------------------------------------------------- +!!! set up input2d layer +!!!----------------------------------------------------------------------------- + deallocate(input_layer) + input_layer = input2d_layer_type( & + input_shape=shape(input_3d), & + batch_size=batch_size) + if(input_layer%batch_size.ne.batch_size)then + write(0,*) 'input2d_layer batch_size failed' + success = .false. + end if + select type(input_layer) + type is (input2d_layer_type) + call input_layer%set(input_1d) + call input_layer%forward(input_1d) + call input_layer%get_output(output_1d) + if(any(abs(output_1d-1.E0).gt.1.E-6))then + write(0,*) 'input2d_layer forward 1d failed' + write(*,*) input_layer%output + success = .false. + end if + call input_layer%set(input_2d) + call input_layer%forward(input_2d) + call input_layer%get_output(output_2d) + if(any(abs(output_2d-2.E0).gt.1.E-6))then + write(0,*) 'input2d_layer forward 2d failed' + write(*,*) input_layer%output + success = .false. + end if + call input_layer%set(input_3d) + call input_layer%forward(input_3d) + if(any(abs(input_layer%output-3.E0).gt.1.E-6))then + write(0,*) 'input2d_layer forward 3d failed' + write(*,*) input_layer%output + success = .false. + end if + call input_layer%set(input_4d) + call input_layer%forward(input_4d) + if(any(abs(input_layer%output-4.E0).gt.1.E-6))then + write(0,*) 'input2d_layer forward 4d failed' + write(*,*) input_layer%output + success = .false. + end if + call input_layer%set(input_5d) + call input_layer%forward(input_5d) + if(any(abs(input_layer%output-5.E0).gt.1.E-6))then + write(0,*) 'input2d_layer forward 5d failed' + write(*,*) input_layer%output + success = .false. + end if + call input_layer%set(input_6d) + call input_layer%forward(input_6d) + if(any(abs(input_layer%output-6.E0).gt.1.E-6))then + write(0,*) 'input2d_layer forward 6d failed' write(*,*) input_layer%output success = .false. end if @@ -131,14 +193,14 @@ program test_input_layer call input_layer%set(input_5d) call input_layer%forward(input_5d) if(any(abs(input_layer%output-5.E0).gt.1.E-6))then - write(0,*) 'input3d_layer forward 4d failed' + write(0,*) 'input3d_layer forward 5d failed' write(*,*) input_layer%output success = .false. end if call input_layer%set(input_6d) call input_layer%forward(input_6d) if(any(abs(input_layer%output-6.E0).gt.1.E-6))then - write(0,*) 'input3d_layer forward 4d failed' + write(0,*) 'input3d_layer forward 6d failed' write(*,*) input_layer%output success = .false. end if @@ -192,14 +254,14 @@ program test_input_layer call input_layer%set(input_5d) call input_layer%forward(input_5d) if(any(abs(input_layer%output-5.E0).gt.1.E-6))then - write(0,*) 'input4d_layer forward 4d failed' + write(0,*) 'input4d_layer forward 5d failed' write(*,*) input_layer%output success = .false. end if call input_layer%set(input_6d) call input_layer%forward(input_6d) if(any(abs(input_layer%output-6.E0).gt.1.E-6))then - write(0,*) 'input4d_layer forward 4d failed' + write(0,*) 'input4d_layer forward 6d failed' write(*,*) input_layer%output success = .false. end if diff --git a/test/test_maxpool1d_layer.f90 b/test/test_maxpool1d_layer.f90 new file mode 100644 index 0000000..57dd22b --- /dev/null +++ b/test/test_maxpool1d_layer.f90 @@ -0,0 +1,204 @@ +program test_maxpool1d_layer + use athena, only: & + maxpool1d_layer_type, & + base_layer_type, & + learnable_layer_type + implicit none + + class(base_layer_type), allocatable :: pool_layer + integer, parameter :: num_channels = 3, pool = 3, stride = 2, width = 18 + real, allocatable, dimension(:,:,:) :: input_data, output, gradient + real, allocatable, dimension(:) :: output_1d + real, allocatable, dimension(:,:) :: output_2d + real, parameter :: tol = 1.E-7 + logical :: success = .true. + + integer :: i, j, output_width, max_loc + integer :: num_windows + real, parameter :: max_value = 3.0 + + +!!!----------------------------------------------------------------------------- +!!! set up layer +!!!----------------------------------------------------------------------------- + pool_layer = maxpool1d_layer_type( & + pool_size = pool, & + stride = stride & + ) + + !! check layer name + if(.not. pool_layer%name .eq. 'maxpool1d')then + success = .false. + write(0,*) 'maxpool1d layer has wrong name' + end if + + !! check layer type + select type(pool_layer) + type is(maxpool1d_layer_type) + !! check pool size + if(any(pool_layer%pool .ne. pool))then + success = .false. + write(0,*) 'maxpool1d layer has wrong pool size' + end if + + !! check stride size + if(any(pool_layer%strd .ne. stride))then + success = .false. + write(0,*) 'maxpool1d layer has wrong stride size' + end if + + !! check input shape allocated + if(allocated(pool_layer%input_shape))then + success = .false. + write(0,*) 'maxpool1d layer shape should not be allocated yet' + end if + class default + success = .false. + write(0,*) 'maxpool1d layer has wrong type' + end select + + +!!!----------------------------------------------------------------------------- +!!! check layer input and output shape based on input layer +!!!----------------------------------------------------------------------------- + !! initialise width and output width + output_width = floor( (width - pool)/real(stride)) + 1 + max_loc = width / 2 + mod(width, 2) + + !! initialise sample input + allocate(input_data(width, num_channels, 1), source = 0.0) + input_data(max_loc, :, 1) = max_value + pool_layer = maxpool1d_layer_type( & + pool_size = pool, & + stride = stride & + ) + + !! check layer input and output shape based on input data + call pool_layer%init(shape(input_data(:,:,1)), batch_size=1) + select type(pool_layer) + type is(maxpool1d_layer_type) + if(any(pool_layer%input_shape .ne. [width,num_channels]))then + success = .false. + write(0,*) 'maxpool1d layer has wrong input_shape' + end if + if(any(pool_layer%output_shape .ne. & + [output_width,num_channels]))then + success = .false. + write(0,*) 'maxpool1d layer has wrong output_shape', & + pool_layer%output_shape + write(0,*) 'expected', [output_width,num_channels] + end if + end select + + +!!!----------------------------------------------------------------------------- +!!! test forward pass and check expected output +!!! use existing layer +!!!----------------------------------------------------------------------------- + !! run forward pass + call pool_layer%forward(input_data) + call pool_layer%get_output(output) + + !! check outputs have expected value + do i = 1, output_width + if( max_loc .ge. (i-1)*stride + 1 .and. & + max_loc .le. (i-1)*stride + pool )then + if(output(i, 1, 1) .ne. max_value)then + success = .false. + write(*,*) 'maxpool1d layer forward pass failed' + end if + else if(output(i, 1, 1) .ne. 0.0) then + success = .false. + write(*,*) 'maxpool1d layer forward pass failed' + end if + end do + + +!!!----------------------------------------------------------------------------- +!!! check output request using rank 1 and rank 2 arrays is consistent +!!!----------------------------------------------------------------------------- + call pool_layer%get_output(output_1d) + call pool_layer%get_output(output_2d) + if(any(abs(output_1d - reshape(output_2d, [size(output_2d)])) .gt. 1.E-6))then + success = .false. + write(0,*) 'output_1d and output_2d are not consistent' + end if + + +!!!----------------------------------------------------------------------------- +!!! test backward pass and check expected output +!!!----------------------------------------------------------------------------- + !! run backward pass + allocate(gradient, source = output) + call pool_layer%backward(input_data, gradient) + + !! check gradient has expected value + select type(current => pool_layer) + type is(maxpool1d_layer_type) + do i = 1, width + num_windows = pool - stride + 1 - mod((stride+1)*(i-1),2) + if(i.eq.maxloc(input_data(:,1,1),dim=1))then + if(current%di(i, 1, 1) .ne. maxval(output)*num_windows)then + success = .false. + write(*,*) num_windows + write(*,*) 'maxpool1d layer backward pass failed' + end if + else + if(current%di(i, 1, 1) .ne. 0.0) then + success = .false. + write(*,*) 'maxpool1d layer backward pass failed' + end if + end if + end do + end select + + +!!!----------------------------------------------------------------------------- +!!! check expected initialisation of pool and stride +!!!----------------------------------------------------------------------------- + pool_layer = maxpool1d_layer_type( & + pool_size = [2], & + stride = [2] & + ) + select type(pool_layer) + type is (maxpool1d_layer_type) + if(any(pool_layer%pool .ne. [2]))then + success = .false. + write(0,*) 'maxpool1d layer has wrong pool size' + end if + if(any(pool_layer%strd .ne. [2]))then + success = .false. + write(0,*) 'maxpool1d layer has wrong stride size' + end if + end select + + !! check expected initialisation of pool and stride + pool_layer = maxpool1d_layer_type( & + pool_size = [4], & + stride = [4] & + ) + select type(pool_layer) + type is (maxpool1d_layer_type) + if(any(pool_layer%pool .ne. 4))then + success = .false. + write(0,*) 'maxpool1d layer has wrong pool size' + end if + if(any(pool_layer%strd .ne. 4))then + success = .false. + write(0,*) 'maxpool1d layer has wrong stride size' + end if + end select + + +!!!----------------------------------------------------------------------------- +!!! check for any failed tests +!!!----------------------------------------------------------------------------- + write(*,*) "----------------------------------------" + if(success)then + write(*,*) 'test_maxpool1d_layer passed all tests' + else + write(0,*) 'test_maxpool1d_layer failed one or more tests' + stop 1 + end if + +end program test_maxpool1d_layer \ No newline at end of file diff --git a/test/test_network.f90 b/test/test_network.f90 index df35083..65d578b 100644 --- a/test/test_network.f90 +++ b/test/test_network.f90 @@ -40,9 +40,8 @@ program test_network !!!----------------------------------------------------------------------------- !!! Initialize random number generator with a seed !!!----------------------------------------------------------------------------- - seed_size = 8 call random_seed(size=seed_size) - seed = [1,1,1,1,1,1,1,1] + allocate(seed(seed_size), source=1) call random_seed(put=seed) diff --git a/test/test_optimiser.f90 b/test/test_optimiser.f90 index 24ccb14..cc9d15a 100644 --- a/test/test_optimiser.f90 +++ b/test/test_optimiser.f90 @@ -4,7 +4,7 @@ program test_optimiser use learning_rate_decay implicit none - class(base_optimiser_type), allocatable :: optimiser + class(base_optimiser_type), allocatable :: optimiser_var integer, parameter :: num_params = 10 @@ -15,36 +15,36 @@ program test_optimiser !!!----------------------------------------------------------------------------- -!!! test empty base optimiser +!!! test empty base optimiser_var !!!----------------------------------------------------------------------------- - allocate(optimiser, source=base_optimiser_type()) + allocate(optimiser_var, source=base_optimiser_type()) !!!----------------------------------------------------------------------------- !!! test base optimiser !!!----------------------------------------------------------------------------- - deallocate(optimiser) - allocate(optimiser, source=base_optimiser_type( & + deallocate(optimiser_var) + allocate(optimiser_var, source=base_optimiser_type( & learning_rate = 0.1, & num_params = num_params, & lr_decay = base_lr_decay_type(), & regulariser=l1_regulariser_type() & )) - select type(optimiser) + select type(optimiser_var) type is (base_optimiser_type) write(*,*) "base_optimiser_type" class default write(0,*) "Failed to allocate base optimiser" success = .false. end select - call optimiser%minimise(param, gradient) + call optimiser_var%minimise(param, gradient) !!!----------------------------------------------------------------------------- !!! test sgd optimiser !!!----------------------------------------------------------------------------- - deallocate(optimiser) - allocate(optimiser, source=sgd_optimiser_type( & + deallocate(optimiser_var) + allocate(optimiser_var, source=sgd_optimiser_type( & learning_rate = 0.1, & num_params = num_params, & momentum = 0.9, & @@ -52,21 +52,21 @@ program test_optimiser lr_decay = base_lr_decay_type(), & regulariser=l1_regulariser_type() & )) - select type(optimiser) + select type(optimiser_var) type is (sgd_optimiser_type) write(*,*) "sgd_optimiser_type" class default write(0,*) "Failed to allocate sgd optimiser" success = .false. end select - call optimiser%minimise(param, gradient) + call optimiser_var%minimise(param, gradient) !!!----------------------------------------------------------------------------- !!! test rmsprop optimiser !!!----------------------------------------------------------------------------- - deallocate(optimiser) - allocate(optimiser, source=rmsprop_optimiser_type( & + deallocate(optimiser_var) + allocate(optimiser_var, source=rmsprop_optimiser_type( & learning_rate = 0.1, & num_params = num_params, & beta = 0.9, & @@ -74,42 +74,42 @@ program test_optimiser lr_decay = base_lr_decay_type(), & regulariser=l1_regulariser_type() & )) - select type(optimiser) + select type(optimiser_var) type is (rmsprop_optimiser_type) write(*,*) "rmsprop_optimiser_type" class default write(0,*) "Failed to allocate rmsprop optimiser" success = .false. end select - call optimiser%minimise(param, gradient) + call optimiser_var%minimise(param, gradient) !!!----------------------------------------------------------------------------- !!! test adagrad optimiser !!!----------------------------------------------------------------------------- - deallocate(optimiser) - allocate(optimiser, source=adagrad_optimiser_type( & + deallocate(optimiser_var) + allocate(optimiser_var, source=adagrad_optimiser_type( & learning_rate = 0.1, & num_params = num_params, & epsilon = 1e-8, & lr_decay = base_lr_decay_type(), & regulariser=l1_regulariser_type() & )) - select type(optimiser) + select type(optimiser_var) type is (adagrad_optimiser_type) write(*,*) "adagrad_optimiser_type" class default write(0,*) "Failed to allocate adagrad optimiser" success = .false. end select - call optimiser%minimise(param, gradient) + call optimiser_var%minimise(param, gradient) !!!----------------------------------------------------------------------------- !!! test adam optimiser with l1 regulariser !!!----------------------------------------------------------------------------- - deallocate(optimiser) - allocate(optimiser, source=adam_optimiser_type( & + deallocate(optimiser_var) + allocate(optimiser_var, source=adam_optimiser_type( & learning_rate = 0.1, & num_params = num_params, & beta1 = 0.9, & @@ -118,21 +118,21 @@ program test_optimiser lr_decay = base_lr_decay_type(), & regulariser=l1_regulariser_type() & )) - select type(optimiser) + select type(optimiser_var) type is (adam_optimiser_type) write(*,*) "adam_optimiser_type" class default write(0,*) "Failed to allocate adam optimiser" success = .false. end select - call optimiser%minimise(param, gradient) + call optimiser_var%minimise(param, gradient) !!!----------------------------------------------------------------------------- !!! test adam optimiser with l2 regulariser !!!----------------------------------------------------------------------------- - deallocate(optimiser) - allocate(optimiser, source=adam_optimiser_type( & + deallocate(optimiser_var) + allocate(optimiser_var, source=adam_optimiser_type( & learning_rate = 0.1, & num_params = num_params, & beta1 = 0.9, & @@ -141,21 +141,21 @@ program test_optimiser lr_decay = base_lr_decay_type(), & regulariser=l2_regulariser_type() & )) - select type(optimiser) + select type(optimiser_var) type is (adam_optimiser_type) write(*,*) "adam_optimiser_type" class default write(0,*) "Failed to allocate adam optimiser" success = .false. end select - call optimiser%minimise(param, gradient) + call optimiser_var%minimise(param, gradient) !!!----------------------------------------------------------------------------- !!! test adam optimiser with l2 decoupled regulariser !!!----------------------------------------------------------------------------- - deallocate(optimiser) - allocate(optimiser, source=adam_optimiser_type( & + deallocate(optimiser_var) + allocate(optimiser_var, source=adam_optimiser_type( & learning_rate = 0.1, & num_params = num_params, & beta1 = 0.9, & @@ -164,14 +164,14 @@ program test_optimiser lr_decay = base_lr_decay_type(), & regulariser=l2_regulariser_type(decoupled=.false.) & )) - select type(optimiser) + select type(optimiser_var) type is (adam_optimiser_type) write(*,*) "adam_optimiser_type" class default write(0,*) "Failed to allocate adam optimiser" success = .false. end select - call optimiser%minimise(param, gradient) + call optimiser_var%minimise(param, gradient) !!!-----------------------------------------------------------------------------