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 <[email protected]>

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 = "[email protected]"

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