safegaurd Box-Muller normal random number generation against u=0.0 (#158)

* safegaurd Box-Muller normal random number generation against u=0.0 * Docstrings and formatting * Bump patch version --------- Co-authored-by: milancurcic <caomaco@gmail.com>

Author: dacarnazzola

CMakeLists.txt

@@ -56,7 +56,6 @@ add_library(neural
  src/nf/nf_parallel.f90
  src/nf/nf_parallel_submodule.f90
  src/nf/nf_random.f90
- src/nf/nf_random_submodule.f90
  src/nf/nf_reshape_layer.f90
  src/nf/nf_reshape_layer_submodule.f90
  src/nf/io/nf_io_binary.f90

fpm.toml

@@ -1,5 +1,5 @@
 name = "neural-fortran"
-version = "0.15.0"
+version = "0.15.1"
 license = "MIT"
 author = "Milan Curcic"
 maintainer = "milancurcic@hey.com"

src/nf/nf_conv2d_layer_submodule.f90

@@ -1,7 +1,7 @@
 submodule(nf_conv2d_layer) nf_conv2d_layer_submodule
 
  use nf_activation, only: activation_function
- use nf_random, only: randn
+ use nf_random, only: random_normal
 
  implicit none
 
@@ -40,9 +40,8 @@ module subroutine init(self, input_shape)
  self % kernel_size, self % kernel_size))
 
  ! Initialize the kernel with random values with a normal distribution.
- self % kernel = randn(self % filters, self % channels, &
- self % kernel_size, self % kernel_size) &
- / self % kernel_size**2 !TODO kernel_width * kernel_height
+ call random_normal(self % kernel)
+ self % kernel = self % kernel / self % kernel_size**2
 
  allocate(self % biases(self % filters))
  self % biases = 0

src/nf/nf_dense_layer_submodule.f90

@@ -2,7 +2,7 @@
 
  use nf_activation, only: activation_function
  use nf_base_layer, only: base_layer
- use nf_random, only: randn
+ use nf_random, only: random_normal
 
  implicit none
 
@@ -114,8 +114,8 @@ module subroutine init(self, input_shape)
  ! Weights are a 2-d array of shape previous layer size
  ! times this layer size.
  allocate(self % weights(self % input_size, self % output_size))
- self % weights = randn(self % input_size, self % output_size) &
-   / self % input_size
+ call random_normal(self % weights)
+ self % weights = self % weights / self % input_size
 
  ! Broadcast weights to all other images, if any.
  call co_broadcast(self % weights, 1)

src/nf/nf_random.f90

@@ -6,34 +6,21 @@ module nf_random
  implicit none
 
  private
- public :: randn
-
- interface randn
-
- module function randn_1d(i) result(r)
- !! Generates i random numbers with a normal distribution,
- !! using the Box-Muller method.
- implicit none
- integer, intent(in) :: i
- real :: r(i)
- end function randn_1d
-
- module function randn_2d(i, j) result(r)
- !! Generates i x j random numbers with a normal distribution,
- !! using the Box-Muller method.
- implicit none
- integer, intent(in) :: i, j
- real :: r(i,j)
- end function randn_2d
-
- module function randn_4d(i, j, k, l) result(r)
- !! Generates i x j x k x l random numbers with a normal distribution,
- !! using the Box-Muller method.
- implicit none
- integer, intent(in) :: i, j, k, l
- real :: r(i,j,k,l)
- end function randn_4d
-
- end interface randn
+ public :: random_normal
+
+ real, parameter :: pi = 4 * atan(1.d0)
+
+contains
+
+ impure elemental subroutine random_normal(x)
+ !! Sample random numbers from a normal distribution using a Box-Muller
+ !! formula.
+ real, intent(out) :: x
+ !! Scalar or array to be filled with random numbers
+ real :: u(2)
+ call random_number(u)
+ u(1) = 1 - u(1)
+ x = sqrt(- 2 * log(u(1))) * cos(2 * pi * u(2))
+ end subroutine random_normal
 
 end module nf_random