Enabling double precision real number

KerasWeightsProcessing/convert_weights.py

@@ -14,7 +14,7 @@
 from keras import optimizers
 
 INPUT = ['input']
-ACTIVATIONS = ['relu', 'linear', 'leakyrelu', 'sigmoid']
+ACTIVATIONS = ['swish', 'elu', 'exponential', 'relu', 'linear', 'leakyrelu', 'softplus', 'sigmoid', 'tanh']
 SUPPORTED_LAYERS = ['dense', 'dropout', 'batchnormalization'] + ACTIVATIONS + INPUT
 
 def txt_to_h5(weights_file_name, output_file_name=''):
@@ -70,8 +70,9 @@ def txt_to_h5(weights_file_name, output_file_name=''):
                 elif layer_type == 'batchnormalization':
                     batchnorm_count += 4
                     x = BatchNormalization(name='batch_normalization_{}'.format(batchnorm_count // 4))(x)
-                elif layer_type == 'linear':
-                    x = Activation('linear')(x)
+
+                elif layer_type == 'swish':
+                    x = Activation('swish')(x)
 
             elif not layer_type.isalpha():
                 if lr == False:
@@ -156,7 +157,10 @@ def h5_to_txt(weights_file_name, output_file_name=''):
         keras_version = weights_file.attrs['keras_version']
 
         if 'training_config' in weights_file.attrs:
-            training_config = weights_file.attrs['training_config'].decode('utf-8')
+            try:
+                training_config = weights_file.attrs['training_config'].decode('utf-8')
+            except AttributeError: 
+                training_config = weights_file.attrs['training_config']
             training_config = training_config.replace('true','True')
             training_config = training_config.replace('false','False')
             training_config = training_config.replace('null','None')
@@ -169,7 +173,10 @@ def h5_to_txt(weights_file_name, output_file_name=''):
             learning_rate = 0.001
 
         # Decode using the utf-8 encoding; change values for eval
-        model_config = weights_file.attrs['model_config'].decode('utf-8')
+        try:
+            model_config = weights_file.attrs['model_config'].decode('utf-8')
+        except AttributeError:
+            model_config = weights_file.attrs['model_config']
         model_config = model_config.replace('true','True')
         model_config = model_config.replace('false','False')
         model_config = model_config.replace('null','None')
@@ -197,7 +204,7 @@ def h5_to_txt(weights_file_name, output_file_name=''):
                 input_layers = model_config['config'].get('input_layers',[])
                 input_names = [layer[0] for layer in input_layers]
 
-        else:
+        else: # 'Sequential'
             layer_config = model_config['config']['layers']
 
         for idx,layer in enumerate(layer_config):
@@ -242,12 +249,21 @@ def h5_to_txt(weights_file_name, output_file_name=''):
                     )
                 )
                 # add information about the activation
-                layer_info.append(
-                    info_str.format(
-                        name = activation,
-                        info = 0
+                if (activation == 'elu'):
+                    layer_info.append(
+                        info_str.format(
+                            name = activation,
+                            info = 1
+                        )
                     )
-                )
+                else:
+                    layer_info.append(
+                        info_str.format(
+                            name = activation,
+                            info = 0
+                        )
+                    )
+
             elif class_name == 'batchnormalization':
                 # get beta, gamma, moving_mean, moving_variance from dictionary
                 for key in sorted(model_weights[name][name].keys()):
@@ -276,9 +292,18 @@ def h5_to_txt(weights_file_name, output_file_name=''):
 
             elif class_name in ACTIVATIONS:
                 # replace previous dense layer with the advanced activation function (LeakyReLU)
+
+                try:
+                    _tmp = layer['config']['alpha']
+                except KeyError:
+                    if class_name == 'relu':
+                        _tmp = 0 # alpha is currently not used in mod_activation
+                    else:
+                        raise KeyError("alpha")
+
                 layer_info[-1] = info_str.format(
                     name = class_name,
-                    info = layer['config']['alpha']
+                    info = _tmp
                 )
 
             # if there are multiple outputs, remove what was just added

README.md

@@ -51,7 +51,7 @@ Check out an example in the [getting started notebook](https://github.com/scient
 Get the code:
 
 ```
-git clone https://github.com/scientific-computing/FKB
+git clone https://github.com/sungdukyu/FKB64
 ```
 
 Dependencies:

build_steps.sh

@@ -2,7 +2,12 @@ rm -rf build
 mkdir build
 cd build
 
-FC=gfortran cmake .. -DSERIAL=1
+# perlmutter
+module load PrgEnv-gnu/8.3.3
+module load gcc/11.2.0
+
+# FC=gfortran cmake .. -DSERIAL=1
+FC=gfortran cmake .. -DSERIAL=1 -DREAL=64    # enable double precision real number (note that the default precision for real number is single. See src/lib/mod_kinds.F90)
 # FC='mpif90 -qopenmp' cmake .. -DSERIAL=1
 
 make

src/lib/mod_activation.F90

@@ -16,6 +16,10 @@ module mod_activation
   public :: tanhf, tanh_prime
   public :: linear, linear_prime
   public :: leaky_relu, leaky_relu_prime
+  public :: elu, elu_prime
+  public :: exponential, exponential_prime
+  public :: softplus, softplus_prime
+  public :: swish, swish_prime
 
   interface
     pure function activation_function(x, alpha)
@@ -28,6 +32,64 @@ end function activation_function
 
 contains
 
+  pure function swish(x, alpha) result(res)
+    !! swish (or silu) activation function.
+    !! beta is fixed (i.e., beta==1.).
+    real(rk), intent(in) :: x(:)
+    real(rk), intent(in) :: alpha
+    real(rk) :: res(size(x))
+    res = x * sigmoid(x, -999._rk)
+  end function swish
+
+  pure function swish_prime(x, alpha) result(res)
+    ! First derivative of swish activation function.
+    ! beta is fixed (i.e., beta==1.).
+    real(rk), intent(in) :: x(:)
+    real(rk), intent(in) :: alpha
+    real(rk) :: res(size(x))
+    res = sigmoid(x, -999._rk) + x * sigmoid_prime(x, -999._rk)
+  end function swish_prime
+
+  pure function elu(x, alpha) result(res)
+    !! Exponential Linear Unit (ELU) activation function.
+    real(rk), intent(in) :: x(:)
+    real(rk), intent(in) :: alpha
+    real(rk) :: res(size(x))
+    where (x > 0)
+      res = x
+    elsewhere
+      res = alpha * (exp(x) - 1)
+    end where
+  end function elu
+
+  pure function elu_prime(x, alpha) result(res)
+    ! First derivative of the Exponential Linear Unit (ELU) activation function.
+    real(rk), intent(in) :: x(:)
+    real(rk), intent(in) :: alpha
+    real(rk) :: res(size(x))
+    where (x > 0)
+      res = 1
+    elsewhere
+      res = alpha * exp(x)
+    end where
+  end function elu_prime
+
+  pure function exponential(x, alpha) result(res)
+    !! Exponential activation function.
+    real(rk), intent(in) :: x(:)
+    real(rk), intent(in) :: alpha
+    real(rk) :: res(size(x))
+    res = exp(x)
+  end function exponential
+
+  pure function exponential_prime(x, alpha) result(res)
+    !! First derivative of the exponential activation function.
+    real(rk), intent(in) :: x(:)
+    real(rk), intent(in) :: alpha
+    real(rk) :: res(size(x))
+    res = exp(x)
+  end function exponential_prime
+
   pure function gaussian(x, alpha) result(res)
     ! Gaussian activation function.
     real(rk), intent(in) :: x(:)
@@ -63,7 +125,7 @@ pure function leaky_relu_prime(x, alpha) result(res)
     where (0.3 * x > 0)
       res = 1
     elsewhere
-      res = 0
+      res = alpha
     end where
   end function leaky_relu_prime
 
@@ -103,12 +165,19 @@ pure function linear_prime(x, alpha) result(res)
     res = 1
   end function linear_prime
 
+  ! Balwinder's version (July 5th, 2022)
+  ! - addressing floating point overflow
   pure function sigmoid(x, alpha) result(res)
     ! Sigmoid activation function.
     real(rk), intent(in) :: x(:)
     real(rk), intent(in) :: alpha
-    real(rk) :: res(size(x))
-    res = 1 / (1 + exp(-x))
+    real(rk) :: res(size(x)),y(size(x))
+    real(rk), parameter :: exp_lim = log(tiny(x)) !precision based limiter
+    y(:) = x(:)
+    where(y < exp_lim)
+       y = exp_lim
+    end where
+    res = 1. / (1. + exp(-y))
   endfunction sigmoid
 
   pure function sigmoid_prime(x, alpha) result(res)
@@ -122,6 +191,22 @@ pure function sigmoid_prime(x, alpha) result(res)
     res = sigmoid(x, tmp_alpha) * (1 - sigmoid(x, tmp_alpha))
   end function sigmoid_prime
 
+  pure function softplus(x, alpha) result(res)
+    !! Softplus activation function.
+    real(rk), intent(in) :: x(:)
+    real(rk), intent(in) :: alpha
+    real(rk) :: res(size(x))
+    res = log(exp(x) + 1)
+  end function softplus
+
+  pure function softplus_prime(x, alpha) result(res)
+    ! First derivative of the Softplus activation function.
+    real(rk), intent(in) :: x(:)
+    real(rk), intent(in) :: alpha
+    real(rk) :: res(size(x))
+    res = 1 / (1 + exp(-x))
+  end function softplus_prime
+
   pure function step(x, alpha) result(res)
     ! Step activation function.
     real(rk), intent(in) :: x(:)

src/lib/mod_dense_layer.F90

@@ -59,6 +59,18 @@ type(Dense) function constructor(this_size, next_size, activation, alpha) result
 
     ! assign activation function
     select case(trim(activation))
+      case('swish')
+        layer % activation => swish
+        layer % activation_prime => swish_prime
+      case('silu')
+        layer % activation => swish
+        layer % activation_prime => swish_prime
+      case('elu')
+        layer % activation => elu
+        layer % activation_prime => elu_prime
+      case('exponential')
+        layer % activation => exponential
+        layer % activation_prime => exponential_prime
       case('gaussian')
         layer % activation => gaussian
         layer % activation_prime => gaussian_prime
@@ -71,6 +83,9 @@ type(Dense) function constructor(this_size, next_size, activation, alpha) result
       case('sigmoid')
         layer % activation => sigmoid
         layer % activation_prime => sigmoid_prime
+      case('softplus')
+        layer % activation => softplus
+        layer % activation_prime => softplus_prime
       case('step')
         layer % activation => step
         layer % activation_prime => step_prime

src/lib/mod_kinds.F90

@@ -11,6 +11,8 @@ module mod_kinds
   integer,parameter :: rk = real64
 #elif REAL128
   integer,parameter :: rk = real128
+#elif REAL32
+  integer,parameter :: rk = real32
 #else
   integer,parameter :: rk = real32
 #endif
@@ -19,6 +21,8 @@ module mod_kinds
 
 #ifdef INT64
   integer, parameter :: ik = int64
+#elif INT32
+  integer, parameter :: ik = int32
 #else
   integer, parameter :: ik = int32
 #endif

src/lib/mod_layer.F90

@@ -68,15 +68,15 @@ end subroutine layer_backward
 
   pure type(array1d) function array1d_constructor(length) result(a)
     ! Overloads the default type constructor.
-    integer, intent(in) :: length
+    integer(ik), intent(in) :: length
     allocate(a % array(length))
     a % array = 0
   end function array1d_constructor
 
 
   pure type(array2d) function array2d_constructor(dims) result(a)
     ! Overloads the default type constructor.
-    integer, intent(in) :: dims(2)
+    integer(ik), intent(in) :: dims(2)
     allocate(a % array(dims(1), dims(2)))
     a % array = 0
   end function array2d_constructor
@@ -106,7 +106,7 @@ pure subroutine dw_init(dw, dims)
     do n = 1, nm - 1
       dw(n) = array2d(dims(n:n+1))
     end do
-    dw(n) = array2d([dims(n), 1])
+    dw(n) = array2d([dims(n), 1_ik])
   end subroutine dw_init
 
 

src/lib/mod_network.F90

@@ -56,7 +56,7 @@ type(network_type) function net_constructor(layer_names, layer_info) result(net)
 
     call net % init(layer_names, layer_info)
 
-    call net % sync(1)
+    call net % sync(1_ik)
 
   end function net_constructor
 

src/tests/test_ensembles.F90

@@ -23,7 +23,7 @@ program test_ensembles
   call get_command_argument(1,args(1))
 
   ! build ensemble from members in specified directory
-  ensemble = ensemble_type(args(1), 0.0)
+  ensemble = ensemble_type(args(1), 0.0_rk)
 
   input = [1, 2, 3, 4, 5]