From 9958bb781fd99a8f08dc981dcd997a2ad68189f4 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 13 Oct 2024 17:00:36 -0700
Subject: [PATCH 1/6] feat(trainable_network): extend inference_engine_t

---
 src/inference_engine/inference_engine_m_.f90 |  77 +++++--
 src/inference_engine/inference_engine_s.F90  | 205 ++++++++++++++++++-
 src/inference_engine/trainable_engine_s.F90  |   1 -
 src/inference_engine/trainable_network_m.f90 |  42 ++++
 src/inference_engine/trainable_network_s.f90 |  15 ++
 src/inference_engine/workspace_s.f90         | 117 +++++++++++
 src/inference_engine_m.f90                   |   1 +
 7 files changed, 434 insertions(+), 24 deletions(-)
 create mode 100644 src/inference_engine/trainable_network_m.f90
 create mode 100644 src/inference_engine/trainable_network_s.f90
 create mode 100644 src/inference_engine/workspace_s.f90

diff --git a/src/inference_engine/inference_engine_m_.f90 b/src/inference_engine/inference_engine_m_.f90
index daeeccf8e..a196fd6ee 100644
--- a/src/inference_engine/inference_engine_m_.f90
+++ b/src/inference_engine/inference_engine_m_.f90
@@ -7,6 +7,7 @@ module inference_engine_m_
   use kind_parameters_m, only : default_real, double_precision
   use julienne_m, only : file_t, string_t
   use metadata_m, only : metadata_t
+  use mini_batch_m, only : mini_batch_t
   use tensor_m, only : tensor_t
   use tensor_map_m, only : tensor_map_t
   implicit none
@@ -15,6 +16,7 @@ module inference_engine_m_
   public :: inference_engine_t
   public :: unmapped_engine_t
   public :: exchange_t
+  public :: workspace_t
 
   type inference_engine_t(k)
     !! Encapsulate the minimal information needed to perform inference
@@ -38,19 +40,33 @@ module inference_engine_m_
     generic :: skip                     => default_real_skip,                    double_precision_skip
     generic :: activation_function_name => default_real_activation_name,         double_precision_activation_name
     generic :: to_exchange              => default_real_to_exchange,             double_precision_to_exchange
-    procedure, private :: default_real_approximately_equal,     double_precision_approximately_equal
+    generic :: learn                    => default_real_learn
+    procedure, private, non_overridable :: default_real_approximately_equal,     double_precision_approximately_equal
     procedure, private, non_overridable :: default_real_infer,                   double_precision_infer
-    procedure, private :: default_real_to_json,                 double_precision_to_json
-    procedure, private :: default_real_map_to_input_range,      double_precision_map_to_input_range
-    procedure, private :: default_real_map_from_output_range,   double_precision_map_from_output_range
-    procedure, private :: default_real_num_hidden_layers,       double_precision_num_hidden_layers
-    procedure, private :: default_real_num_inputs,              double_precision_num_inputs
-    procedure, private :: default_real_num_outputs,             double_precision_num_outputs
-    procedure, private :: default_real_nodes_per_layer,         double_precision_nodes_per_layer
-    procedure, private :: default_real_assert_conformable_with, double_precision_assert_conformable_with
-    procedure, private :: default_real_skip,                    double_precision_skip
-    procedure, private :: default_real_activation_name,         double_precision_activation_name
-    procedure, private :: default_real_to_exchange,             double_precision_to_exchange
+    procedure, private, non_overridable :: default_real_learn
+    procedure, private, non_overridable :: default_real_to_json,                 double_precision_to_json
+    procedure, private, non_overridable :: default_real_map_to_input_range,      double_precision_map_to_input_range
+    procedure, private, non_overridable :: default_real_map_from_output_range,   double_precision_map_from_output_range
+    procedure, private, non_overridable :: default_real_num_hidden_layers,       double_precision_num_hidden_layers
+    procedure, private, non_overridable :: default_real_num_inputs,              double_precision_num_inputs
+    procedure, private, non_overridable :: default_real_num_outputs,             double_precision_num_outputs
+    procedure, private, non_overridable :: default_real_nodes_per_layer,         double_precision_nodes_per_layer
+    procedure, private, non_overridable :: default_real_assert_conformable_with, double_precision_assert_conformable_with
+    procedure, private, non_overridable :: default_real_skip,                    double_precision_skip
+    procedure, private, non_overridable :: default_real_activation_name,         double_precision_activation_name
+    procedure, private, non_overridable :: default_real_to_exchange,             double_precision_to_exchange
+  end type
+
+  type workspace_t(k)
+    integer, kind :: k = default_real
+    real(k), allocatable, dimension(:,:) :: a
+    real(k), allocatable, dimension(:,:,:) :: dcdw, vdw, sdw, vdwc, sdwc
+    real(k), allocatable, dimension(:,:) :: z, delta, dcdb, vdb, sdb, vdbc, sdbc
+  contains
+    generic :: fully_allocated          => default_real_allocated
+    generic :: allocate_if_necessary    => default_real_allocate
+    procedure, non_overridable, private :: default_real_allocated
+    procedure, non_overridable, private :: default_real_allocate
   end type
 
   type, extends(inference_engine_t) :: unmapped_engine_t
@@ -115,8 +131,35 @@ impure elemental module function double_precision_unmapped_from_json(file) resul
   end interface
 
 
+  interface workspace_t
+
+    pure module function default_real_workspace(inference_engine) result(workspace)
+      implicit none
+      type(inference_engine_t), intent(in) :: inference_engine
+      type(workspace_t) workspace
+    end function
+
+  end interface
+
   interface
 
+    module subroutine default_real_allocate(self, inference_engine)
+      implicit none
+      class(workspace_t), intent(inout) :: self
+      type(inference_engine_t), intent(in) :: inference_engine
+    end subroutine
+
+    pure module function default_real_allocated(self) result(all_allocated)
+      implicit none
+      class(workspace_t), intent(in) :: self
+      logical all_allocated
+    end function
+
+  end interface
+
+
+  interface ! inference_engine_t type-bound procedures
+
     elemental module function default_real_approximately_equal(lhs, rhs) result(lhs_eq_rhs)
       !! The result is true if lhs and rhs are the same to within a tolerance
       implicit none
@@ -299,6 +342,16 @@ pure module function double_precision_skip(self) result(use_skip_connections)
       logical use_skip_connections
     end function
 
+    pure module subroutine default_real_learn(self, mini_batches_arr, cost, adam, learning_rate, workspace)
+      implicit none
+      class(inference_engine_t), intent(inout) :: self
+      type(mini_batch_t), intent(in) :: mini_batches_arr(:)
+      real, intent(out), allocatable, optional :: cost(:)
+      logical, intent(in) :: adam
+      real, intent(in) :: learning_rate
+      type(workspace_t), intent(inout) :: workspace
+    end subroutine
+
   end interface
 
 end module inference_engine_m_
diff --git a/src/inference_engine/inference_engine_s.F90 b/src/inference_engine/inference_engine_s.F90
index 9050888d3..ce8b5c7ea 100644
--- a/src/inference_engine/inference_engine_s.F90
+++ b/src/inference_engine/inference_engine_s.F90
@@ -14,6 +14,7 @@
   end interface
 
   character(len=*), parameter :: acceptable_engine_tag = "0.13.0" ! git tag capable of reading the current json file format
+  integer, parameter :: input_layer = 0 
 
 contains
 
@@ -60,7 +61,6 @@
   module procedure default_real_infer_unmapped
 
     real, allocatable :: a(:,:)
-    integer, parameter :: input_layer = 0 
     integer l
 
     call assert_consistency(self)
@@ -91,7 +91,6 @@
   module procedure double_precision_infer_unmapped
 
     double precision, allocatable :: a(:,:)
-    integer, parameter :: input_layer = 0 
     integer l
 
     call assert_consistency(self)
@@ -122,7 +121,6 @@
   module procedure default_real_infer
 
     real, allocatable :: a(:,:)
-    integer, parameter :: input_layer = 0
     integer l
 
     call assert_consistency(self)
@@ -169,7 +167,6 @@
   module procedure double_precision_infer
 
     double precision, allocatable :: a(:,:)
-    integer, parameter :: input_layer = 0
     integer l
 
     call assert_consistency(self)
@@ -219,8 +216,6 @@ pure subroutine default_real_consistency(self)
 
     class(inference_engine_t), intent(in) :: self
 
-    integer, parameter :: input_layer=0
-
     associate( &
       all_allocated=>[allocated(self%weights_),allocated(self%biases_),allocated(self%nodes_)]&
     )   
@@ -244,7 +239,6 @@ pure subroutine double_precision_consistency(self)
 
     class(inference_engine_t(double_precision)), intent(in) :: self
 
-    integer, parameter :: input_layer=0
 
     associate( &
       all_allocated=>[allocated(self%weights_),allocated(self%biases_),allocated(self%nodes_)]&
@@ -906,18 +900,18 @@ pure subroutine double_precision_consistency(self)
   end procedure
 
   module procedure default_real_num_hidden_layers
-    integer, parameter :: input_layer = 1, output_layer = 1
+    integer, parameter :: num_non_hidden_layers = 2
     call assert_consistency(self)
     associate(num_layers => size(self%nodes_))
-      hidden_layer_count =  num_layers - (input_layer + output_layer)
+      hidden_layer_count =  num_layers - num_non_hidden_layers
     end associate
   end procedure
 
   module procedure double_precision_num_hidden_layers
-    integer, parameter :: input_layer = 1, output_layer = 1
+    integer, parameter :: num_non_hidden_layers = 2
     call assert_consistency(self)
     associate(num_layers => size(self%nodes_))
-      hidden_layer_count =  num_layers - (input_layer + output_layer)
+      hidden_layer_count =  num_layers - num_non_hidden_layers
     end associate
   end procedure
 
@@ -965,4 +959,193 @@ pure subroutine double_precision_consistency(self)
     end associate
   end procedure
 
+  module procedure default_real_learn
+    integer l, batch, mini_batch_size, pair
+    type(tensor_t), allocatable :: inputs(:), expected_outputs(:)
+
+    call assert_consistency(self)
+    call assert(workspace%fully_allocated(), "inference_engine_s(default_real_learn): workspace%fully_allocated()")
+
+    associate(output_layer => ubound(self%nodes_,1))
+
+      associate( &
+        dcdw => workspace%dcdw, vdw => workspace%vdw, sdw   => workspace%sdw  , vdwc => workspace%vdwc, sdwc => workspace%sdwc &
+       ,dcdb => workspace%dcdb, vdb => workspace%vdb, sdb   => workspace%sdb  , vdbc => workspace%vdbc, sdbc => workspace%sdbc &
+       ,a    => workspace%a   , z   => workspace%z  , delta => workspace%delta &
+      )
+        vdw = 0.; sdw = 1.; vdb = 0.; sdb = 1.
+
+        associate(w => self%weights_, b => self%biases_, n => self%nodes_, num_mini_batches => size(mini_batches_arr))
+
+          if (present(cost)) allocate(cost(num_mini_batches))
+        
+          iterate_across_batches: &
+          do batch = 1, num_mini_batches
+
+            dcdw = 0.; dcdb = 0.
+
+#ifndef _CRAYFTN
+            associate(input_output_pairs => mini_batches_arr(batch)%input_output_pairs())
+#else
+            block
+              type(input_output_pair_t), allocatable :: input_output_pairs(:)
+              input_output_pairs = mini_batches_arr(batch)%input_output_pairs()
+#endif  
+              inputs = input_output_pairs%inputs()
+              expected_outputs = input_output_pairs%expected_outputs()
+              mini_batch_size = size(input_output_pairs)
+#ifndef _CRAYFTN
+            end associate
+#else
+            end block
+#endif  
+            sum_cost: &
+            block
+              real, allocatable :: pair_cost(:)
+              if (present(cost)) allocate(pair_cost(mini_batch_size))
+
+#if F2023_LOCALITY
+              iterate_through_batch: &
+              do concurrent (pair = 1:mini_batch_size) local(a,z,delta) reduce(+: dcdb, dcdw)
+
+#elif F2018_LOCALITY
+
+              reduce_gradients: &
+              block
+                real reduce_dcdb(size(dcdb,1),size(dcdb,2),mini_batch_size)
+                real reduce_dcdw(size(dcdw,1),size(dcdw,2),size(dcdw,3),mini_batch_size)
+                reduce_dcdb = 0.
+                reduce_dcdw = 0.
+
+                iterate_through_batch: &
+                do concurrent (pair = 1:mini_batch_size) local(a,z,delta)
+
+#else
+
+              reduce_gradients: &
+              block
+                real reduce_dcdb(size(dcdb,1),size(dcdb,2),mini_batch_size)
+                real reduce_dcdw(size(dcdw,1),size(dcdw,2),size(dcdw,3),mini_batch_size)
+                reduce_dcdb = 0.
+                reduce_dcdw = 0.
+              
+                iterate_through_batch: &
+                do concurrent (pair = 1:mini_batch_size)
+
+                  iteration: &
+                  block
+
+                    real a(maxval(self%nodes_), input_layer:output_layer) ! Activations
+                    real z(size(b,1),size(b,2)), delta(size(b,1),size(b,2))
+#endif
+
+                    a(1:self%num_inputs(), input_layer) = inputs(pair)%values()
+
+                    feed_forward: &
+                    do l = 1,output_layer
+                      z(1:n(l),l) = matmul(w(1:n(l),1:n(l-1),l), a(1:n(l-1),l-1)) + b(1:n(l),l) ! z_j^l =  sum_k(w_jk^{l} a_k^{l-1}) + b_j^l
+                      a(1:n(l),l) = self%activation_%evaluate(z(1:n(l),l))
+                    end do feed_forward
+
+                    associate(y => expected_outputs(pair)%values())
+                      if (present(cost)) pair_cost(pair) = sum((y(1:n(output_layer))-a(1:n(output_layer),output_layer))**2)
+              
+                      delta(1:n(output_layer),output_layer) = (a(1:n(output_layer),output_layer) - y(1:n(output_layer))) &
+                                                             * self%activation_%differentiate(z(1:n(output_layer),output_layer))
+                    end associate
+                
+                    associate(n_hidden => self%num_hidden_layers())
+                      back_propagate_error: &
+                      do l = n_hidden,1,-1
+                        delta(1:n(l),l) = matmul(transpose(w(1:n(l+1),1:n(l),l+1)), delta(1:n(l+1),l+1)) &
+                                         * self%activation_%differentiate(z(1:n(l),l))
+                      end do back_propagate_error
+                    end associate
+
+                    
+
+                    block
+                      integer j
+                      sum_gradients: &
+                      do l = 1,output_layer
+#if F2023_LOCALITY
+                        dcdb(1:n(l),l) = dcdb(1:n(l),l) + delta(1:n(l),l)
+                        do concurrent(j = 1:n(l)) reduce(+: dcdw)
+                          dcdw(j,1:n(l-1),l) = dcdw(j,1:n(l-1),l) + a(1:n(l-1),l-1)*delta(j,l)
+                        end do
+#else
+                        reduce_dcdb(1:n(l),l,pair) = reduce_dcdb(1:n(l),l,pair) + delta(1:n(l),l)
+                        do j = 1,n(l)
+                          reduce_dcdw(j,1:n(l-1),l,pair) = reduce_dcdw(j,1:n(l-1),l,pair) + a(1:n(l-1),l-1)*delta(j,l)
+                        end do
+#endif
+                      end do sum_gradients
+                    end block
+    
+#if F2023_LOCALITY
+              end do iterate_through_batch
+#elif F2018_LOCALITY
+
+                end do iterate_through_batch
+                dcdb = sum(reduce_dcdb,dim=3)
+                dcdw = sum(reduce_dcdw,dim=4)
+
+              end block reduce_gradients
+#else
+                  end block iteration
+                end do iterate_through_batch
+                dcdb = sum(reduce_dcdb,dim=3)
+                dcdw = sum(reduce_dcdw,dim=4)
+  
+              end block reduce_gradients
+#endif
+
+              if (present(cost)) cost(batch) = sum(pair_cost)/(2*mini_batch_size)
+            end block sum_cost
+          
+            if (adam) then
+              adam: &
+              block
+                ! Adam parameters  
+                real, parameter :: beta(*) = [.9, .999]
+                real, parameter :: obeta(*) = [1.- beta(1), 1.- beta(2)]
+                real, parameter :: epsilon = 1.E-08
+
+                associate(alpha => learning_rate)
+                  adam_adjust_weights_and_biases: &
+                  do concurrent(l = 1:output_layer)
+                    dcdw(1:n(l),1:n(l-1),l) = dcdw(1:n(l),1:n(l-1),l)/(mini_batch_size)
+                    vdw(1:n(l),1:n(l-1),l)  = beta(1)*vdw(1:n(l),1:n(l-1),l) + obeta(1)*dcdw(1:n(l),1:n(l-1),l)
+                    sdw (1:n(l),1:n(l-1),l) = beta(2)*sdw(1:n(l),1:n(l-1),l) + obeta(2)*(dcdw(1:n(l),1:n(l-1),l)**2)
+                    vdwc(1:n(l),1:n(l-1),l) = vdw(1:n(l),1:n(l-1),l)/(1.- beta(1)**num_mini_batches)
+                    sdwc(1:n(l),1:n(l-1),l) = sdw(1:n(l),1:n(l-1),l)/(1.- beta(2)**num_mini_batches)
+                    w(1:n(l),1:n(l-1),l) = w(1:n(l),1:n(l-1),l) &
+                      - alpha*vdwc(1:n(l),1:n(l-1),l)/(sqrt(sdwc(1:n(l),1:n(l-1),l))+epsilon) ! Adjust weights
+
+                    dcdb(1:n(l),l) = dcdb(1:n(l),l)/mini_batch_size
+                    vdb(1:n(l),l) = beta(1)*vdb(1:n(l),l) + obeta(1)*dcdb(1:n(l),l)
+                    sdb(1:n(l),l) = beta(2)*sdb(1:n(l),l) + obeta(2)*(dcdb(1:n(l),l)**2)
+                    vdbc(1:n(l),l) = vdb(1:n(l),l)/(1. - beta(1)**num_mini_batches)
+                    sdbc(1:n(l),l) = sdb(1:n(l),l)/(1. - beta(2)**num_mini_batches)
+                    b(1:n(l),l) = b(1:n(l),l) - alpha*vdbc(1:n(l),l)/(sqrt(sdbc(1:n(l),l))+epsilon) ! Adjust weights
+                  end do adam_adjust_weights_and_biases
+                end associate
+              end block adam
+            else
+              associate(eta => learning_rate)
+                adjust_weights_and_biases: &
+                do concurrent(l = 1:output_layer)
+                  dcdb(1:n(l),l) = dcdb(1:n(l),l)/mini_batch_size
+                  b(1:n(l),l) = b(1:n(l),l) - eta*dcdb(1:n(l),l) ! Adjust biases
+                  dcdw(1:n(l),1:n(l-1),l) = dcdw(1:n(l),1:n(l-1),l)/mini_batch_size
+                  w(1:n(l),1:n(l-1),l) = w(1:n(l),1:n(l-1),l) - eta*dcdw(1:n(l),1:n(l-1),l) ! Adjust weights
+                end do adjust_weights_and_biases
+              end associate
+            end if
+          end do iterate_across_batches
+        end associate
+      end associate
+    end associate
+  end procedure default_real_learn
+
 end submodule inference_engine_s
diff --git a/src/inference_engine/trainable_engine_s.F90 b/src/inference_engine/trainable_engine_s.F90
index 155a084cc..f72096fe6 100644
--- a/src/inference_engine/trainable_engine_s.F90
+++ b/src/inference_engine/trainable_engine_s.F90
@@ -412,6 +412,5 @@ pure function e(j,n) result(unit_vector)
   module procedure default_real_map_from_output_training_range
     unnormalized_tensor = self%output_map_%map_from_training_range(tensor)
   end procedure
-  
 
 end submodule trainable_engine_s
diff --git a/src/inference_engine/trainable_network_m.f90 b/src/inference_engine/trainable_network_m.f90
new file mode 100644
index 000000000..3338f0f69
--- /dev/null
+++ b/src/inference_engine/trainable_network_m.f90
@@ -0,0 +1,42 @@
+module trainable_network_m
+  use inference_engine_m_, only : inference_engine_t, workspace_t
+  use mini_batch_m, only : mini_batch_t
+  use kind_parameters_m, only : default_real
+  implicit none
+
+  private
+  public :: trainable_network_t 
+
+  type, extends(inference_engine_t) ::  trainable_network_t(k)
+    integer, kind :: k = default_real 
+    private
+    type(workspace_t), private :: workspace_
+  contains
+    generic :: train           => default_real_train                   
+    procedure, non_overridable :: default_real_train
+  end type
+
+  interface trainable_network_t 
+
+    pure module function default_real_network(inference_engine) result(trainable_network)
+      implicit none
+      type(inference_engine_t), intent(in) :: inference_engine
+      type(trainable_network_t) trainable_network
+    end function 
+ 
+  end interface
+
+  interface
+
+    pure module subroutine default_real_train(self, mini_batches_arr, cost, adam, learning_rate)
+      implicit none
+      class(trainable_network_t), intent(inout) :: self
+      type(mini_batch_t), intent(in) :: mini_batches_arr(:)
+      real, intent(out), allocatable, optional :: cost(:)
+      logical, intent(in) :: adam
+      real, intent(in) :: learning_rate
+    end subroutine
+
+  end interface
+
+end module trainable_network_m
diff --git a/src/inference_engine/trainable_network_s.f90 b/src/inference_engine/trainable_network_s.f90
new file mode 100644
index 000000000..d6e1d83a6
--- /dev/null
+++ b/src/inference_engine/trainable_network_s.f90
@@ -0,0 +1,15 @@
+submodule(trainable_network_m) trainable_network_s
+  implicit none 
+
+contains
+
+  module procedure default_real_network
+    trainable_network%inference_engine_t = inference_engine
+    trainable_network%workspace_ = workspace_t(inference_engine)
+  end procedure
+
+  module procedure default_real_train
+    call self%inference_engine_t%default_real_learn(mini_batches_arr, cost, adam, learning_rate, self%workspace_)
+  end procedure
+
+end submodule trainable_network_s
diff --git a/src/inference_engine/workspace_s.f90 b/src/inference_engine/workspace_s.f90
new file mode 100644
index 000000000..f3d80139c
--- /dev/null
+++ b/src/inference_engine/workspace_s.f90
@@ -0,0 +1,117 @@
+submodule(inference_engine_m_) workspace_s
+  use assert_m, only : assert
+  implicit none
+
+  integer, parameter :: input_layer = 0
+
+contains
+
+  module procedure default_real_workspace
+
+    allocate(workspace%dcdw, mold=inference_engine%weights_) ! Gradient of cost function with respect to weights
+    allocate(workspace%vdw , mold=inference_engine%weights_)
+    allocate(workspace%sdw , mold=inference_engine%weights_)
+    allocate(workspace%vdwc, mold=inference_engine%weights_)
+    allocate(workspace%sdwc, mold=inference_engine%weights_)
+
+    allocate(workspace%dcdb, mold=inference_engine%biases_ ) ! Gradient of cost function with respect with biases
+    allocate(workspace%vdb , mold=inference_engine%biases_ )
+    allocate(workspace%sdb , mold=inference_engine%biases_ )
+    allocate(workspace%vdbc, mold=inference_engine%biases_ )
+    allocate(workspace%sdbc, mold=inference_engine%biases_ )
+
+    ! TODO: #if ! (F2023_LOCALITY || F2018_LOCALITY)
+    !          then don't allocate a, z, and delta
+
+    allocate(workspace%z    , mold=inference_engine%biases_)
+    allocate(workspace%delta, mold=inference_engine%biases_)
+
+    associate(output_layer => ubound(inference_engine%nodes_,1))
+      allocate(workspace%a(maxval(inference_engine%nodes_), input_layer:output_layer)) ! Activations
+    end associate
+
+    call assert(workspace%fully_allocated(), "workspace_s(defalt_real_workspace): workspace allocated")
+
+  end procedure
+
+  module procedure default_real_allocated
+
+    ! TODO: #if ! (F2023_LOCALITY || F2018_LOCALITY)
+    !          then don't check a, z, and delta allocations
+
+    all_allocated = all( [ & 
+      allocated(self%a), allocated(self%dcdw), allocated(self%vdw), allocated(self%sdw), allocated(self%vdwc), allocated(self%sdwc)&
+     ,allocated(self%z), allocated(self%dcdb), allocated(self%vdb), allocated(self%sdb), allocated(self%vdbc), allocated(self%sdbc)&
+     ,allocated(self%delta) &
+    ])
+  end procedure
+
+
+  module procedure default_real_allocate
+
+    call allocate_if_necessary(self%dcdw,  mold=inference_engine%weights_) ! Gradient of cost function with respect to weights
+    call allocate_if_necessary(self%vdw ,  mold=inference_engine%weights_) 
+    call allocate_if_necessary(self%sdw ,  mold=inference_engine%weights_) 
+    call allocate_if_necessary(self%vdwc,  mold=inference_engine%weights_) 
+    call allocate_if_necessary(self%sdwc,  mold=inference_engine%weights_) 
+
+    call allocate_if_necessary(self%dcdb , mold=inference_engine%biases_) ! Gradient of cost function with respect with biases
+    call allocate_if_necessary(self%vdb  , mold=inference_engine%biases_) 
+    call allocate_if_necessary(self%sdb  , mold=inference_engine%biases_) 
+    call allocate_if_necessary(self%vdbc , mold=inference_engine%biases_) 
+    call allocate_if_necessary(self%sdbc , mold=inference_engine%biases_) 
+
+    ! TODO: #if ! (F2023_LOCALITY || F2018_LOCALITY)
+    !          then don't allocate a, z, and delta
+
+    call allocate_if_necessary(self%z    , mold=inference_engine%biases_) 
+    call allocate_if_necessary(self%delta, mold=inference_engine%biases_) 
+
+    associate(output_layer => ubound(inference_engine%nodes_,1))
+      allocate(self%a(maxval(inference_engine%nodes_), input_layer:output_layer)) ! Activations
+    end associate
+
+  contains
+
+    subroutine allocate_if_necessary(array, mold)
+      real, allocatable, intent(inout) :: array(..)
+      real, intent(in) :: mold(..)
+
+      select rank(array)
+        rank(2)
+          select rank(mold)
+            rank(2)
+              if (.not. allocated(array)) then
+                allocate(array, mold=mold)
+              else
+                if (any(shape(array) /= shape(mold))) then
+                  deallocate(array)
+                  allocate(array, mold=mold)
+                end if
+              end if
+            rank default
+              error stop "workspace_s(allocate_if_necessary): mold-rank mismatch with rank-2 'array'"
+          end select
+        rank(3)
+          select rank(mold)
+            rank(3)
+              if (.not. allocated(array)) then
+                allocate(array, mold=mold)
+              else
+                if (any(shape(array) /= shape(mold))) then
+                  deallocate(array)
+                  allocate(array, mold=mold)
+                end if
+              end if
+            rank default
+              error stop "workspace_s(allocate_if_necessary): mold-rank mismatch with rank-3 'array'"
+          end select
+        rank default
+          error stop "workspace_s(allocate_if_necessary): unsupported 'array' rank"
+      end select
+
+    end subroutine allocate_if_necessary
+
+  end procedure default_real_allocate
+
+end submodule workspace_s
diff --git a/src/inference_engine_m.f90 b/src/inference_engine_m.f90
index 22f505671..732b20ca0 100644
--- a/src/inference_engine_m.f90
+++ b/src/inference_engine_m.f90
@@ -13,6 +13,7 @@ module inference_engine_m
  use network_configuration_m, only : network_configuration_t
  use tensor_m, only : tensor_t
  use tensor_map_m, only : tensor_map_t
+ use trainable_network_m, only : trainable_network_t
  use trainable_engine_m, only : trainable_engine_t
  use training_configuration_m, only : training_configuration_t
  use ubounds_m, only : ubounds_t

From 5ef6625dc3cb747f7e96394fd91bdacb4e2a39ae Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 13 Oct 2024 20:43:20 -0700
Subject: [PATCH 2/6] feat(sat-mix-rat):adjust to use trainale_network_t

---
 example/learn-saturated-mixing-ratio.F90 | 45 +++++++++++-------------
 1 file changed, 21 insertions(+), 24 deletions(-)

diff --git a/example/learn-saturated-mixing-ratio.F90 b/example/learn-saturated-mixing-ratio.F90
index ba8bd599c..641f33f3e 100644
--- a/example/learn-saturated-mixing-ratio.F90
+++ b/example/learn-saturated-mixing-ratio.F90
@@ -3,7 +3,7 @@
 program train_saturated_mixture_ratio
   !! This program trains a neural network to learn the saturated mixing ratio function of ICAR.
   use inference_engine_m, only : &
-    inference_engine_t, trainable_engine_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
+    inference_engine_t, trainable_network_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
   use julienne_m, only : string_t, file_t, command_line_t, bin_t, csv
   use assert_m, only : assert, intrinsic_array_t
   use saturated_mixing_ratio_m, only : y, T, p
@@ -30,29 +30,29 @@ program train_saturated_mixture_ratio
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(input_output_pair_t), allocatable :: input_output_pairs(:)
     type(tensor_t), allocatable :: inputs(:), desired_outputs(:)
-    type(trainable_engine_t)  trainable_engine
+    type(trainable_network_t)  trainable_network
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:), random_numbers(:)
     integer io_status, network_unit, plot_unit
-    integer, parameter :: io_success=0, diagnostics_print_interval = 1000, network_save_interval = 10000
+    integer, parameter :: io_success=0, diagnostics_print_interval = 100, network_save_interval = 1000
     integer, parameter :: nodes_per_layer(*) = [2, 4, 72, 2, 1]
-    real, parameter :: cost_tolerance = 1.E-08
+    real, parameter :: cost_tolerance = 1.E-06
 
     call random_init(image_distinct=.true., repeatable=.true.)
     open(newunit=network_unit, file=network_file%string(), form='formatted', status='old', iostat=io_status, action='read')
 
     if (io_status == io_success) then
       print *,"Reading network from file " // network_file%string()
-      trainable_engine = trainable_engine_t(inference_engine_t(file_t(network_file)))
+      trainable_network = trainable_network_t(inference_engine_t(file_t(network_file)))
       close(network_unit)
     else
       close(network_unit)
       print *,"Initializing a new network"
-      trainable_engine = perturbed_identity_network(perturbation_magnitude=0.05, n = nodes_per_layer)
+      trainable_network = perturbed_identity_network(perturbation_magnitude=0.05, n = nodes_per_layer)
     end if
-    call output(trainable_engine%to_inference_engine(), string_t("initial-network.json"))
+    call output(trainable_network, string_t("initial-network.json"))
 
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
+    associate(num_inputs => trainable_network%num_inputs(), num_outputs => trainable_network%num_outputs())
 
       block
         integer i, j
@@ -84,7 +84,7 @@ program train_saturated_mixture_ratio
           call random_number(random_numbers)
           call shuffle(input_output_pairs)
           mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))]
-          call trainable_engine%train(mini_batches, cost, adam=.true., learning_rate=1.5)
+          call trainable_network%train(mini_batches, cost, adam=.true., learning_rate=1.5)
           call system_clock(counter_end, clock_rate)
 
           associate( &
@@ -95,7 +95,7 @@ program train_saturated_mixture_ratio
             write_and_exit_if_converged: &
             if (cost_avg < cost_tolerance) then
               call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer)
-              call output(trainable_engine%to_inference_engine(), network_file)
+              call output(trainable_network, network_file)
               exit
             end if write_and_exit_if_converged
 
@@ -104,13 +104,13 @@ program train_saturated_mixture_ratio
             write_and_exit_if_stop_file_exists: &
             if (io_status==0) then
               call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer)
-              call output(trainable_engine%to_inference_engine(), network_file)
+              call output(trainable_network, network_file)
               exit
             end if write_and_exit_if_stop_file_exists
 
             if (mod(e,diagnostics_print_interval)==0 .or. loop_ending) &
               call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer)
-            if (mod(e,network_save_interval)==0 .or. loop_ending) call output(trainable_engine%to_inference_engine(), network_file)
+            if (mod(e,network_save_interval)==0 .or. loop_ending) call output(trainable_network, network_file)
           end associate
         end do
 
@@ -121,9 +121,9 @@ program train_saturated_mixture_ratio
           integer p
 #if defined _CRAYFTN || __GFORTRAN__
           type(tensor_t), allocatable :: network_outputs(:)
-          network_outputs = trainable_engine%infer(inputs)
+          network_outputs = trainable_network%infer(inputs)
 #else
-          associate(network_outputs => trainable_engine%infer(inputs))
+          associate(network_outputs => trainable_network%infer(inputs))
 #endif
             print *,"Inputs (normalized)          | Outputs      | Desired outputs"
             do p = 1, num_pairs
@@ -139,7 +139,7 @@ program train_saturated_mixture_ratio
 
    end associate
 
-   call output(trainable_engine%to_inference_engine(), network_file)
+   call output(trainable_network, network_file)
 
   end block
 
@@ -156,7 +156,7 @@ subroutine print_diagnostics(plot_file_unit, epoch, cost, clock, nodes)
   end subroutine
 
   subroutine output(inference_engine, file_name)
-    type(inference_engine_t), intent(in) :: inference_engine
+    class(inference_engine_t), intent(in) :: inference_engine
     type(string_t), intent(in) :: file_name
     type(file_t) json_file
     json_file = inference_engine%to_json()
@@ -170,8 +170,8 @@ pure function e(j,n) result(unit_vector)
     unit_vector = real([(merge(1,0,j==k),k=1,n)])
   end function
 
-  function perturbed_identity_network(perturbation_magnitude, n) result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function perturbed_identity_network(perturbation_magnitude, n) result(trainable_network)
+    type(trainable_network_t) trainable_network
     real, intent(in) :: perturbation_magnitude
     integer, intent(in)  :: n(:)
     integer k, l
@@ -187,12 +187,9 @@ function perturbed_identity_network(perturbation_magnitude, n) result(trainable_
       call random_number(b_harvest)
 
       associate(w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, b => perturbation_magnitude*(b_harvest-0.5)/0.5)
-
-        trainable_engine = trainable_engine_t( &
-          nodes = n, weights = w, biases = b, metadata = &
-            [string_t("Saturated Mixing Ratio"), string_t("Damian Rouson"), string_t("2023-09-23"), string_t("relu"), &
-             string_t("false")] &
-        )
+        trainable_network = trainable_network_t( inference_engine_t(nodes=n, weights=w, biases=b, &
+          metadata=[string_t("Saturated Mixing Ratio"),string_t("Rouson"),string_t("20241013"),string_t("relu"),string_t("false")] &
+        ))
       end associate
     end associate
   end function

From d45c57eebf4ff5f1b3ff8716b0e528da04af11a8 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 13 Oct 2024 20:56:46 -0700
Subject: [PATCH 3/6] fix(activation_t):add :: in default initialization

---
 src/inference_engine/activation_m.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/inference_engine/activation_m.f90 b/src/inference_engine/activation_m.f90
index 3566fbff8..966f6beea 100644
--- a/src/inference_engine/activation_m.f90
+++ b/src/inference_engine/activation_m.f90
@@ -15,7 +15,7 @@ module activation_m
 
   type activation_t
     private
-    integer(c_int) selection_ = sigmoid
+    integer(c_int) :: selection_ = sigmoid
   contains
     procedure, non_overridable :: function_name
     generic :: operator(==) => equals

From 4c76b7a85a94f3dc031ad1f80f2236963dbc908c Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 13 Oct 2024 21:50:51 -0700
Subject: [PATCH 4/6] fix(unmapped_network): rm ambiguous specific proc

also adjust demo/app/infer-aerosol.f90 to reflect the
renaming of unmapped_engine_t to unmapped_network_t
---
 demo/app/infer-aerosol.f90                   |   6 +-
 src/inference_engine/inference_engine_m_.f90 | 158 +++++++++++--------
 src/inference_engine/inference_engine_s.F90  | 142 ++++-------------
 src/inference_engine/unmapped_network_s.f90  |  88 +++++++++++
 src/inference_engine_m.f90                   |   2 +-
 5 files changed, 216 insertions(+), 180 deletions(-)
 create mode 100644 src/inference_engine/unmapped_network_s.f90

diff --git a/demo/app/infer-aerosol.f90 b/demo/app/infer-aerosol.f90
index 6810af0cc..d2211f2e9 100644
--- a/demo/app/infer-aerosol.f90
+++ b/demo/app/infer-aerosol.f90
@@ -7,7 +7,7 @@ program infer_aerosol
   use iso_fortran_env, only : int64, real64
 
   ! External dependencies:
-  use inference_engine_m, only : unmapped_engine_t, tensor_t, double_precision, double_precision_file_t
+  use inference_engine_m, only : unmapped_network_t, tensor_t, double_precision, double_precision_file_t
   use julienne_m, only : string_t, command_line_t
   use omp_lib
 
@@ -49,7 +49,7 @@ subroutine read_stats_and_perform_inference(path)
     double precision cube_root
     double precision, allocatable, dimension(:,:) :: aerosol_data, input_components, output_components
     type(tensor_statistics_t) input_stats, output_stats
-    type(unmapped_engine_t(double_precision)) inference_engine
+    type(unmapped_network_t(double_precision)) inference_engine
     integer i, j
 
     input_stats  = read_tensor_statistics(path // "meanxp.txt", path // "stdxp.txt", num_inputs)  !for pre-processing normalization
@@ -76,7 +76,7 @@ subroutine read_stats_and_perform_inference(path)
     !$omp end parallel do   
       
     print *, "Reading the neural network from " // network_file_name
-    inference_engine = unmapped_engine_t(double_precision_file_t(path // network_file_name))
+    inference_engine = unmapped_network_t(double_precision_file_t(path // network_file_name))
 
     time_inference: &
     block
diff --git a/src/inference_engine/inference_engine_m_.f90 b/src/inference_engine/inference_engine_m_.f90
index a196fd6ee..40121b4cf 100644
--- a/src/inference_engine/inference_engine_m_.f90
+++ b/src/inference_engine/inference_engine_m_.f90
@@ -14,7 +14,7 @@ module inference_engine_m_
 
   private
   public :: inference_engine_t
-  public :: unmapped_engine_t
+  public :: unmapped_network_t
   public :: exchange_t
   public :: workspace_t
 
@@ -41,6 +41,8 @@ module inference_engine_m_
     generic :: activation_function_name => default_real_activation_name,         double_precision_activation_name
     generic :: to_exchange              => default_real_to_exchange,             double_precision_to_exchange
     generic :: learn                    => default_real_learn
+    generic :: assert_consistency       => default_real_consistency,             double_precision_consistency
+    procedure, private, non_overridable :: default_real_consistency,             double_precision_consistency
     procedure, private, non_overridable :: default_real_approximately_equal,     double_precision_approximately_equal
     procedure, private, non_overridable :: default_real_infer,                   double_precision_infer
     procedure, private, non_overridable :: default_real_learn
@@ -69,11 +71,31 @@ module inference_engine_m_
     procedure, non_overridable, private :: default_real_allocate
   end type
 
-  type, extends(inference_engine_t) :: unmapped_engine_t
-  contains
-    generic :: infer =>   default_real_infer_unmapped, double_precision_infer_unmapped
-    procedure, private :: default_real_infer_unmapped, double_precision_infer_unmapped
-  end type
+  interface workspace_t
+
+    pure module function default_real_workspace(inference_engine) result(workspace)
+      implicit none
+      type(inference_engine_t), intent(in) :: inference_engine
+      type(workspace_t) workspace
+    end function
+
+  end interface
+
+  interface
+
+    module subroutine default_real_allocate(self, inference_engine)
+      implicit none
+      class(workspace_t), intent(inout) :: self
+      type(inference_engine_t), intent(in) :: inference_engine
+    end subroutine
+
+    pure module function default_real_allocated(self) result(all_allocated)
+      implicit none
+      class(workspace_t), intent(in) :: self
+      logical all_allocated
+    end function
+
+  end interface
 
   type exchange_t(k)
     integer, kind :: k = default_real
@@ -84,6 +106,22 @@ module inference_engine_m_
     type(activation_t) activation_
   end type
 
+  interface
+
+    module function default_real_to_exchange(self) result(exchange)
+      implicit none
+      class(inference_engine_t), intent(in) :: self
+      type(exchange_t) exchange
+    end function
+
+    module function double_precision_to_exchange(self) result(exchange)
+      implicit none
+      class(inference_engine_t(double_precision)), intent(in) :: self
+      type(exchange_t(double_precision)) exchange
+    end function
+
+  end interface
+
   interface inference_engine_t
 
     module function default_real_construct_from_components(metadata, weights, biases, nodes, input_map, output_map) &
@@ -120,42 +158,7 @@ impure elemental module function double_precision_from_json(file) result(inferen
 
   end interface
 
-  interface unmapped_engine_t
-
-    impure elemental module function double_precision_unmapped_from_json(file) result(unmapped_engine)
-      implicit none
-      type(double_precision_file_t), intent(in) :: file
-      type(unmapped_engine_t(double_precision)) unmapped_engine
-    end function
-
-  end interface
-
-
-  interface workspace_t
-
-    pure module function default_real_workspace(inference_engine) result(workspace)
-      implicit none
-      type(inference_engine_t), intent(in) :: inference_engine
-      type(workspace_t) workspace
-    end function
-
-  end interface
-
-  interface
-
-    module subroutine default_real_allocate(self, inference_engine)
-      implicit none
-      class(workspace_t), intent(inout) :: self
-      type(inference_engine_t), intent(in) :: inference_engine
-    end subroutine
-
-    pure module function default_real_allocated(self) result(all_allocated)
-      implicit none
-      class(workspace_t), intent(in) :: self
-      logical all_allocated
-    end function
 
-  end interface
 
 
   interface ! inference_engine_t type-bound procedures
@@ -206,18 +209,6 @@ elemental module function double_precision_map_from_output_range(self, normalize
       type(tensor_t(double_precision)) tensor
     end function
 
-    module function default_real_to_exchange(self) result(exchange)
-      implicit none
-      class(inference_engine_t), intent(in) :: self
-      type(exchange_t) exchange
-    end function
-
-    module function double_precision_to_exchange(self) result(exchange)
-      implicit none
-      class(inference_engine_t(double_precision)), intent(in) :: self
-      type(exchange_t(double_precision)) exchange
-    end function
-
     impure elemental module function default_real_to_json(self) result(json_file)
       implicit none
       class(inference_engine_t), intent(in) :: self
@@ -249,20 +240,6 @@ elemental module function default_real_infer(self, inputs) result(outputs)
       type(tensor_t) outputs
     end function
 
-    elemental module function default_real_infer_unmapped(self, inputs) result(outputs)
-      implicit none
-      class(unmapped_engine_t), intent(in) :: self
-      type(tensor_t), intent(in) :: inputs
-      type(tensor_t) outputs
-    end function
-
-    elemental module function double_precision_infer_unmapped(self, inputs) result(outputs)
-      implicit none
-      class(unmapped_engine_t(double_precision)), intent(in) :: self
-      type(tensor_t(double_precision)), intent(in) :: inputs
-      type(tensor_t(double_precision)) outputs
-    end function
-
     elemental module function double_precision_infer(self, inputs) result(outputs)
       implicit none
       class(inference_engine_t(double_precision)), intent(in) :: self
@@ -352,6 +329,53 @@ pure module subroutine default_real_learn(self, mini_batches_arr, cost, adam, le
       type(workspace_t), intent(inout) :: workspace
     end subroutine
 
+    pure module subroutine default_real_consistency(self)
+      implicit none
+      class(inference_engine_t), intent(in) :: self
+    end subroutine
+
+    pure module subroutine double_precision_consistency(self)
+      implicit none
+      class(inference_engine_t(double_precision)), intent(in) :: self
+    end subroutine
+
+  end interface
+
+  type unmapped_network_t(k)
+    integer, kind :: k = default_real
+    private
+    type(inference_engine_t(k)) inference_engine_
+  contains
+    generic :: infer                    => default_real_infer_unmapped, double_precision_infer_unmapped
+    procedure, private, non_overridable :: default_real_infer_unmapped, double_precision_infer_unmapped
+  end type
+
+  interface unmapped_network_t
+
+    impure elemental module function double_precision_unmapped_from_json(file) result(unmapped_network)
+      implicit none
+      type(double_precision_file_t), intent(in) :: file
+      type(unmapped_network_t(double_precision)) unmapped_network
+    end function
+
+  end interface
+
+  interface
+
+    elemental module function default_real_infer_unmapped(self, inputs) result(outputs)
+      implicit none
+      class(unmapped_network_t), intent(in) :: self
+      type(tensor_t), intent(in) :: inputs
+      type(tensor_t) outputs
+    end function
+
+    elemental module function double_precision_infer_unmapped(self, inputs) result(outputs)
+      implicit none
+      class(unmapped_network_t(double_precision)), intent(in) :: self
+      type(tensor_t(double_precision)), intent(in) :: inputs
+      type(tensor_t(double_precision)) outputs
+    end function
+
   end interface
 
 end module inference_engine_m_
diff --git a/src/inference_engine/inference_engine_s.F90 b/src/inference_engine/inference_engine_s.F90
index ce8b5c7ea..ac31642f6 100644
--- a/src/inference_engine/inference_engine_s.F90
+++ b/src/inference_engine/inference_engine_s.F90
@@ -8,11 +8,6 @@
   use neuron_m, only : neuron_t
   implicit none
 
-  interface assert_consistency
-    procedure default_real_consistency
-    procedure double_precision_consistency
-  end interface
-
   character(len=*), parameter :: acceptable_engine_tag = "0.13.0" ! git tag capable of reading the current json file format
   integer, parameter :: input_layer = 0 
 
@@ -58,72 +53,12 @@
     exchange%activation_ = self%activation_
   end procedure
 
-  module procedure default_real_infer_unmapped
-
-    real, allocatable :: a(:,:)
-    integer l
-
-    call assert_consistency(self)
-
-    associate(w => self%weights_, b => self%biases_, n => self%nodes_, output_layer => ubound(self%nodes_,1))
-
-      allocate(a(maxval(n), input_layer:output_layer))
-
-      a(1:n(input_layer),input_layer) = inputs%values()
-
-      feed_forward: &
-      do l = input_layer+1, output_layer
-      associate(z => matmul(w(1:n(l),1:n(l-1),l), a(1:n(l-1),l-1)) + b(1:n(l),l))
-          if (l .lt. output_layer) then
-             a(1:n(l),l) = self%activation_%evaluate(z)
-          else
-             a(1:n(l),l) = z(1:n(l))
-          end if
-        end associate
-      end do feed_forward
-
-      outputs = tensor_t(a(1:n(output_layer), output_layer))
-
-    end associate
-
-  end procedure
-
-  module procedure double_precision_infer_unmapped
-
-    double precision, allocatable :: a(:,:)
-    integer l
-
-    call assert_consistency(self)
-
-    associate(w => self%weights_, b => self%biases_, n => self%nodes_, output_layer => ubound(self%nodes_,1))
-
-      allocate(a(maxval(n), input_layer:output_layer))
-
-      a(1:n(input_layer),input_layer) = inputs%values()
-
-      feed_forward: &
-      do l = input_layer+1, output_layer
-      associate(z => matmul(w(1:n(l),1:n(l-1),l), a(1:n(l-1),l-1)) + b(1:n(l),l))
-          if (l .lt. output_layer) then
-             a(1:n(l),l) = self%activation_%evaluate(z)
-          else
-             a(1:n(l),l) = z(1:n(l))
-          end if
-        end associate
-      end do feed_forward
-
-      outputs = tensor_t(a(1:n(output_layer), output_layer))
-
-    end associate
-
-  end procedure
-
   module procedure default_real_infer
 
     real, allocatable :: a(:,:)
     integer l
 
-    call assert_consistency(self)
+    call self%assert_consistency()
 
     associate(w => self%weights_, b => self%biases_, n => self%nodes_, output_layer => ubound(self%nodes_,1))
 
@@ -169,7 +104,7 @@
     double precision, allocatable :: a(:,:)
     integer l
 
-    call assert_consistency(self)
+    call self%assert_consistency()
 
     associate(w => self%weights_, b => self%biases_, n => self%nodes_, output_layer => ubound(self%nodes_,1))
 
@@ -212,52 +147,47 @@
 
   end procedure
 
-  pure subroutine default_real_consistency(self)
-
-    class(inference_engine_t), intent(in) :: self
+  module procedure default_real_consistency
 
     associate( &
       all_allocated=>[allocated(self%weights_),allocated(self%biases_),allocated(self%nodes_)]&
     )   
-      call assert(all(all_allocated),"inference_engine_s(inference_engine_consistency): fully_allocated", &
+      call assert(all(all_allocated),"inference_engine_s(default_real_consistency): fully_allocated", &
         intrinsic_array_t(all_allocated))
     end associate
 
     associate(max_width=>maxval(self%nodes_), component_dims=>[size(self%biases_,1), size(self%weights_,1), size(self%weights_,2)])
-      call assert(all(component_dims == max_width), "inference_engine_s(inference_engine_consistency): conformable arrays", &
+      call assert(all(component_dims == max_width), "inference_engine_s(default_real_consistency): conformable arrays", &
         intrinsic_array_t([max_width,component_dims]))
     end associate
 
     associate(input_subscript => lbound(self%nodes_,1))
-      call assert(input_subscript == input_layer, "inference_engine_s(inference_engine_consistency): n base subsscript", &
+      call assert(input_subscript == input_layer, "inference_engine_s(default_real_consistency): n base subsscript", &
         input_subscript)
     end associate
 
-  end subroutine
-
-  pure subroutine double_precision_consistency(self)
-
-    class(inference_engine_t(double_precision)), intent(in) :: self
+  end procedure
 
+  module procedure double_precision_consistency
 
     associate( &
       all_allocated=>[allocated(self%weights_),allocated(self%biases_),allocated(self%nodes_)]&
     )   
-      call assert(all(all_allocated),"inference_engine_s(inference_engine_consistency): fully_allocated", &
+      call assert(all(all_allocated),"inference_engine_s(default_real_consistency): fully_allocated", &
         intrinsic_array_t(all_allocated))
     end associate
 
     associate(max_width=>maxval(self%nodes_), component_dims=>[size(self%biases_,1), size(self%weights_,1), size(self%weights_,2)])
-      call assert(all(component_dims == max_width), "inference_engine_s(inference_engine_consistency): conformable arrays", &
+      call assert(all(component_dims == max_width), "inference_engine_s(default_real_consistency): conformable arrays", &
         intrinsic_array_t([max_width,component_dims]))
     end associate
 
     associate(input_subscript => lbound(self%nodes_,1))
-      call assert(input_subscript == input_layer, "inference_engine_s(inference_engine_consistency): n base subsscript", &
+      call assert(input_subscript == input_layer, "inference_engine_s(default_real_consistency): n base subsscript", &
         input_subscript)
     end associate
 
-  end subroutine
+  end procedure
 
   module procedure default_real_construct_from_components
 
@@ -292,7 +222,7 @@ pure subroutine double_precision_consistency(self)
 
     inference_engine%activation_ = activation_t(metadata(4)%string())
 
-    call assert_consistency(inference_engine)
+    call inference_engine%assert_consistency()
 
   end procedure default_real_construct_from_components
 
@@ -331,7 +261,7 @@ pure subroutine double_precision_consistency(self)
       inference_engine%activation_ = activation_t(function_name%string())
     end associate
 
-    call assert_consistency(inference_engine)
+    call inference_engine%assert_consistency()
 
   end procedure double_precision_construct_from_components
 
@@ -346,7 +276,7 @@ pure subroutine double_precision_consistency(self)
     proto_neuron = neuron_t([0.],1.)
 #endif
 
-    call assert_consistency(self)
+    call self%assert_consistency()
 
     associate( &
        num_hidden_layers => self%num_hidden_layers() &
@@ -462,7 +392,7 @@ pure subroutine double_precision_consistency(self)
     proto_neuron = neuron_t([0D0],1D0)
 #endif
 
-    call assert_consistency(self)
+    call self%assert_consistency()
 
     associate( &
        num_hidden_layers => self%num_hidden_layers() &
@@ -665,14 +595,10 @@ pure subroutine double_precision_consistency(self)
       end associate
     end associate read_metadata
 
-    call assert_consistency(inference_engine)
+    call inference_engine%assert_consistency()
 
   end procedure default_real_from_json
 
-  module procedure double_precision_unmapped_from_json
-    unmapped_engine%inference_engine_t = double_precision_from_json(file)
-  end procedure
-
   module procedure double_precision_from_json
 
     character(len=:), allocatable :: justified_line
@@ -769,14 +695,13 @@ pure subroutine double_precision_consistency(self)
       end associate
     end associate read_metadata
 
-    call assert_consistency(inference_engine)
+    call inference_engine%assert_consistency()
 
   end procedure double_precision_from_json
 
   module procedure default_real_assert_conformable_with
 
-    call assert_consistency(self)
-    call assert_consistency(inference_engine)
+    call self%assert_consistency()
 
     associate(equal_shapes => [ &
       shape(self%weights_) == shape(inference_engine%weights_), &
@@ -792,8 +717,7 @@ pure subroutine double_precision_consistency(self)
 
   module procedure double_precision_assert_conformable_with
 
-    call assert_consistency(self)
-    call assert_consistency(inference_engine)
+    call self%assert_consistency()
 
     associate(equal_shapes => [ &
       shape(self%weights_) == shape(inference_engine%weights_), &
@@ -813,8 +737,8 @@ pure subroutine double_precision_consistency(self)
 
     nodes_eq = all(lhs%nodes_ == rhs%nodes_)
 
-    call assert_consistency(lhs)
-    call assert_consistency(rhs)
+    call lhs%assert_consistency()
+    call rhs%assert_consistency()
     call lhs%assert_conformable_with(rhs)
 
     block
@@ -854,8 +778,8 @@ pure subroutine double_precision_consistency(self)
 
     nodes_eq = all(lhs%nodes_ == rhs%nodes_)
 
-    call assert_consistency(lhs)
-    call assert_consistency(rhs)
+    call lhs%assert_consistency()
+    call rhs%assert_consistency()
     call lhs%assert_conformable_with(rhs)
 
     block
@@ -890,18 +814,18 @@ pure subroutine double_precision_consistency(self)
   end procedure
 
   module procedure default_real_num_outputs
-    call assert_consistency(self)
+    call self%assert_consistency()
     output_count = self%nodes_(ubound(self%nodes_,1))
   end procedure
 
   module procedure double_precision_num_outputs
-    call assert_consistency(self)
+    call self%assert_consistency()
     output_count = self%nodes_(ubound(self%nodes_,1))
   end procedure
 
   module procedure default_real_num_hidden_layers
     integer, parameter :: num_non_hidden_layers = 2
-    call assert_consistency(self)
+    call self%assert_consistency()
     associate(num_layers => size(self%nodes_))
       hidden_layer_count =  num_layers - num_non_hidden_layers
     end associate
@@ -909,29 +833,29 @@ pure subroutine double_precision_consistency(self)
 
   module procedure double_precision_num_hidden_layers
     integer, parameter :: num_non_hidden_layers = 2
-    call assert_consistency(self)
+    call self%assert_consistency()
     associate(num_layers => size(self%nodes_))
       hidden_layer_count =  num_layers - num_non_hidden_layers
     end associate
   end procedure
 
   module procedure default_real_num_inputs
-    call assert_consistency(self)
+    call self%assert_consistency()
     input_count = self%nodes_(lbound(self%nodes_,1))
   end procedure
 
   module procedure double_precision_num_inputs
-    call assert_consistency(self)
+    call self%assert_consistency()
     input_count = self%nodes_(lbound(self%nodes_,1))
   end procedure
 
   module procedure default_real_nodes_per_layer
-    call assert_consistency(self)
+    call self%assert_consistency()
     node_count = self%nodes_
   end procedure
 
   module procedure double_precision_nodes_per_layer
-    call assert_consistency(self)
+    call self%assert_consistency()
     node_count = self%nodes_
   end procedure
 
@@ -963,7 +887,7 @@ pure subroutine double_precision_consistency(self)
     integer l, batch, mini_batch_size, pair
     type(tensor_t), allocatable :: inputs(:), expected_outputs(:)
 
-    call assert_consistency(self)
+    call self%assert_consistency()
     call assert(workspace%fully_allocated(), "inference_engine_s(default_real_learn): workspace%fully_allocated()")
 
     associate(output_layer => ubound(self%nodes_,1))
diff --git a/src/inference_engine/unmapped_network_s.f90 b/src/inference_engine/unmapped_network_s.f90
new file mode 100644
index 000000000..9ced289e5
--- /dev/null
+++ b/src/inference_engine/unmapped_network_s.f90
@@ -0,0 +1,88 @@
+! Copyright (c), The Regents of the University of California
+! Terms of use are as specified in LICENSE.txt
+submodule(inference_engine_m_) unmapped_network_s
+  implicit none
+
+  integer, parameter :: input_layer = 0 
+
+contains
+
+  module procedure double_precision_unmapped_from_json
+    unmapped_network%inference_engine_ = double_precision_from_json(file)
+  end procedure
+
+  module procedure default_real_infer_unmapped
+
+    real, allocatable :: a(:,:)
+    integer l
+
+    associate(inference_engine => self%inference_engine_)
+
+      call inference_engine%assert_consistency()
+
+      associate( &
+        w => inference_engine%weights_ &
+       ,b => inference_engine%biases_ &
+       ,n => inference_engine%nodes_ &
+       ,output_layer => ubound(inference_engine%nodes_,1) &
+      )
+        allocate(a(maxval(n), input_layer:output_layer))
+
+        a(1:n(input_layer),input_layer) = inputs%values()
+
+        feed_forward: &
+        do l = input_layer+1, output_layer
+          associate(z => matmul(w(1:n(l),1:n(l-1),l), a(1:n(l-1),l-1)) + b(1:n(l),l))
+            if (l .lt. output_layer) then
+               a(1:n(l),l) = inference_engine%activation_%evaluate(z)
+            else
+               a(1:n(l),l) = z(1:n(l))
+            end if
+          end associate
+        end do feed_forward
+
+        outputs = tensor_t(a(1:n(output_layer), output_layer))
+
+      end associate
+    end associate
+
+  end procedure
+
+  module procedure double_precision_infer_unmapped
+
+    double precision, allocatable :: a(:,:)
+    integer l
+
+    associate(inference_engine => self%inference_engine_)
+
+      call inference_engine%assert_consistency()
+
+      associate( &
+        w => inference_engine%weights_ &
+       ,b => inference_engine%biases_ &
+       ,n => inference_engine%nodes_ &
+       ,output_layer => ubound(inference_engine%nodes_,1) &
+      )
+
+        allocate(a(maxval(n), input_layer:output_layer))
+
+        a(1:n(input_layer),input_layer) = inputs%values()
+
+        feed_forward: &
+        do l = input_layer+1, output_layer
+          associate(z => matmul(w(1:n(l),1:n(l-1),l), a(1:n(l-1),l-1)) + b(1:n(l),l))
+            if (l .lt. output_layer) then
+               a(1:n(l),l) = inference_engine%activation_%evaluate(z)
+            else
+               a(1:n(l),l) = z(1:n(l))
+            end if
+          end associate
+        end do feed_forward
+
+        outputs = tensor_t(a(1:n(output_layer), output_layer))
+      end associate
+    end associate
+
+  end procedure
+
+end submodule unmapped_network_s
diff --git a/src/inference_engine_m.f90 b/src/inference_engine_m.f90
index 732b20ca0..5b6820d1e 100644
--- a/src/inference_engine_m.f90
+++ b/src/inference_engine_m.f90
@@ -6,7 +6,7 @@ module inference_engine_m
  use double_precision_string_m, only : double_precision_string_t
  use hyperparameters_m, only : hyperparameters_t
  use input_output_pair_m, only : input_output_pair_t, shuffle, write_to_stdout
- use inference_engine_m_, only : inference_engine_t, unmapped_engine_t
+ use inference_engine_m_, only : inference_engine_t, unmapped_network_t
  use kind_parameters_m, only : default_real, double_precision
  use metadata_m, only : metadata_t
  use mini_batch_m, only : mini_batch_t

From e4f8ba9c8784c373d98529b43c2966491c085582 Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Sun, 13 Oct 2024 23:35:30 -0700
Subject: [PATCH 5/6] chore(trainable_network):replaces trainable_engine

---
 demo/app/train-cloud-microphysics.F90         |  20 +-
 example/learn-addition.F90                    |  29 +-
 example/learn-exponentiation.F90              |  29 +-
 example/learn-microphysics-procedures.F90     | 236 ----------
 example/learn-multiplication.F90              |  29 +-
 example/learn-power-series.F90                |  29 +-
 example/learn-saturated-mixing-ratio.F90      |   4 +-
 example/supporting-modules/mp_thompson.f90    | 102 -----
 .../supporting-modules/thompson_tensors_m.f90 |  33 --
 example/train-and-write.F90                   |  39 +-
 src/inference_engine/tmp                      |   0
 src/inference_engine/trainable_engine_m.F90   | 179 --------
 src/inference_engine/trainable_engine_s.F90   | 416 ------------------
 src/inference_engine/trainable_network_m.f90  |  29 +-
 src/inference_engine/trainable_network_s.f90  |  48 +-
 src/inference_engine_m.f90                    |   1 -
 test/main.F90                                 |   6 +-
 ...est_m.F90 => trainable_network_test_m.F90} | 106 ++---
 18 files changed, 216 insertions(+), 1119 deletions(-)
 delete mode 100644 example/learn-microphysics-procedures.F90
 delete mode 100644 example/supporting-modules/mp_thompson.f90
 delete mode 100644 example/supporting-modules/thompson_tensors_m.f90
 create mode 100644 src/inference_engine/tmp
 delete mode 100644 src/inference_engine/trainable_engine_m.F90
 delete mode 100644 src/inference_engine/trainable_engine_s.F90
 rename test/{trainable_engine_test_m.F90 => trainable_network_test_m.F90} (86%)

diff --git a/demo/app/train-cloud-microphysics.F90 b/demo/app/train-cloud-microphysics.F90
index fcb952056..e96038f5d 100644
--- a/demo/app/train-cloud-microphysics.F90
+++ b/demo/app/train-cloud-microphysics.F90
@@ -13,7 +13,7 @@ program train_on_flat_distribution
   use julienne_m, only : string_t, file_t, command_line_t, bin_t
   use assert_m, only : assert, intrinsic_array_t
   use inference_engine_m, only : &
-    inference_engine_t, mini_batch_t, input_output_pair_t, tensor_t, trainable_engine_t, tensor_map_t, &
+    inference_engine_t, mini_batch_t, input_output_pair_t, tensor_t, trainable_network_t, tensor_map_t, &
     training_configuration_t, shuffle
 
   !! Internal dependencies:
@@ -271,7 +271,7 @@ subroutine read_train_write(training_configuration, args, plot_file)
 
       train_network: &
       block
-        type(trainable_engine_t) trainable_engine
+        type(trainable_network_t) trainable_network
         type(mini_batch_t), allocatable :: mini_batches(:)
         type(bin_t), allocatable :: bins(:)
         type(input_output_pair_t), allocatable :: input_output_pairs(:)
@@ -319,7 +319,7 @@ subroutine read_train_write(training_configuration, args, plot_file)
             read_or_initialize_engine: &
             if (io_status==0) then
               print *,"Reading network from file " // network_file
-              trainable_engine = trainable_engine_t(inference_engine_t(file_t(string_t(network_file))))
+              trainable_network = trainable_network_t(inference_engine_t(file_t(string_t(network_file))))
               close(network_unit)
             else
               close(network_unit)
@@ -339,9 +339,9 @@ subroutine read_train_write(training_configuration, args, plot_file)
                     maxima = [maxval(pressure_in), maxval(potential_temperature_in), maxval(temperature_in), &
                       maxval(qv_in), maxval(qc_in), maxval(qr_in), maxval(qs_in)] &
                 ) )
-                  associate(activation => training_configuration%differentiable_activation_strategy())
+                  associate(activation => training_configuration%differentiable_activation())
                     associate(residual_network=> string_t(trim(merge("true ", "false", training_configuration%skip_connections()))))
-                      trainable_engine = trainable_engine_t( &
+                      trainable_network = trainable_network_t( &
                         training_configuration,  &
                         perturbation_magnitude = 0.05, &
                         metadata = [ &
@@ -379,7 +379,7 @@ subroutine read_train_write(training_configuration, args, plot_file)
         end associate output_extrema
 
         print *,"Normalizing the remaining input and output tensors"
-        input_output_pairs = trainable_engine%map_to_training_ranges(input_output_pairs)
+        input_output_pairs = trainable_network%map_to_training_ranges(input_output_pairs)
 
         associate( &
           num_pairs => size(input_output_pairs), &
@@ -417,7 +417,7 @@ subroutine read_train_write(training_configuration, args, plot_file)
                 if (size(bins)>1) call shuffle(input_output_pairs) ! set up for stochastic gradient descent
                 mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))]
 
-                call trainable_engine%train(mini_batches, cost, adam, learning_rate)
+                call trainable_network%train(mini_batches, cost, adam, learning_rate)
 
                 associate(average_cost => sum(cost)/size(cost))
                   associate(converged => average_cost <= args%cost_tolerance)
@@ -432,10 +432,8 @@ subroutine read_train_write(training_configuration, args, plot_file)
                         integer net_unit
 
                         open(newunit=net_unit, file=network_file, form='formatted', status='unknown', iostat=io_status, action='write')
-                        associate(inference_engine => trainable_engine%to_inference_engine())
-                          associate(json_file => inference_engine%to_json())
-                            call json_file%write_lines(string_t(network_file))
-                          end associate
+                        associate(json_file => trainable_network%to_json())
+                          call json_file%write_lines(string_t(network_file))
                         end associate
                         close(net_unit)
                       end block
diff --git a/example/learn-addition.F90 b/example/learn-addition.F90
index 0e2efdf3f..b6174650f 100644
--- a/example/learn-addition.F90
+++ b/example/learn-addition.F90
@@ -21,7 +21,7 @@ elemental function y(x_tensor) result(a_tensor)
 program learn_addition
   !! This trains a neural network to learn the following six polynomial functions of its eight inputs.
   use inference_engine_m, only : &
-    inference_engine_t, trainable_engine_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
+    inference_engine_t, trainable_network_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
   use julienne_m, only : string_t, file_t, command_line_t, bin_t
   use assert_m, only : assert, intrinsic_array_t
   use addition_m, only : y
@@ -43,15 +43,15 @@ program learn_addition
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(input_output_pair_t), allocatable :: input_output_pairs(:)
     type(tensor_t), allocatable :: inputs(:), desired_outputs(:)
-    type(trainable_engine_t)  trainable_engine
+    type(trainable_network_t)  trainable_network
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:), random_numbers(:)
 
     call random_init(image_distinct=.true., repeatable=.true.)
-    trainable_engine = perturbed_identity_network(perturbation_magnitude=0.05)
-    call output(trainable_engine%to_inference_engine(), string_t("initial-network.json"))
+    trainable_network = perturbed_identity_network(perturbation_magnitude=0.05)
+    call output(trainable_network, string_t("initial-network.json"))
 
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
+    associate(num_inputs => trainable_network%num_inputs(), num_outputs => trainable_network%num_outputs())
 
       block
         integer i, j
@@ -76,19 +76,18 @@ program learn_addition
           call random_number(random_numbers)
           call shuffle(input_output_pairs)
           mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))]
-          call trainable_engine%train(mini_batches, cost, adam=.true., learning_rate=1.5)
+          call trainable_network%train(mini_batches, cost, adam=.true., learning_rate=1.5)
           print *,sum(cost)/size(cost)
         end do
       end block
 
       block
-        real, parameter :: tolerance = 1.E-06
         integer p
 #if defined _CRAYFTN || __GFORTRAN__
         type(tensor_t), allocatable :: network_outputs(:)
-        network_outputs = trainable_engine%infer(inputs)
+        network_outputs = trainable_network%infer(inputs)
 #else
-        associate(network_outputs => trainable_engine%infer(inputs))
+        associate(network_outputs => trainable_network%infer(inputs))
 #endif
           print "(a,69x,a)","  Outputs", "| Desired outputs"
           do p = 1, num_pairs
@@ -102,14 +101,14 @@ program learn_addition
 
    end associate
 
-   call output(trainable_engine%to_inference_engine(), final_network_file)
+   call output(trainable_network, final_network_file)
 
   end block
 
 contains
 
   subroutine output(inference_engine, file_name)
-    type(inference_engine_t), intent(in) :: inference_engine
+    class(inference_engine_t), intent(in) :: inference_engine
     type(string_t), intent(in) :: file_name
     type(file_t) json_file
     json_file = inference_engine%to_json()
@@ -123,8 +122,8 @@ pure function e(j,n) result(unit_vector)
     unit_vector = real([(merge(1,0,j==k),k=1,n)])
   end function
 
-  function perturbed_identity_network(perturbation_magnitude) result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function perturbed_identity_network(perturbation_magnitude) result(trainable_network)
+    type(trainable_network_t) trainable_network
     real, intent(in) :: perturbation_magnitude
     integer, parameter :: n(*) = [8, 64, 64, 64, 6]
     integer, parameter :: n_max = maxval(n), layers = size(n)
@@ -141,10 +140,10 @@ function perturbed_identity_network(perturbation_magnitude) result(trainable_eng
 
     associate(w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, b => perturbation_magnitude*(b_harvest-0.5)/0.5)
 
-      trainable_engine = trainable_engine_t( &
+      trainable_network = trainable_network_t( inference_engine_t( &
         nodes = n, weights = w, biases = b, metadata = &
           [string_t("Perturbed Identity"), string_t("Damian Rouson"), string_t("2023-09-23"), string_t("relu"), string_t("false")] &
-      )
+      ))
 
     end associate
   end function
diff --git a/example/learn-exponentiation.F90 b/example/learn-exponentiation.F90
index 6b8ba0624..a09b9359d 100644
--- a/example/learn-exponentiation.F90
+++ b/example/learn-exponentiation.F90
@@ -21,7 +21,7 @@ elemental function y(x_tensor) result(a_tensor)
 program learn_exponentiation
   !! This trains a neural network to learn the following six polynomial functions of its eight inputs.
   use inference_engine_m, only : &
-    inference_engine_t, trainable_engine_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
+    inference_engine_t, trainable_network_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
   use julienne_m, only : string_t, file_t, command_line_t, bin_t
   use assert_m, only : assert, intrinsic_array_t
   use exponentiation_m, only : y
@@ -43,15 +43,15 @@ program learn_exponentiation
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(input_output_pair_t), allocatable :: input_output_pairs(:)
     type(tensor_t), allocatable :: inputs(:), desired_outputs(:)
-    type(trainable_engine_t)  trainable_engine
+    type(trainable_network_t)  trainable_network
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:), random_numbers(:)
 
     call random_init(image_distinct=.true., repeatable=.true.)
-    trainable_engine = perturbed_identity_network(perturbation_magnitude=0.05)
-    call output(trainable_engine%to_inference_engine(), string_t("initial-network.json"))
+    trainable_network = perturbed_identity_network(perturbation_magnitude=0.05)
+    call output(trainable_network, string_t("initial-network.json"))
 
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
+    associate(num_inputs => trainable_network%num_inputs(), num_outputs => trainable_network%num_outputs())
 
       block
         integer i, j
@@ -76,19 +76,18 @@ program learn_exponentiation
           call random_number(random_numbers)
           call shuffle(input_output_pairs)
           mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))]
-          call trainable_engine%train(mini_batches, cost, adam=.true., learning_rate=1.5)
+          call trainable_network%train(mini_batches, cost, adam=.true., learning_rate=1.5)
           print *,sum(cost)/size(cost)
         end do
       end block
 
       block
-        real, parameter :: tolerance = 1.E-06
         integer p
 #if defined _CRAYFTN || __GFORTRAN__
         type(tensor_t), allocatable :: network_outputs(:)
-        network_outputs = trainable_engine%infer(inputs)
+        network_outputs = trainable_network%infer(inputs)
 #else
-        associate(network_outputs => trainable_engine%infer(inputs))
+        associate(network_outputs => trainable_network%infer(inputs))
 #endif
           print "(a,69x,a)","  Outputs", "| Desired outputs"
           do p = 1, num_pairs
@@ -102,14 +101,14 @@ program learn_exponentiation
 
    end associate
 
-   call output(trainable_engine%to_inference_engine(), final_network_file)
+   call output(trainable_network, final_network_file)
 
   end block
 
 contains
 
   subroutine output(inference_engine, file_name)
-    type(inference_engine_t), intent(in) :: inference_engine
+    class(inference_engine_t), intent(in) :: inference_engine
     type(string_t), intent(in) :: file_name
     type(file_t) json_file
     json_file = inference_engine%to_json()
@@ -123,8 +122,8 @@ pure function e(j,n) result(unit_vector)
     unit_vector = real([(merge(1,0,j==k),k=1,n)])
   end function
 
-  function perturbed_identity_network(perturbation_magnitude) result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function perturbed_identity_network(perturbation_magnitude) result(trainable_network)
+    type(trainable_network_t) trainable_network
     real, intent(in) :: perturbation_magnitude
     integer, parameter :: n(*) = [8, 64, 64, 64, 6] ! nodes per layer (first layer = input, last layer = output)
     integer, parameter :: n_max = maxval(n), layers = size(n)
@@ -141,10 +140,10 @@ function perturbed_identity_network(perturbation_magnitude) result(trainable_eng
 
     associate(w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, b => perturbation_magnitude*(b_harvest-0.5)/0.5)
 
-      trainable_engine = trainable_engine_t( &
+      trainable_network = trainable_network_t( inference_engine_t( &
         nodes = n, weights = w, biases = b, metadata = &
           [string_t("Perturbed Identity"), string_t("Damian Rouson"), string_t("2023-09-23"), string_t("relu"), string_t("false")] &
-      )
+      ))
 
     end associate
   end function
diff --git a/example/learn-microphysics-procedures.F90 b/example/learn-microphysics-procedures.F90
deleted file mode 100644
index 552bd9ac0..000000000
--- a/example/learn-microphysics-procedures.F90
+++ /dev/null
@@ -1,236 +0,0 @@
-! Copyright (c), The Regents of the University of California
-! Terms of use are as specified in LICENSE.txt
-program learn_microphysics_procedures
-  !! Train a neural network proxies for procedures in the Thompson microphysics model
-  !! in of ICAR (https://github.com/BerkeleyLab/icar).
-  use inference_engine_m, only : &
-    inference_engine_t, trainable_engine_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
-  use julienne_m, only : string_t, file_t, command_line_t, bin_t, csv
-  use assert_m, only : assert, intrinsic_array_t
-  use thompson_tensors_m, only : y, T, p
-  use iso_fortran_env, only : int64, output_unit
-  implicit none
-
-  type(string_t) network_file
-  type(command_line_t) command_line
-  integer(int64) counter_start, counter_end, clock_rate
-
-  network_file = string_t(command_line%flag_value("--output-file"))
-
-  if (len(network_file%string())==0) then
-    error stop new_line('a') // new_line('a') // &
-      'Usage: fpm run learn-microphysics-procedures --profile release --flag "-fopenmp" -- --output-file "<file-name>"'
-  end if
-
-  call system_clock(counter_start, clock_rate)
-
-  block
-    integer, parameter :: max_num_epochs = 10000000, num_mini_batches = 10
-    integer num_pairs ! number of input/output pairs
-
-    type(mini_batch_t), allocatable :: mini_batches(:)
-    type(input_output_pair_t), allocatable :: input_output_pairs(:)
-    type(tensor_t), allocatable :: inputs(:), desired_outputs(:)
-    type(trainable_engine_t)  trainable_engine
-    type(bin_t), allocatable :: bins(:)
-    real, allocatable :: cost(:), random_numbers(:)
-    integer io_status, network_unit, plot_unit
-    integer, parameter :: io_success=0, diagnostics_print_interval = 1000, network_save_interval = 10000
-    integer, parameter :: nodes_per_layer(*) = [2, 72, 2]
-    real, parameter :: cost_tolerance = 1.E-08
-
-    call random_init(image_distinct=.true., repeatable=.true.)
-    open(newunit=network_unit, file=network_file%string(), form='formatted', status='old', iostat=io_status, action='read')
-
-    if (io_status == io_success) then
-      print *,"Reading network from file " // network_file%string()
-      trainable_engine = trainable_engine_t(inference_engine_t(file_t(network_file)))
-      close(network_unit)
-    else
-      close(network_unit)
-      print *,"Initializing a new network"
-      trainable_engine = perturbed_identity_network(perturbation_magnitude=0.05, n = nodes_per_layer)
-    end if
-    call output(trainable_engine%to_inference_engine(), string_t("initial-network.json"))
-
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
-
-      block
-        integer i, j
-        integer, allocatable :: output_sizes(:)
-        inputs = [( [(tensor_t([T(i), p(j)]), j=1,size(p))], i = 1,size(T))]
-        num_pairs = size(inputs)
-        call assert(num_pairs == size(T)*size(p), "train_cloud_microphysics: inputs tensor array complete")
-        desired_outputs = y(inputs)
-        output_sizes = [(size(desired_outputs(i)%values()),i=1,size(desired_outputs))]
-        call assert(all([num_outputs==output_sizes]), "fit-polynomials: # outputs", intrinsic_array_t([num_outputs,output_sizes])) 
-      end block
-
-      input_output_pairs = input_output_pair_t(inputs, desired_outputs)
-
-      block 
-        integer b
-        bins = [(bin_t(num_items=num_pairs, num_bins=num_mini_batches, bin_number=b), b = 1, num_mini_batches)]
-      end block
-
-      block
-        integer e, b, stop_unit, previous_epoch
-        real previous_clock_time
-
-        call open_plot_file_for_appending("cost.plt", plot_unit, previous_epoch, previous_clock_time)
-        print *, "        Epoch | Cost Function| System_Clock | Nodes per Layer"
-        allocate(random_numbers(2:size(input_output_pairs)))
-
-        do e = previous_epoch + 1, previous_epoch + max_num_epochs
-          call random_number(random_numbers)
-          call shuffle(input_output_pairs)
-          mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))]
-          call trainable_engine%train(mini_batches, cost, adam=.true., learning_rate=1.5)
-          call system_clock(counter_end, clock_rate)
-
-          associate( &
-            cost_avg => sum(cost)/size(cost), &
-            cumulative_clock_time => previous_clock_time +  real(counter_end - counter_start) / real(clock_rate), &
-            loop_ending => e == previous_epoch + max_num_epochs &
-          )
-            write_and_exit_if_converged: &
-            if (cost_avg < cost_tolerance) then
-              call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer)
-              call output(trainable_engine%to_inference_engine(), network_file)
-              exit
-            end if write_and_exit_if_converged
-
-            open(newunit=stop_unit, file="stop", form='formatted', status='old', iostat=io_status)
-
-            write_and_exit_if_stop_file_exists: &
-            if (io_status==0) then
-              call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer)
-              call output(trainable_engine%to_inference_engine(), network_file)
-              exit
-            end if write_and_exit_if_stop_file_exists
-
-            if (mod(e,diagnostics_print_interval)==0 .or. loop_ending) &
-              call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer)
-            if (mod(e,network_save_interval)==0 .or. loop_ending) call output(trainable_engine%to_inference_engine(), network_file)
-          end associate
-        end do
-
-        close(plot_unit)
-
-        report_network_performance: &
-        block
-          integer p
-#if defined _CRAYFTN || __GFORTRAN__
-          type(tensor_t), allocatable :: network_outputs(:)
-          network_outputs = trainable_engine%infer(inputs)
-#else
-          associate(network_outputs => trainable_engine%infer(inputs))
-#endif
-            print *," Inputs (normalized)          | Outputs                    | Desired outputs"
-            do p = 1, num_pairs
-              print "(6(G13.5,2x))", inputs(p)%values(), network_outputs(p)%values(), desired_outputs(p)%values()
-            end do
-#if defined _CRAYFTN || __GFORTRAN__
-#else
-          end associate
-#endif
-        end block report_network_performance
-
-      end block
-
-   end associate
-
-   call output(trainable_engine%to_inference_engine(), network_file)
-
-  end block
-
-contains
-
-  subroutine print_diagnostics(plot_file_unit, epoch, cost, clock, nodes)
-     integer, intent(in) :: plot_file_unit, epoch, nodes(:)
-     real, intent(in) :: cost, clock
-
-     write(unit=output_unit, fmt='(3(g13.5,2x))', advance='no') epoch, cost, clock
-     write(unit=output_unit, fmt=csv) nodes
-     write(unit=plot_file_unit, fmt='(3(g13.5,2x))', advance='no') epoch, cost, clock
-     write(unit=plot_file_unit, fmt=csv) nodes
-  end subroutine
-
-  subroutine output(inference_engine, file_name)
-    type(inference_engine_t), intent(in) :: inference_engine
-    type(string_t), intent(in) :: file_name
-    type(file_t) json_file
-    json_file = inference_engine%to_json()
-    call json_file%write_lines(file_name)
-  end subroutine
-
-  pure function e(j,n) result(unit_vector)
-    integer, intent(in) :: j, n
-    integer k
-    real, allocatable :: unit_vector(:)
-    unit_vector = real([(merge(1,0,j==k),k=1,n)])
-  end function
-
-  function perturbed_identity_network(perturbation_magnitude, n) result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
-    real, intent(in) :: perturbation_magnitude
-    integer, intent(in)  :: n(:)
-    integer j, k, l
-    real, allocatable :: identity(:,:,:), w_harvest(:,:,:), b_harvest(:,:)
-
-    associate(n_max => maxval(n), layers => size(n))
-      identity =  reshape( [( [(e(k,n_max), k=1,n_max)], l = 1, layers-1 )], [n_max, n_max, layers-1])
-
-      allocate(w_harvest, mold = identity)
-      allocate(b_harvest(size(identity,1), size(identity,3)))
-
-      call random_number(w_harvest)
-      call random_number(b_harvest)
-
-      associate(w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, b => perturbation_magnitude*(b_harvest-0.5)/0.5)
-
-        trainable_engine = trainable_engine_t( &
-          nodes = n, weights = w, biases = b, metadata = &
-            [string_t("Thompson microphysics procedures"), string_t("Damian Rouson"), string_t("2023-09-23"), string_t("sigmoid"), &
-             string_t("false")] &
-        )
-      end associate
-    end associate
-  end function
-
-  subroutine open_plot_file_for_appending(plot_file_name, plot_unit, previous_epoch, previous_clock)
-    character(len=*), intent(in) :: plot_file_name
-    integer, intent(out) :: plot_unit, previous_epoch
-    real, intent(out) :: previous_clock
-
-    type(file_t) plot_file
-    type(string_t), allocatable :: lines(:)
-    character(len=:), allocatable :: last_line
-    integer io_status
-    integer, parameter :: io_success = 0
-    logical preexisting_plot_file
-    real cost
-
-    inquire(file=plot_file_name, exist=preexisting_plot_file)
-    open(newunit=plot_unit,file="cost.plt",status="unknown",position="append")
-
-    associate(header => "        Epoch | Cost Function| System_Clock | Nodes per Layer")
-      if (.not. preexisting_plot_file) then
-        write(plot_unit,*) header
-        previous_epoch = 0
-        previous_clock = 0
-      else
-        plot_file = file_t(string_t(plot_file_name))
-        lines = plot_file%lines()
-        last_line = lines(size(lines))%string()
-        read(last_line,*, iostat=io_status) previous_epoch, cost, previous_clock
-        if ((io_status /= io_success .and. last_line == header) .or. len(trim(last_line))==0) then
-          previous_epoch = 0
-          previous_clock = 0
-        end if
-      end if
-    end associate
-
-  end subroutine
-
-end program
diff --git a/example/learn-multiplication.F90 b/example/learn-multiplication.F90
index 797335dd9..81d72e88d 100644
--- a/example/learn-multiplication.F90
+++ b/example/learn-multiplication.F90
@@ -21,7 +21,7 @@ elemental function y(x_tensor) result(a_tensor)
 program learn_multiplication
   !! This trains a neural network to learn the following six polynomial functions of its eight inputs.
   use inference_engine_m, only : &
-    inference_engine_t, trainable_engine_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
+    inference_engine_t, trainable_network_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
   use julienne_m, only : string_t, file_t, command_line_t, bin_t
   use assert_m, only : assert, intrinsic_array_t
   use multiply_inputs, only : y
@@ -43,15 +43,15 @@ program learn_multiplication
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(input_output_pair_t), allocatable :: input_output_pairs(:)
     type(tensor_t), allocatable :: inputs(:), desired_outputs(:)
-    type(trainable_engine_t)  trainable_engine
+    type(trainable_network_t)  trainable_network
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:), random_numbers(:)
 
     call random_init(image_distinct=.true., repeatable=.true.)
-    trainable_engine = perturbed_identity_network(perturbation_magnitude=0.05)
-    call output(trainable_engine%to_inference_engine(), string_t("initial-network.json"))
+    trainable_network = perturbed_identity_network(perturbation_magnitude=0.05)
+    call output(trainable_network, string_t("initial-network.json"))
 
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
+    associate(num_inputs => trainable_network%num_inputs(), num_outputs => trainable_network%num_outputs())
 
       block
         integer i, j
@@ -76,19 +76,18 @@ program learn_multiplication
           call random_number(random_numbers)
           call shuffle(input_output_pairs)
           mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))]
-          call trainable_engine%train(mini_batches, cost, adam=.true., learning_rate=1.5)
+          call trainable_network%train(mini_batches, cost, adam=.true., learning_rate=1.5)
           print *,sum(cost)/size(cost)
         end do
       end block
 
       block
-        real, parameter :: tolerance = 1.E-06
         integer p
 #if defined _CRAYFTN || __GFORTRAN__
         type(tensor_t), allocatable :: network_outputs(:)
-        network_outputs = trainable_engine%infer(inputs)
+        network_outputs = trainable_network%infer(inputs)
 #else
-        associate(network_outputs => trainable_engine%infer(inputs))
+        associate(network_outputs => trainable_network%infer(inputs))
 #endif
           print "(a,69x,a)","  Outputs", "| Desired outputs"
           do p = 1, num_pairs
@@ -102,14 +101,14 @@ program learn_multiplication
 
    end associate
 
-   call output(trainable_engine%to_inference_engine(), final_network_file)
+   call output(trainable_network, final_network_file)
 
   end block
 
 contains
 
   subroutine output(inference_engine, file_name)
-    type(inference_engine_t), intent(in) :: inference_engine
+    class(inference_engine_t), intent(in) :: inference_engine
     type(string_t), intent(in) :: file_name
     type(file_t) json_file
     json_file = inference_engine%to_json()
@@ -123,8 +122,8 @@ pure function e(j,n) result(unit_vector)
     unit_vector = real([(merge(1,0,j==k),k=1,n)])
   end function
 
-  function perturbed_identity_network(perturbation_magnitude) result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function perturbed_identity_network(perturbation_magnitude) result(trainable_network)
+    type(trainable_network_t) trainable_network
     real, intent(in) :: perturbation_magnitude
     integer, parameter :: n(*) = [8, 64, 64, 64, 6]
     integer, parameter :: n_max = maxval(n), layers = size(n)
@@ -141,10 +140,10 @@ function perturbed_identity_network(perturbation_magnitude) result(trainable_eng
 
     associate(w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, b => perturbation_magnitude*(b_harvest-0.5)/0.5)
 
-      trainable_engine = trainable_engine_t( &
+      trainable_network = trainable_network_t( inference_engine_t( &
         nodes = n, weights = w, biases = b, metadata = &
           [string_t("Perturbed Identity"), string_t("Damian Rouson"), string_t("2023-09-23"), string_t("relu"), string_t("false")] &
-      )
+      ))
 
     end associate
   end function
diff --git a/example/learn-power-series.F90 b/example/learn-power-series.F90
index c0490669e..e8a46e6dc 100644
--- a/example/learn-power-series.F90
+++ b/example/learn-power-series.F90
@@ -21,7 +21,7 @@ elemental function y(x_in) result(a)
 program learn_power_series
   !! This trains a neural network to learn the following six polynomial functions of its eight inputs.
   use inference_engine_m, only : &
-    inference_engine_t, trainable_engine_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
+    inference_engine_t, trainable_network_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
   use julienne_m, only : string_t, file_t, command_line_t, bin_t
   use assert_m, only : assert, intrinsic_array_t
   use power_series, only : y
@@ -43,15 +43,15 @@ program learn_power_series
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(input_output_pair_t), allocatable :: input_output_pairs(:)
     type(tensor_t), allocatable :: inputs(:), desired_outputs(:)
-    type(trainable_engine_t)  trainable_engine
+    type(trainable_network_t)  trainable_network
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:), random_numbers(:)
 
     call random_init(image_distinct=.true., repeatable=.true.)
-    trainable_engine = perturbed_identity_network(perturbation_magnitude=0.05)
-    call output(trainable_engine%to_inference_engine(), string_t("initial-network.json"))
+    trainable_network = perturbed_identity_network(perturbation_magnitude=0.05)
+    call output(trainable_network, string_t("initial-network.json"))
 
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
+    associate(num_inputs => trainable_network%num_inputs(), num_outputs => trainable_network%num_outputs())
 
       block
         integer i, j
@@ -78,19 +78,18 @@ program learn_power_series
           call random_number(random_numbers)
           call shuffle(input_output_pairs)
           mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))]
-          call trainable_engine%train(mini_batches, cost, adam=.true., learning_rate=1.5)
+          call trainable_network%train(mini_batches, cost, adam=.true., learning_rate=1.5)
           print *,sum(cost)/size(cost)
         end do
       end block
 
       block
-        real, parameter :: tolerance = 1.E-06
         integer p
 #if defined _CRAYFTN || __GFORTRAN__
         type(tensor_t), allocatable :: network_outputs(:)
-        network_outputs = trainable_engine%infer(inputs)
+        network_outputs = trainable_network%infer(inputs)
 #else
-        associate(network_outputs => trainable_engine%infer(inputs))
+        associate(network_outputs => trainable_network%infer(inputs))
 #endif
           print "(a,69x,a)","  Outputs", "| Desired outputs"
           do p = 1, num_pairs
@@ -104,14 +103,14 @@ program learn_power_series
 
    end associate
 
-   call output(trainable_engine%to_inference_engine(), final_network_file)
+   call output(trainable_network, final_network_file)
 
   end block
 
 contains
 
   subroutine output(inference_engine, file_name)
-    type(inference_engine_t), intent(in) :: inference_engine
+    class(inference_engine_t), intent(in) :: inference_engine
     type(string_t), intent(in) :: file_name
     type(file_t) json_file
     json_file = inference_engine%to_json()
@@ -125,8 +124,8 @@ pure function e(j,n) result(unit_vector)
     unit_vector = real([(merge(1,0,j==k),k=1,n)])
   end function
 
-  function perturbed_identity_network(perturbation_magnitude) result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function perturbed_identity_network(perturbation_magnitude) result(trainable_network)
+    type(trainable_network_t) trainable_network
     real, intent(in) :: perturbation_magnitude
     integer, parameter :: n(*) = [8, 196, 196, 196, 196, 6]
     integer, parameter :: n_max = maxval(n), layers = size(n)
@@ -143,10 +142,10 @@ function perturbed_identity_network(perturbation_magnitude) result(trainable_eng
 
     associate(w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, b => perturbation_magnitude*(b_harvest-0.5)/0.5)
 
-      trainable_engine = trainable_engine_t( &
+      trainable_network = trainable_network_t( inference_engine_t( &
         nodes = n, weights = w, biases = b, metadata = &
           [string_t("Perturbed Identity"), string_t("Damian Rouson"), string_t("2023-09-23"), string_t("relu"), string_t("false")] &
-      )
+      ))
 
     end associate
   end function
diff --git a/example/learn-saturated-mixing-ratio.F90 b/example/learn-saturated-mixing-ratio.F90
index 641f33f3e..bf49da58c 100644
--- a/example/learn-saturated-mixing-ratio.F90
+++ b/example/learn-saturated-mixing-ratio.F90
@@ -34,9 +34,9 @@ program train_saturated_mixture_ratio
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:), random_numbers(:)
     integer io_status, network_unit, plot_unit
-    integer, parameter :: io_success=0, diagnostics_print_interval = 100, network_save_interval = 1000
+    integer, parameter :: io_success=0, diagnostics_print_interval = 1000, network_save_interval = 10000
     integer, parameter :: nodes_per_layer(*) = [2, 4, 72, 2, 1]
-    real, parameter :: cost_tolerance = 1.E-06
+    real, parameter :: cost_tolerance = 1.E-08
 
     call random_init(image_distinct=.true., repeatable=.true.)
     open(newunit=network_unit, file=network_file%string(), form='formatted', status='old', iostat=io_status, action='read')
diff --git a/example/supporting-modules/mp_thompson.f90 b/example/supporting-modules/mp_thompson.f90
deleted file mode 100644
index 2e765471d..000000000
--- a/example/supporting-modules/mp_thompson.f90
+++ /dev/null
@@ -1,102 +0,0 @@
-! Copyright (c), The Regents of the University of California
-! Terms of use are as specified in LICENSE.txt
- 
-! MIT License
-! 
-! Copyright (c) 2017 National Center for Atmospheric Research
-! 
-! Permission is hereby granted, free of charge, to any person obtaining a copy
-! of this software and associated documentation files (the "Software"), to deal
-! in the Software without restriction, including without limitation the rights
-! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-! copies of the Software, and to permit persons to whom the Software is
-! furnished to do so, subject to the following conditions:
-! 
-! The above copyright notice and this permission notice shall be included in all
-! copies or substantial portions of the Software.
-! 
-! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-! SOFTWARE.
-
-      MODULE module_mp_thompson
-
-      ! Adapted from https://github.com/BerkeleyLab/icar
-
-      IMPLICIT NONE
-
-      CONTAINS
-
-!+---+-----------------------------------------------------------------+
-! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
-! A FUNCTION OF TEMPERATURE AND PRESSURE
-!
-      REAL FUNCTION RSLF(P,T)
-
-      IMPLICIT NONE
-      REAL, INTENT(IN):: P, T
-      REAL:: ESL,X
-      REAL, PARAMETER:: C0= .611583699E03
-      REAL, PARAMETER:: C1= .444606896E02
-      REAL, PARAMETER:: C2= .143177157E01
-      REAL, PARAMETER:: C3= .264224321E-1
-      REAL, PARAMETER:: C4= .299291081E-3
-      REAL, PARAMETER:: C5= .203154182E-5
-      REAL, PARAMETER:: C6= .702620698E-8
-      REAL, PARAMETER:: C7= .379534310E-11
-      REAL, PARAMETER:: C8=-.321582393E-13
-
-      X=MAX(-80.,T-273.16)
-
-!      ESL=612.2*EXP(17.67*X/(T-29.65))
-      ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
-      RSLF=.622*ESL/(P-ESL)
-
-!    ALTERNATIVE
-!  ; Source: Murphy and Koop, Review of the vapour pressure of ice and
-!             supercooled water for atmospheric applications, Q. J. R.
-!             Meteorol. Soc (2005), 131, pp. 1539-1565.
-!    ESL = EXP(54.842763 - 6763.22 / T - 4.210 * ALOG(T) + 0.000367 * T
-!        + TANH(0.0415 * (T - 218.8)) * (53.878 - 1331.22
-!        / T - 9.44523 * ALOG(T) + 0.014025 * T))
-
-      END FUNCTION RSLF
-!+---+-----------------------------------------------------------------+
-! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
-! FUNCTION OF TEMPERATURE AND PRESSURE
-!
-      REAL FUNCTION RSIF(P,T)
-
-      IMPLICIT NONE
-      REAL, INTENT(IN):: P, T
-      REAL:: ESI,X
-      REAL, PARAMETER:: C0= .609868993E03
-      REAL, PARAMETER:: C1= .499320233E02
-      REAL, PARAMETER:: C2= .184672631E01
-      REAL, PARAMETER:: C3= .402737184E-1
-      REAL, PARAMETER:: C4= .565392987E-3
-      REAL, PARAMETER:: C5= .521693933E-5
-      REAL, PARAMETER:: C6= .307839583E-7
-      REAL, PARAMETER:: C7= .105785160E-9
-      REAL, PARAMETER:: C8= .161444444E-12
-
-      X=MAX(-80.,T-273.16)
-      ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
-      RSIF=.622*ESI/(P-ESI)
-
-!    ALTERNATIVE
-!  ; Source: Murphy and Koop, Review of the vapour pressure of ice and
-!             supercooled water for atmospheric applications, Q. J. R.
-!             Meteorol. Soc (2005), 131, pp. 1539-1565.
-!     ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T)
-
-      END FUNCTION RSIF
-!+---+-----------------------------------------------------------------+
-
-!+---+-----------------------------------------------------------------+
-END MODULE module_mp_thompson
-!+---+-----------------------------------------------------------------+
diff --git a/example/supporting-modules/thompson_tensors_m.f90 b/example/supporting-modules/thompson_tensors_m.f90
deleted file mode 100644
index 548ab5d52..000000000
--- a/example/supporting-modules/thompson_tensors_m.f90
+++ /dev/null
@@ -1,33 +0,0 @@
-! Copyright (c), The Regents of the University of California
-! Terms of use are as specified in LICENSE.txt
-module thompson_tensors_m
-  !! This module supports the program in the file example/learn-microphysics-procedures.f90.
-  use module_mp_thompson, only : rslf, rsif
-  use inference_engine_m, only : tensor_t
-  use assert_m, only : assert
-  implicit none
- 
-  private
-  public :: T, p, y
-
-  real, parameter :: T_min = 236.352524, T_max = 307.610779
-  real, parameter :: p_min = 29671.1348, p_max = 98596.7578
-  integer, parameter :: resolution = 10
-  integer i
-  real, parameter :: T(*) = [(real(i)/real(resolution), i=0,resolution)]
-  real, parameter :: p(*) = [(real(i)/real(resolution), i=0,resolution)]
-
-contains
-
-  elemental impure function y(x_in) result(a)
-    type(tensor_t), intent(in) :: x_in
-    type(tensor_t) a
-    associate(x => x_in%values())
-      call assert(lbound(x,1)==1 .and. ubound(x,1)==2,"y(x) :: sufficient input")
-      associate(temperature => T_min + (T_max - T_min)*x(1), pressure => p_min + (p_max - p_min)*x(2) )
-        a = tensor_t([rslf(pressure, temperature), rsif(pressure, temperature)])
-      end associate
-    end associate
-  end function
- 
-end module
diff --git a/example/train-and-write.F90 b/example/train-and-write.F90
index 04df61a55..7105afde0 100644
--- a/example/train-and-write.F90
+++ b/example/train-and-write.F90
@@ -9,7 +9,7 @@ program train_and_write
   !! The initial condition corresponds to the desired network with all weights and biases perturbed by a random variable
   !! that is uniformly distributed on the range [0,0.1].
   use inference_engine_m, only : &
-    inference_engine_t, trainable_engine_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
+    inference_engine_t, trainable_network_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle
   use julienne_m, only : string_t, file_t, command_line_t, bin_t
   use assert_m, only : assert, intrinsic_array_t
   implicit none
@@ -30,17 +30,17 @@ program train_and_write
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(input_output_pair_t), allocatable :: input_output_pairs(:)
     type(tensor_t), allocatable :: inputs(:)
-    type(trainable_engine_t)  trainable_engine
+    type(trainable_network_t)  trainable_network
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:), random_numbers(:)
 
     call random_init(image_distinct=.true., repeatable=.true.)
-    trainable_engine = perturbed_identity_network(perturbation_magnitude=0.2)
-    call output(trainable_engine%to_inference_engine(), string_t("initial-network.json"))
+    trainable_network = perturbed_identity_network(perturbation_magnitude=0.2)
+    call output(trainable_network, string_t("initial-network.json"))
 
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
+    associate(num_inputs => trainable_network%num_inputs(), num_outputs => trainable_network%num_outputs())
 
-      call assert(num_inputs == num_outputs,"trainable_engine_test_m(identity_mapping): # inputs == # outputs", &
+      call assert(num_inputs == num_outputs,"trainable_network_test_m(identity_mapping): # inputs == # outputs", &
         intrinsic_array_t([num_inputs, num_outputs]) &
       )
       block
@@ -64,7 +64,7 @@ program train_and_write
           call random_number(random_numbers)
           call shuffle(input_output_pairs)
           mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))]
-          call trainable_engine%train(mini_batches, cost, adam=.true., learning_rate=1.5)
+          call trainable_network%train(mini_batches, cost, adam=.true., learning_rate=1.5)
           print *,sum(cost)/size(cost)
         end do
       end block
@@ -74,14 +74,15 @@ program train_and_write
         integer p
 #if defined _CRAYFTN || __GFORTRAN__
         type(tensor_t), allocatable :: network_outputs(:)
-        network_outputs = trainable_engine%infer(inputs)
+        network_outputs = trainable_network%infer(inputs)
 #else
-        associate(network_outputs => trainable_engine%infer(inputs))
+        associate(network_outputs => trainable_network%infer(inputs))
 #endif
-          print "(a,62x,a,53x,a)", " Output", "| Desired outputs", "| Errors"
-
+          print "(a,43x,a,35x,a)", " Output", "| Desired outputs", "| Errors"
+          
           do p = 1, num_pairs
-            print *,network_outputs(p)%values(),"|", inputs(p)%values(), "|",  network_outputs(p)%values() - inputs(p)%values()
+            print "(2(3(G11.5,', '), G11.5, a), 3(G11.5,', '), G11.5)", &
+              network_outputs(p)%values(),"| ", inputs(p)%values(), "| ",  network_outputs(p)%values() - inputs(p)%values()
           end do
 #if defined _CRAYFTN || __GFORTRAN__
 #else
@@ -91,14 +92,14 @@ program train_and_write
 
    end associate
 
-   call output(trainable_engine%to_inference_engine(), final_network_file)
+   call output(trainable_network, final_network_file)
 
   end block
 
 contains
 
   subroutine output(inference_engine, file_name)
-    type(inference_engine_t), intent(in) :: inference_engine
+    class(inference_engine_t), intent(in) :: inference_engine
     type(string_t), intent(in) :: file_name
     type(file_t) json_file
     json_file = inference_engine%to_json()
@@ -111,8 +112,8 @@ pure function e(m,n) result(e_mn)
     e_mn = real(merge(1,0,m==n))
   end function
 
-  function perturbed_identity_network(perturbation_magnitude) result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function perturbed_identity_network(perturbation_magnitude) result(trainable_network)
+    type(trainable_network_t) trainable_network
     real, intent(in) :: perturbation_magnitude
     integer, parameter :: nodes_per_layer(*) = [4, 4, 4, 4]
     integer, parameter :: max_n = maxval(nodes_per_layer), layers = size(nodes_per_layer)
@@ -129,10 +130,10 @@ function perturbed_identity_network(perturbation_magnitude) result(trainable_eng
 
     associate(w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, b => perturbation_magnitude*(b_harvest-0.5)/0.5)
 
-      trainable_engine = trainable_engine_t( &
+      trainable_network = trainable_network_t( inference_engine_t( &
         nodes = nodes_per_layer, weights = w, biases = b, metadata = &
-          [string_t("Perturbed Identity"), string_t("Damian Rouson"), string_t("2023-09-23"), string_t("relu"), string_t("false")] &
-      )
+          [string_t("Perturbed Identity"), string_t("Damian Rouson"), string_t("2024-10-13"), string_t("relu"), string_t("false")] &
+      ))
 
     end associate
 
diff --git a/src/inference_engine/tmp b/src/inference_engine/tmp
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/inference_engine/trainable_engine_m.F90 b/src/inference_engine/trainable_engine_m.F90
deleted file mode 100644
index b7acc6e83..000000000
--- a/src/inference_engine/trainable_engine_m.F90
+++ /dev/null
@@ -1,179 +0,0 @@
-! Copyright (c), The Regents of the University of California
-! Terms of use are as specified in LICENSE.txt
-module trainable_engine_m
-  !! Define an abstraction that supports training a neural network
-  use activation_m, only : activation_t
-  use julienne_string_m, only : string_t
-  use inference_engine_m_, only : inference_engine_t
-  use kind_parameters_m, only : default_real
-  use metadata_m, only : metadata_t
-  use tensor_m, only :  tensor_t
-  use tensor_map_m, only :  tensor_map_t
-  use mini_batch_m, only : mini_batch_t
-  use training_configuration_m, only : training_configuration_t
-  use input_output_pair_m, only : input_output_pair_t
-  implicit none
-
-  private
-  public :: trainable_engine_t
-
-  type trainable_engine_t(k)
-    !! Encapsulate the information needed to perform training
-    integer, kind :: k = default_real
-    type(tensor_map_t(k)), private :: input_map_, output_map_ !! mappings to/from data ranges used during training
-    type(activation_t), private :: activation_
-    type(metadata_t) metadata_       !! metadata_ encapsulates strings for which default-kind suffices
-    real(k), allocatable :: w(:,:,:) !! weights
-    real(k), allocatable :: b(:,:)   !! biases
-    integer, allocatable :: n(:)     !! nodes per layer
-    real(k), allocatable, dimension(:,:) :: a
-    real(k), allocatable, dimension(:,:,:) :: dcdw, vdw, sdw, vdwc, sdwc
-    real(k), allocatable, dimension(:,:) :: z, delta, dcdb, vdb, sdb, vdbc, sdbc
-  contains
-    generic :: train                    => default_real_train
-    procedure, private, non_overridable :: default_real_train
-    generic :: predict => default_real_predict
-    generic :: infer =>   default_real_predict
-    procedure, private :: default_real_predict
-    generic   :: assert_consistent => default_real_assert_consistent
-    procedure, private ::             default_real_assert_consistent
-    generic :: num_layers => default_real_num_layers
-    procedure, private ::    default_real_num_layers
-    generic :: num_inputs => default_real_num_inputs
-    procedure, private, non_overridable :: default_real_num_inputs
-    generic :: num_outputs              => default_real_num_outputs
-    procedure, private ::     default_real_num_outputs
-    generic :: to_inference_engine => default_real_to_inference_engine
-    procedure, private ::             default_real_to_inference_engine
-    generic :: map_to_input_training_range => default_real_map_to_input_training_range
-    procedure, private ::                     default_real_map_to_input_training_range
-    generic :: map_from_input_training_range => default_real_map_from_input_training_range
-    procedure, private ::                       default_real_map_from_input_training_range
-    generic :: map_to_output_training_range => default_real_map_to_output_training_range
-    procedure, private ::                      default_real_map_to_output_training_range
-    generic :: map_from_output_training_range => default_real_map_from_output_training_range
-    procedure, private ::                        default_real_map_from_output_training_range
-    generic :: map_to_training_ranges => default_real_map_to_training_ranges
-    procedure, private ::                default_real_map_to_training_ranges
-  end type
-
-  integer, parameter :: input_layer = 0
-
-  interface trainable_engine_t
-#ifdef __INTEL_COMPILER
-     module function construct_trainable_engine_from_padded_arrays(nodes, weights, biases, metadata, input_map, output_map) &
-#else
-     module function construct_from_padded_arrays( nodes, weights, biases, metadata, input_map, output_map) &
-#endif
-      result(trainable_engine)
-      implicit none
-      integer, intent(in) :: nodes(input_layer:)
-      real, intent(in)  :: weights(:,:,:), biases(:,:)
-      type(string_t), intent(in) :: metadata(:)
-      type(tensor_map_t), intent(in), optional :: input_map, output_map
-      type(trainable_engine_t) trainable_engine
-    end function
-
-    module function construct_from_inference_engine(inference_engine) result(trainable_engine)
-      implicit none
-      type(inference_engine_t), intent(in) :: inference_engine
-      type(trainable_engine_t) trainable_engine
-    end function
-
-    module function perturbed_identity_network(training_configuration, perturbation_magnitude, metadata, input_map, output_map)&
-      result(trainable_engine)
-      implicit none
-      type(training_configuration_t), intent(in) :: training_configuration
-      type(string_t), intent(in) :: metadata(:)
-      real, intent(in) :: perturbation_magnitude
-      type(tensor_map_t) input_map, output_map
-      type(trainable_engine_t) trainable_engine
-    end function
-
-  end interface
-
-  interface
-
-    pure module subroutine default_real_assert_consistent(self)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-    end subroutine
-
-    pure module subroutine default_real_train(self, mini_batches_arr, cost, adam, learning_rate)
-      implicit none
-      class(trainable_engine_t), intent(inout) :: self
-      type(mini_batch_t), intent(in) :: mini_batches_arr(:)
-      real, intent(out), allocatable, optional :: cost(:)
-      logical, intent(in) :: adam
-      real, intent(in) :: learning_rate
-    end subroutine
-
-    elemental module function default_real_predict(self, inputs) result(outputs)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      type(tensor_t), intent(in) :: inputs
-      type(tensor_t) outputs
-    end function
-    
-    elemental module function default_real_num_inputs(self) result(n_in)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      integer n_in
-    end function
-
-    elemental module function default_real_num_outputs(self) result(n_out)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      integer n_out
-    end function
-
-    elemental module function default_real_num_layers(self) result(n_layers)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      integer n_layers
-    end function
-
-    module function default_real_to_inference_engine(self) result(inference_engine)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      type(inference_engine_t) :: inference_engine
-    end function
-
-    elemental module function default_real_map_to_input_training_range(self, tensor) result(normalized_tensor)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      type(tensor_t), intent(in) :: tensor
-      type(tensor_t) normalized_tensor
-    end function
-
-    elemental module function default_real_map_from_input_training_range(self, tensor) result(unnormalized_tensor)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      type(tensor_t), intent(in) :: tensor
-      type(tensor_t) unnormalized_tensor
-    end function
-
-    elemental module function default_real_map_to_output_training_range(self, tensor) result(normalized_tensor)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      type(tensor_t), intent(in) :: tensor
-      type(tensor_t) normalized_tensor
-    end function
-
-    elemental module function default_real_map_from_output_training_range(self, tensor) result(unnormalized_tensor)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      type(tensor_t), intent(in) :: tensor
-      type(tensor_t) unnormalized_tensor
-    end function
-
-    elemental module function default_real_map_to_training_ranges(self, input_output_pair) result(normalized_input_output_pair)
-      implicit none
-      class(trainable_engine_t), intent(in) :: self
-      type(input_output_pair_t), intent(in) :: input_output_pair
-      type(input_output_pair_t) normalized_input_output_pair
-    end function
-
-  end interface
-
-end module trainable_engine_m
diff --git a/src/inference_engine/trainable_engine_s.F90 b/src/inference_engine/trainable_engine_s.F90
deleted file mode 100644
index f72096fe6..000000000
--- a/src/inference_engine/trainable_engine_s.F90
+++ /dev/null
@@ -1,416 +0,0 @@
-! Copyright (c), The Regents of the University of California
-! Terms of use are as specified in LICENSE.txt
-
-#include "language-support.F90"
-
-submodule(trainable_engine_m) trainable_engine_s
-  use assert_m, only : assert
-  use intrinsic_array_m, only : intrinsic_array_t
-  use tensor_m, only : tensor_t
-#ifdef _CRAYFTN
-  use input_output_pair_m, only : input_output_pair_t
-#endif
-  implicit none
-
-  integer, parameter :: input_layer = 0
-
-contains
-
-  module procedure default_real_num_inputs
-    n_in = self%n(input_layer)
-  end procedure
-
-  module procedure default_real_num_layers
-    n_layers = size(self%n,1)
-  end procedure
-
-  module procedure default_real_num_outputs
-    n_out = self%n(ubound(self%n,1))
-  end procedure
-
-  module procedure construct_from_inference_engine
-
-#ifndef _CRAYFTN
-    associate(exchange => inference_engine%to_exchange())
-#else
-    use inference_engine_m_, only: exchange_t
-    type(exchange_t) exchange
-    exchange = inference_engine%to_exchange()
-#endif
-      trainable_engine%input_map_ = exchange%input_map_
-      trainable_engine%output_map_ = exchange%output_map_
-      trainable_engine%metadata_ = exchange%metadata_
-      trainable_engine%w = exchange%weights_
-      trainable_engine%b = exchange%biases_
-      trainable_engine%n = exchange%nodes_
-      trainable_engine%activation_ = exchange%activation_
-#ifndef _CRAYFTN
-    end associate
-#endif
-
-  end procedure
-
-  module procedure default_real_assert_consistent
-
-    associate( &
-      fully_allocated=>[allocated(self%w),allocated(self%b),allocated(self%n)] &
-    )
-      call assert(all(fully_allocated),"trainable_engine_s(assert_consistent): fully_allocated",intrinsic_array_t(fully_allocated))
-    end associate
-
-    associate(max_width => maxval(self%n), component_dims => [size(self%b,1), size(self%w,1), size(self%w,2)])
-      call assert(all(component_dims == max_width), "trainable_engine_s(assert_consistent): conformable arrays", &
-        intrinsic_array_t([max_width,component_dims]))
-    end associate
-
-    call assert(lbound(self%n,1)==input_layer, "trainable_engine_s(assert_consistent): n base subsscript", lbound(self%n,1))
-
-  end procedure
-
-  module procedure default_real_predict
-
-    real, allocatable :: a(:,:)
-    integer l
-
-    call self%assert_consistent
-
-    associate(w => self%w, b => self%b, n => self%n, output_layer => ubound(self%n,1))
-
-      allocate(a(maxval(n), input_layer:output_layer)) ! Activations
-
-#ifndef _CRAYFTN
-      associate(normalized_inputs => self%input_map_%map_to_training_range(inputs))
-#else
-      block
-        type(tensor_t) normalized_inputs
-        normalized_inputs = self%input_map_%map_to_training_range(inputs)
-#endif
-        a(1:n(input_layer),input_layer) = normalized_inputs%values()
-#ifndef _CRAYFTN
-      end associate
-#else
-      end block
-#endif
-
-      feed_forward: &
-      do l = 1,output_layer
-        associate(z=>matmul(w(1:n(l),1:n(l-1),l), a(1:n(l-1),l-1)) + b(1:n(l),l))
-          a(1:n(l),l) = self%activation_%evaluate(z)
-        end associate
-      end do feed_forward
- 
-      associate(normalized_outputs => tensor_t(a(1:n(output_layer), output_layer)))
-        outputs = self%output_map_%map_from_training_range(normalized_outputs)
-      end associate
-
-    end associate
-
-  end procedure
-
-  module procedure default_real_train
-    integer l, batch, mini_batch_size, pair
-    type(tensor_t), allocatable :: inputs(:), expected_outputs(:)
-
-    call self%assert_consistent
-
-    if (.not. allocated(self%dcdw)) allocate(self%dcdw, mold=self%w) ! Gradient of cost function with respect to weights
-    if (.not. allocated(self%vdw)) allocate(self%vdw, mold=self%w) 
-    if (.not. allocated(self%sdw)) allocate(self%sdw, mold=self%w) 
-    if (.not. allocated(self%vdwc)) allocate(self%vdwc,  mold=self%w) 
-    if (.not. allocated(self%sdwc)) allocate(self%sdwc,  mold=self%w) 
-
-    if (.not. allocated(self%dcdb)) allocate(self%dcdb,  mold=self%b) ! Gradient of cost function with respect with biases
-    if (.not. allocated(self%vdb)) allocate(self%vdb,   mold=self%b) 
-    if (.not. allocated(self%sdb)) allocate(self%sdb,   mold=self%b) 
-    if (.not. allocated(self%vdbc)) allocate(self%vdbc,  mold=self%b) 
-    if (.not. allocated(self%sdbc)) allocate(self%sdbc,  mold=self%b) 
-
-    associate(output_layer => ubound(self%n,1))
-      
-#if F2023_LOCALITY || F2018_LOCALITY
-      if (.not. allocated(self%z)) allocate(self%z,  mold=self%b) ! z-values: Sum z_j^l = w_jk^{l} a_k^{l-1} + b_j^l
-      if (.not. allocated(self%delta)) allocate(self%delta, mold=self%b)
-      if (.not. allocated(self%a)) allocate(self%a(maxval(self%n), input_layer:output_layer)) ! Activations
-#endif
-
-      associate( &
-        a => self%a, dcdw => self%dcdw, vdw => self%vdw, sdw => self%sdw, vdwc => self%vdwc, sdwc => self%sdwc, &
-        z => self%z, delta => self%delta, dcdb => self%dcdb, vdb => self%vdb, sdb => self%sdb, vdbc => self%vdbc, sdbc=> self%sdbc &
-      )
-        vdw = 0.; sdw = 1.; vdb = 0.; sdb = 1.
-
-        associate(w => self%w, b => self%b, n => self%n, num_mini_batches => size(mini_batches_arr))
-
-          if (present(cost)) allocate(cost(num_mini_batches))
-        
-          iterate_across_batches: &
-          do batch = 1, num_mini_batches
-
-            dcdw = 0.; dcdb = 0.
-
-#ifndef _CRAYFTN
-            associate(input_output_pairs => mini_batches_arr(batch)%input_output_pairs())
-#else
-            block
-              type(input_output_pair_t), allocatable :: input_output_pairs(:)
-              input_output_pairs = mini_batches_arr(batch)%input_output_pairs()
-#endif  
-              inputs = input_output_pairs%inputs()
-              expected_outputs = input_output_pairs%expected_outputs()
-              mini_batch_size = size(input_output_pairs)
-#ifndef _CRAYFTN
-            end associate
-#else
-            end block
-#endif  
-            block
-              real, allocatable :: pair_cost(:)
-              if (present(cost)) allocate(pair_cost(mini_batch_size))
-
-#if F2023_LOCALITY
-              iterate_through_batch: &
-              do concurrent (pair = 1:mini_batch_size) local(a,z,delta) reduce(+: dcdb, dcdw)
-
-#elif F2018_LOCALITY
-
-              reduce_gradients: &
-              block
-                real reduce_dcdb(size(dcdb,1),size(dcdb,2),mini_batch_size)
-                real reduce_dcdw(size(dcdw,1),size(dcdw,2),size(dcdw,3),mini_batch_size)
-                reduce_dcdb = 0.
-                reduce_dcdw = 0.
-
-              iterate_through_batch: &
-              do concurrent (pair = 1:mini_batch_size) local(a,z,delta)
-
-#else
-
-              reduce_gradients: &
-              block
-                real reduce_dcdb(size(dcdb,1),size(dcdb,2),mini_batch_size)
-                real reduce_dcdw(size(dcdw,1),size(dcdw,2),size(dcdw,3),mini_batch_size)
-                reduce_dcdb = 0.
-                reduce_dcdw = 0.
-              
-              iterate_through_batch: &
-              do concurrent (pair = 1:mini_batch_size)
-
-                iteration: &
-                block
-
-                real a(maxval(self%n), input_layer:output_layer) ! Activations
-                real z(size(b,1),size(b,2)), delta(size(b,1),size(b,2))
-#endif
-
-                a(1:self%num_inputs(), input_layer) = inputs(pair)%values()
-
-                feed_forward: &
-                do l = 1,output_layer
-                  z(1:n(l),l) = matmul(w(1:n(l),1:n(l-1),l), a(1:n(l-1),l-1)) + b(1:n(l),l)
-                  a(1:n(l),l) = self%activation_%evaluate(z(1:n(l),l))
-                end do feed_forward
-
-                associate(y => expected_outputs(pair)%values())
-                  if (present(cost)) pair_cost(pair) = sum((y(1:n(output_layer))-a(1:n(output_layer),output_layer))**2)
-              
-                  delta(1:n(output_layer),output_layer) = (a(1:n(output_layer),output_layer) - y(1:n(output_layer))) &
-                                                         * self%activation_%differentiate(z(1:n(output_layer),output_layer))
-                end associate
-                
-                associate(n_hidden => self%num_layers()-2)
-                  back_propagate_error: &
-                  do l = n_hidden,1,-1
-                    delta(1:n(l),l) = matmul(transpose(w(1:n(l+1),1:n(l),l+1)), delta(1:n(l+1),l+1)) &
-                                     * self%activation_%differentiate(z(1:n(l),l))
-                  end do back_propagate_error
-                end associate
-
-                block
-                  integer j
-                  sum_gradients: &
-                  do l = 1,output_layer
-#if F2023_LOCALITY
-                    dcdb(1:n(l),l) = dcdb(1:n(l),l) + delta(1:n(l),l)
-                    do concurrent(j = 1:n(l)) reduce(+: dcdw)
-                      dcdw(j,1:n(l-1),l) = dcdw(j,1:n(l-1),l) + a(1:n(l-1),l-1)*delta(j,l)
-                    end do
-#else
-                    reduce_dcdb(1:n(l),l,pair) = reduce_dcdb(1:n(l),l,pair) + delta(1:n(l),l)
-                    do j = 1,n(l)
-                      reduce_dcdw(j,1:n(l-1),l,pair) = reduce_dcdw(j,1:n(l-1),l,pair) + a(1:n(l-1),l-1)*delta(j,l)
-                    end do
-#endif
-                  end do sum_gradients
-                end block
-    
-#if F2023_LOCALITY
-                end do iterate_through_batch
-#elif F2018_LOCALITY
-
-                end do iterate_through_batch
-                dcdb = sum(reduce_dcdb,dim=3)
-                dcdw = sum(reduce_dcdw,dim=4)
-
-                end block reduce_gradients
-#else
-                end block iteration
-                end do iterate_through_batch
-                dcdb = sum(reduce_dcdb,dim=3)
-                dcdw = sum(reduce_dcdw,dim=4)
-  
-                end block reduce_gradients
-#endif
-
-              if (present(cost)) cost(batch) = sum(pair_cost)/(2*mini_batch_size)
-            end block
-          
-            if (adam) then
-              block
-                ! Adam parameters  
-                real, parameter :: beta(*) = [.9, .999]
-                real, parameter :: obeta(*) = [1.- beta(1), 1.- beta(2)]
-                real, parameter :: epsilon = 1.E-08
-
-                associate(alpha => learning_rate)
-                  adam_adjust_weights_and_biases: &
-                  do concurrent(l = 1:output_layer)
-                    dcdw(1:n(l),1:n(l-1),l) = dcdw(1:n(l),1:n(l-1),l)/(mini_batch_size)
-                    vdw(1:n(l),1:n(l-1),l)  = beta(1)*vdw(1:n(l),1:n(l-1),l) + obeta(1)*dcdw(1:n(l),1:n(l-1),l)
-                    sdw (1:n(l),1:n(l-1),l) = beta(2)*sdw(1:n(l),1:n(l-1),l) + obeta(2)*(dcdw(1:n(l),1:n(l-1),l)**2)
-                    vdwc(1:n(l),1:n(l-1),l) = vdw(1:n(l),1:n(l-1),l)/(1.- beta(1)**num_mini_batches)
-                    sdwc(1:n(l),1:n(l-1),l) = sdw(1:n(l),1:n(l-1),l)/(1.- beta(2)**num_mini_batches)
-                    w(1:n(l),1:n(l-1),l) = w(1:n(l),1:n(l-1),l) &
-                      - alpha*vdwc(1:n(l),1:n(l-1),l)/(sqrt(sdwc(1:n(l),1:n(l-1),l))+epsilon) ! Adjust weights
-
-                    dcdb(1:n(l),l) = dcdb(1:n(l),l)/mini_batch_size
-                    vdb(1:n(l),l) = beta(1)*vdb(1:n(l),l) + obeta(1)*dcdb(1:n(l),l)
-                    sdb(1:n(l),l) = beta(2)*sdb(1:n(l),l) + obeta(2)*(dcdb(1:n(l),l)**2)
-                    vdbc(1:n(l),l) = vdb(1:n(l),l)/(1. - beta(1)**num_mini_batches)
-                    sdbc(1:n(l),l) = sdb(1:n(l),l)/(1. - beta(2)**num_mini_batches)
-                    b(1:n(l),l) = b(1:n(l),l) - alpha*vdbc(1:n(l),l)/(sqrt(sdbc(1:n(l),l))+epsilon) ! Adjust weights
-                  end do adam_adjust_weights_and_biases
-                end associate
-              end block
-            else
-              associate(eta => learning_rate)
-                adjust_weights_and_biases: &
-                do concurrent(l = 1:output_layer)
-                  dcdb(1:n(l),l) = dcdb(1:n(l),l)/mini_batch_size
-                  b(1:n(l),l) = b(1:n(l),l) - eta*dcdb(1:n(l),l) ! Adjust biases
-                  dcdw(1:n(l),1:n(l-1),l) = dcdw(1:n(l),1:n(l-1),l)/mini_batch_size
-                  w(1:n(l),1:n(l-1),l) = w(1:n(l),1:n(l-1),l) - eta*dcdw(1:n(l),1:n(l-1),l) ! Adjust weights
-                end do adjust_weights_and_biases
-              end associate
-            end if
-          end do iterate_across_batches
-        end associate
-      end associate
-    end associate
-  end procedure
-
-#ifdef __INTEL_COMPILER
-  module procedure construct_trainable_engine_from_padded_arrays
-#else
-  module procedure construct_from_padded_arrays
-#endif
-
-    trainable_engine%metadata_ = metadata_t(metadata(1),metadata(2),metadata(3),metadata(4),metadata(5))
-    trainable_engine%n = nodes
-    trainable_engine%w = weights
-    trainable_engine%b = biases
-    trainable_engine%activation_ = activation_t(metadata(4)%string())
-
-    block 
-      integer i
-
-      if (present(input_map)) then
-         trainable_engine%input_map_ = input_map
-      else
-        associate(num_inputs => nodes(lbound(nodes,1)))
-          trainable_engine%input_map_ = tensor_map_t("inputs", minima=[(0., i=1,num_inputs)], maxima=[(1., i=1,num_inputs)])
-        end associate
-      end if
-
-      if (present(output_map)) then
-         trainable_engine%output_map_ = output_map
-      else
-        associate(num_outputs => nodes(ubound(nodes,1)))
-          trainable_engine%output_map_ = tensor_map_t("outputs", minima=[(0., i=1,num_outputs)], maxima=[(1., i=1,num_outputs)])
-        end associate
-      end if
-    end block
-
-    call trainable_engine%assert_consistent
-  end procedure
-
-  module procedure default_real_to_inference_engine
-    inference_engine = inference_engine_t(self%metadata_%strings(), self%w, self%b, self%n, self%input_map_, self%output_map_)
-  end procedure
-
-  module procedure perturbed_identity_network
-
-    integer k, l
-    real, allocatable :: identity(:,:,:), w_harvest(:,:,:), b_harvest(:,:)
-
-    associate(n=>training_configuration%nodes_per_layer())
-      associate(n_max => maxval(n), layers => size(n))
-
-        identity = reshape( [( [(e(k,n_max), k=1,n_max)], l = 1, layers-1 )], [n_max, n_max, layers-1])
-        allocate(w_harvest, mold = identity)
-        allocate(b_harvest(size(identity,1), size(identity,3)))
-        call random_number(w_harvest)
-        call random_number(b_harvest)
-
-        associate( &
-          w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, &
-          b => perturbation_magnitude*(b_harvest-0.5)/0.5 &
-        )
-          trainable_engine = trainable_engine_t( &
-            nodes = n, weights = w, biases = b, metadata = metadata, input_map = input_map, output_map = output_map &
-          )
-        end associate
-      end associate
-    end associate
-
-  contains
-
-    pure function e(j,n) result(unit_vector)
-      integer, intent(in) :: j, n
-      integer k
-      real, allocatable :: unit_vector(:)
-      unit_vector = real([(merge(1,0,j==k),k=1,n)])
-    end function
-
-  end procedure
-
-  module procedure default_real_map_to_training_ranges
-    associate( &
-      inputs => input_output_pair%inputs(), &
-      expected_outputs => input_output_pair%expected_outputs() &
-    )
-      associate( &
-         normalized_inputs => self%input_map_%map_to_training_range(inputs), &
-         normalized_outputs => self%output_map_%map_to_training_range(expected_outputs) &
-      )
-        normalized_input_output_pair = input_output_pair_t(normalized_inputs, normalized_outputs)
-      end associate
-    end associate
-  end procedure
-
-  module procedure default_real_map_to_input_training_range
-    normalized_tensor = self%input_map_%map_to_training_range(tensor)
-  end procedure
-
-  module procedure default_real_map_from_input_training_range
-    unnormalized_tensor = self%input_map_%map_from_training_range(tensor)
-  end procedure
-  
-  module procedure default_real_map_to_output_training_range
-    normalized_tensor = self%output_map_%map_to_training_range(tensor)
-  end procedure
-
-  module procedure default_real_map_from_output_training_range
-    unnormalized_tensor = self%output_map_%map_from_training_range(tensor)
-  end procedure
-
-end submodule trainable_engine_s
diff --git a/src/inference_engine/trainable_network_m.f90 b/src/inference_engine/trainable_network_m.f90
index 3338f0f69..86380eb0f 100644
--- a/src/inference_engine/trainable_network_m.f90
+++ b/src/inference_engine/trainable_network_m.f90
@@ -1,7 +1,11 @@
 module trainable_network_m
   use inference_engine_m_, only : inference_engine_t, workspace_t
-  use mini_batch_m, only : mini_batch_t
+  use input_output_pair_m, only : input_output_pair_t
+  use julienne_m, only : string_t
   use kind_parameters_m, only : default_real
+  use mini_batch_m, only : mini_batch_t
+  use training_configuration_m, only : training_configuration_t
+  use tensor_map_m, only : tensor_map_t
   implicit none
 
   private
@@ -13,7 +17,9 @@ module trainable_network_m
     type(workspace_t), private :: workspace_
   contains
     generic :: train           => default_real_train                   
-    procedure, non_overridable :: default_real_train
+    procedure, private, non_overridable :: default_real_train
+    generic ::   map_to_training_ranges => default_real_map_to_training_ranges
+    procedure, private, non_overridable :: default_real_map_to_training_ranges
   end type
 
   interface trainable_network_t 
@@ -23,7 +29,17 @@ pure module function default_real_network(inference_engine) result(trainable_net
       type(inference_engine_t), intent(in) :: inference_engine
       type(trainable_network_t) trainable_network
     end function 
- 
+
+    module function perturbed_identity_network(training_configuration, perturbation_magnitude, metadata, input_map, output_map) &
+      result(trainable_network)
+      implicit none
+      type(training_configuration_t), intent(in) :: training_configuration
+      type(string_t), intent(in) :: metadata(:)
+      real, intent(in) :: perturbation_magnitude
+      type(tensor_map_t) input_map, output_map
+      type(trainable_network_t) trainable_network
+    end function
+
   end interface
 
   interface
@@ -37,6 +53,13 @@ pure module subroutine default_real_train(self, mini_batches_arr, cost, adam, le
       real, intent(in) :: learning_rate
     end subroutine
 
+    elemental module function default_real_map_to_training_ranges(self, input_output_pair) result(normalized_input_output_pair)
+      implicit none
+      class(trainable_network_t), intent(in) :: self
+      type(input_output_pair_t), intent(in) :: input_output_pair
+      type(input_output_pair_t) normalized_input_output_pair
+    end function
+
   end interface
 
 end module trainable_network_m
diff --git a/src/inference_engine/trainable_network_s.f90 b/src/inference_engine/trainable_network_s.f90
index d6e1d83a6..a4e3caa6c 100644
--- a/src/inference_engine/trainable_network_s.f90
+++ b/src/inference_engine/trainable_network_s.f90
@@ -9,7 +9,53 @@
   end procedure
 
   module procedure default_real_train
-    call self%inference_engine_t%default_real_learn(mini_batches_arr, cost, adam, learning_rate, self%workspace_)
+    call self%learn(mini_batches_arr, cost, adam, learning_rate, self%workspace_)
   end procedure
 
+  module procedure default_real_map_to_training_ranges
+    associate(inputs => input_output_pair%inputs(), expected_outputs => input_output_pair%expected_outputs())
+      associate( &
+         normalized_inputs  => self%input_map_%map_to_training_range(inputs), &
+         normalized_outputs => self%output_map_%map_to_training_range(expected_outputs) &
+      )
+        normalized_input_output_pair = input_output_pair_t(normalized_inputs, normalized_outputs)
+      end associate
+    end associate
+  end procedure
+
+  module procedure perturbed_identity_network
+
+    integer k, l
+    real, allocatable :: identity(:,:,:), w_harvest(:,:,:), b_harvest(:,:)
+
+    associate(n=>training_configuration%nodes_per_layer())
+      associate(n_max => maxval(n), layers => size(n))
+
+        identity = reshape( [( [(e(k,n_max), k=1,n_max)], l = 1, layers-1 )], [n_max, n_max, layers-1])
+        allocate(w_harvest, mold = identity)
+        allocate(b_harvest(size(identity,1), size(identity,3)))
+        call random_number(w_harvest)
+        call random_number(b_harvest)
+
+        associate( &
+          w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, &
+          b => perturbation_magnitude*(b_harvest-0.5)/0.5 &
+        )
+          trainable_network = trainable_network_t( &
+            inference_engine_t(nodes=n, weights=w, biases=b, metadata=metadata, input_map=input_map, output_map=output_map) &
+          )
+        end associate
+      end associate
+    end associate
+
+  contains
+
+    pure function e(j,n) result(unit_vector)
+      integer, intent(in) :: j, n
+      integer k
+      real, allocatable :: unit_vector(:)
+      unit_vector = real([(merge(1,0,j==k),k=1,n)])
+    end function
+
+  end procedure
 end submodule trainable_network_s
diff --git a/src/inference_engine_m.f90 b/src/inference_engine_m.f90
index 5b6820d1e..45a803ab7 100644
--- a/src/inference_engine_m.f90
+++ b/src/inference_engine_m.f90
@@ -14,7 +14,6 @@ module inference_engine_m
  use tensor_m, only : tensor_t
  use tensor_map_m, only : tensor_map_t
  use trainable_network_m, only : trainable_network_t
- use trainable_engine_m, only : trainable_engine_t
  use training_configuration_m, only : training_configuration_t
  use ubounds_m, only : ubounds_t
  implicit none
diff --git a/test/main.F90 b/test/main.F90
index 1b7d8d4ec..a21444d24 100644
--- a/test/main.F90
+++ b/test/main.F90
@@ -3,7 +3,7 @@
 program main
   use inference_engine_test_m, only : inference_engine_test_t
   use asymmetric_engine_test_m, only : asymmetric_engine_test_t
-  use trainable_engine_test_m, only : trainable_engine_test_t
+  use trainable_network_test_m, only : trainable_network_test_t
   use metadata_test_m, only : metadata_test_t
   use hyperparameters_test_m, only : hyperparameters_test_t
   use network_configuration_test_m, only : network_configuration_test_t
@@ -15,7 +15,7 @@ program main
 
   type(inference_engine_test_t) inference_engine_test
   type(asymmetric_engine_test_t) asymmetric_engine_test
-  type(trainable_engine_test_t) trainable_engine_test
+  type(trainable_network_test_t) trainable_network_test
   type(hyperparameters_test_t) hyperparameters_test
   type(metadata_test_t) metadata_test
   type(network_configuration_test_t) network_configuration_test
@@ -49,7 +49,7 @@ program main
   call tensor_test%report(passes, tests)
   call asymmetric_engine_test%report(passes, tests)
   call inference_engine_test%report(passes, tests)
-  call trainable_engine_test%report(passes, tests)
+  call trainable_network_test%report(passes, tests)
   call cpu_time(t_finish)
 
   print *
diff --git a/test/trainable_engine_test_m.F90 b/test/trainable_network_test_m.F90
similarity index 86%
rename from test/trainable_engine_test_m.F90
rename to test/trainable_network_test_m.F90
index 929d40ea6..92518d74d 100644
--- a/test/trainable_engine_test_m.F90
+++ b/test/trainable_network_test_m.F90
@@ -1,6 +1,6 @@
 ! Copyright (c), The Regents of the University of California
 ! Terms of use are as specified in LICENSE.txt
-module trainable_engine_test_m
+module trainable_network_test_m
   !! Define inference tests and procedures required for reporting results
 
   ! External dependencies
@@ -12,11 +12,11 @@ module trainable_engine_test_m
 #endif
 
   ! Internal dependencies
-  use inference_engine_m, only : trainable_engine_t, tensor_t, input_output_pair_t, mini_batch_t, shuffle
+  use inference_engine_m, only : trainable_network_t, inference_engine_t, tensor_t, input_output_pair_t, mini_batch_t, shuffle
   implicit none
 
   private
-  public :: trainable_engine_test_t
+  public :: trainable_network_test_t
 
   type, extends(vector_function_strategy_t) :: and_gate_test_function_t
   contains
@@ -38,7 +38,7 @@ module trainable_engine_test_m
     procedure, nopass :: vector_function => xor_gate_with_random_weights
   end type
     
-  type, extends(test_t) :: trainable_engine_test_t
+  type, extends(test_t) :: trainable_network_test_t
   contains
     procedure, nopass :: subject
     procedure, nopass :: results
@@ -60,7 +60,7 @@ function map_i(inputs) result(expected_outputs)
 
   pure function subject() result(specimen)
     character(len=:), allocatable :: specimen
-    specimen = "A trainable_engine_t" 
+    specimen = "A trainable_network_t object"
   end function
 
   function results() result(test_results)
@@ -163,7 +163,7 @@ subroutine print_truth_table(gate_name, gate_function_ptr, test_inputs, actual_o
     integer i
 
     call assert( size(test_inputs) == size(actual_outputs), &
-      "trainable_engine_test_m(print_truth_table): size(test_inputs) == size(actual_outputs)")
+      "trainable_network_test_m(print_truth_table): size(test_inputs) == size(actual_outputs)")
 
     print *,"_______" // gate_name // "_______"
 
@@ -173,8 +173,8 @@ subroutine print_truth_table(gate_name, gate_function_ptr, test_inputs, actual_o
     end do
   end subroutine
 
-  function two_zeroed_hidden_layers() result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function two_zeroed_hidden_layers() result(trainable_network)
+    type(trainable_network_t) trainable_network
     integer, parameter :: inputs = 2, outputs = 1, hidden = 3 ! number of neurons in input, output, and hidden layers
     integer, parameter :: neurons(*) = [inputs, hidden, hidden, outputs] ! neurons per layer
     integer, parameter :: max_neurons = maxval(neurons), layers=size(neurons) ! max layer width, number of layers
@@ -183,14 +183,14 @@ function two_zeroed_hidden_layers() result(trainable_engine)
     w = 0.
     b = 0.
 
-    trainable_engine = trainable_engine_t( &
-      nodes = neurons, weights = w, biases = b, metadata = &
-      [string_t("2-hide|3-wide"), string_t("Damian Rouson"), string_t("2023-06-30"), string_t("sigmoid"), string_t("false")] &
-    )   
+    trainable_network = trainable_network_t( inference_engine_t( &
+      nodes = neurons, weights = w, biases = b &
+     ,metadata = [string_t("2-hide|3-wide"), string_t("Rouson"), string_t("2023-06-30"), string_t("sigmoid"), string_t("false")] &
+    ))
   end function
 
-  function two_random_hidden_layers() result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function two_random_hidden_layers() result(trainable_network)
+    type(trainable_network_t) trainable_network
     integer, parameter :: inputs = 2, outputs = 1, hidden = 3 ! number of neurons in input, output, and hidden layers
     integer, parameter :: neurons(*) = [inputs, hidden, hidden, outputs] ! neurons per layer
     integer, parameter :: max_neurons = maxval(neurons), layers=size(neurons) ! max layer width, number of layers
@@ -201,10 +201,10 @@ function two_random_hidden_layers() result(trainable_engine)
     b = 2*b
     w = 2*w
 
-    trainable_engine = trainable_engine_t( &
-      nodes = neurons, weights = w, biases = b, metadata = &
-      [string_t("2-hide|3-wide"), string_t("Damian Rouson"), string_t("2023-06-30"), string_t("sigmoid"), string_t("false")] &
-    )   
+    trainable_network = trainable_network_t( inference_engine_t( &
+      nodes = neurons, weights = w, biases = b &
+      ,metadata = [string_t("2-hide|3-wide"), string_t("Rouson"), string_t("2023-06-30"), string_t("sigmoid"), string_t("false")] &
+    ))
   end function
 
   function and_gate_with_skewed_training_data() result(test_passes)
@@ -212,7 +212,7 @@ function and_gate_with_skewed_training_data() result(test_passes)
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(tensor_t), allocatable, dimension(:,:) :: training_inputs, training_outputs
     type(tensor_t), allocatable, dimension(:) :: tmp, tmp2, test_inputs, expected_test_outputs, actual_outputs
-    type(trainable_engine_t) trainable_engine
+    type(trainable_network_t) trainable_network
     real, parameter :: tolerance = 1.E-02
     real, allocatable :: harvest(:,:,:)
     integer, parameter :: num_inputs=2, mini_batch_size = 1, num_iterations=20000
@@ -231,13 +231,13 @@ function and_gate_with_skewed_training_data() result(test_passes)
     training_outputs = reshape(tmp2, [mini_batch_size, num_iterations])
 
     mini_batches = [(mini_batch_t(input_output_pair_t(training_inputs(:,iter), training_outputs(:,iter))), iter=1, num_iterations)]        
-    trainable_engine = two_zeroed_hidden_layers()
+    trainable_network = two_zeroed_hidden_layers()
 
-    call trainable_engine%train(mini_batches, adam=.false., learning_rate=1.5)
+    call trainable_network%train(mini_batches, adam=.false., learning_rate=1.5)
 
     test_inputs = [tensor_t([true,true]), tensor_t([false,true]), tensor_t([true,false]), tensor_t([false,false])]
     expected_test_outputs = [(and(test_inputs(i)), i=1, size(test_inputs))]
-    actual_outputs = trainable_engine%infer(test_inputs)
+    actual_outputs = trainable_network%infer(test_inputs)
     test_passes = [(abs(actual_outputs(i)%values() - expected_test_outputs(i)%values()) < tolerance, i=1, size(actual_outputs))]
 
   contains
@@ -255,7 +255,7 @@ function not_and_gate_with_skewed_training_data() result(test_passes)
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(tensor_t), allocatable :: training_inputs(:,:), tmp(:), test_inputs(:)
     type(tensor_t), allocatable :: training_outputs(:,:), expected_test_outputs(:), tmp2(:)
-    type(trainable_engine_t) trainable_engine
+    type(trainable_network_t) trainable_network
     type(tensor_t), allocatable :: actual_outputs(:)
     real, parameter :: tolerance = 1.E-02
     real, allocatable :: harvest(:,:,:)
@@ -275,13 +275,13 @@ function not_and_gate_with_skewed_training_data() result(test_passes)
     training_outputs = reshape(tmp2, [mini_batch_size, num_iterations])
 
     mini_batches = [(mini_batch_t(input_output_pair_t(training_inputs(:,iter), training_outputs(:,iter))), iter=1, num_iterations)]        
-    trainable_engine = two_zeroed_hidden_layers()
+    trainable_network = two_zeroed_hidden_layers()
 
-    call trainable_engine%train(mini_batches, adam=.false., learning_rate=1.5)
+    call trainable_network%train(mini_batches, adam=.false., learning_rate=1.5)
 
     test_inputs = [tensor_t([true,true]), tensor_t([false,true]), tensor_t([true,false]), tensor_t([false,false])]
     expected_test_outputs = [(not_and(test_inputs(i)), i=1, size(test_inputs))]
-    actual_outputs = trainable_engine%infer(test_inputs)
+    actual_outputs = trainable_network%infer(test_inputs)
     test_passes = [(abs(actual_outputs(i)%values() - expected_test_outputs(i)%values()) < tolerance, i=1, size(actual_outputs))]
 
   contains
@@ -299,7 +299,7 @@ function or_gate_with_random_weights() result(test_passes)
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(tensor_t), allocatable :: training_inputs(:,:), test_inputs(:), actual_outputs(:)
     type(tensor_t), allocatable :: training_outputs(:,:), expected_test_outputs(:)
-    type(trainable_engine_t) trainable_engine
+    type(trainable_network_t) trainable_network
     real, parameter :: tolerance = 1.E-02
     real, allocatable :: harvest(:,:,:)
     integer, parameter :: num_inputs=2, mini_batch_size = 1, num_iterations=50000
@@ -320,13 +320,13 @@ function or_gate_with_random_weights() result(test_passes)
       mini_batches(iter) = mini_batch_t(input_output_pair_t(training_inputs(:,iter), training_outputs(:,iter)))
     end do
 
-    trainable_engine = two_random_hidden_layers()
+    trainable_network = two_random_hidden_layers()
 
-    call trainable_engine%train(mini_batches, adam=.false., learning_rate=1.5)
+    call trainable_network%train(mini_batches, adam=.false., learning_rate=1.5)
 
     test_inputs = [tensor_t([true,true]), tensor_t([false,true]), tensor_t([true,false]), tensor_t([false,false])]
     expected_test_outputs = [(or(test_inputs(i)), i=1, size(test_inputs))]
-    actual_outputs = trainable_engine%infer(test_inputs)
+    actual_outputs = trainable_network%infer(test_inputs)
     test_passes = [(abs(actual_outputs(i)%values() - expected_test_outputs(i)%values()) < tolerance, i=1, size(actual_outputs))]
 
   contains
@@ -344,7 +344,7 @@ function xor_gate_with_random_weights() result(test_passes)
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(tensor_t), allocatable, dimension(:,:) :: training_inputs, training_outputs 
     type(tensor_t), allocatable, dimension(:) :: actual_outputs, test_inputs, expected_test_outputs
-    type(trainable_engine_t) trainable_engine
+    type(trainable_network_t) trainable_network
     real, parameter :: tolerance = 1.E-02
     real, allocatable :: harvest(:,:,:)
 #ifdef __flang__
@@ -375,16 +375,16 @@ function xor_gate_with_random_weights() result(test_passes)
       mini_batches(iter) = mini_batch_t(input_output_pair_t(training_inputs(:,iter), training_outputs(:,iter)))
     end do
 
-    trainable_engine = two_random_hidden_layers()
+    trainable_network = two_random_hidden_layers()
 
-    call trainable_engine%train(mini_batches, adam=.true., learning_rate=1.5)
+    call trainable_network%train(mini_batches, adam=.true., learning_rate=1.5)
 
     test_inputs = [tensor_t([true,true]), tensor_t([false,true]), tensor_t([true,false]), tensor_t([false,false])]
     block
       integer i
 
       expected_test_outputs = [(local_xor(test_inputs(i)), i=1, size(test_inputs))]
-      actual_outputs = trainable_engine%infer(test_inputs)
+      actual_outputs = trainable_network%infer(test_inputs)
       test_passes = [(abs(actual_outputs(i)%values() - expected_test_outputs(i)%values()) < tolerance, i=1, size(actual_outputs))]
     end block
 
@@ -400,8 +400,8 @@ pure function local_xor(inputs) result(expected_outputs)
 
   end function
 
-  function perturbed_identity_network(perturbation_magnitude) result(trainable_engine)
-    type(trainable_engine_t) trainable_engine
+  function perturbed_identity_network(perturbation_magnitude) result(trainable_network)
+    type(trainable_network_t) trainable_network
     real, intent(in) :: perturbation_magnitude
     integer, parameter :: nodes_per_layer(*) = [2, 2, 2, 2]
     integer, parameter :: max_n = maxval(nodes_per_layer), layers = size(nodes_per_layer)
@@ -420,12 +420,12 @@ function perturbed_identity_network(perturbation_magnitude) result(trainable_eng
     call random_number(harvest)
     harvest = perturbation_magnitude*harvest
 
-    trainable_engine = trainable_engine_t( &
+    trainable_network = trainable_network_t( inference_engine_t( &
       nodes = nodes_per_layer, &
       weights = identity + harvest , & 
       biases = reshape([real:: [0,0], [0,0], [0,0]], [max_n, layers-1]), &
       metadata = [string_t("Identity"), string_t("Damian Rouson"), string_t("2023-09-18"), string_t("relu"), string_t("false")] &
-    )
+    ))
   end function
 
   function preserves_identity_mapping() result(test_passes)
@@ -433,17 +433,17 @@ function preserves_identity_mapping() result(test_passes)
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(input_output_pair_t), allocatable :: input_output_pairs(:)
     type(tensor_t), allocatable :: inputs(:)
-    type(trainable_engine_t)  trainable_engine
+    type(trainable_network_t)  trainable_network
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:)
     integer, parameter :: num_pairs = 100, num_epochs = 100, n_bins = 3
     integer i, bin, epoch
 
-    trainable_engine = perturbed_identity_network(perturbation_magnitude=0.)
+    trainable_network = perturbed_identity_network(perturbation_magnitude=0.)
 
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
+    associate(num_inputs => trainable_network%num_inputs(), num_outputs => trainable_network%num_outputs())
 
-      call assert(num_inputs == num_outputs,"trainable_engine_test_m(identity_mapping): # inputs == # outputs", &
+      call assert(num_inputs == num_outputs,"trainable_network_test_m(identity_mapping): # inputs == # outputs", &
         intrinsic_array_t([num_inputs, num_outputs]) &
       )
 #ifdef _CRAYFTN
@@ -461,16 +461,16 @@ function preserves_identity_mapping() result(test_passes)
 
       do epoch = 1,num_epochs
         mini_batches = [(mini_batch_t(input_output_pairs(bins(bin)%first():bins(bin)%last())), bin = 1, size(bins))]
-        call trainable_engine%train(mini_batches, cost, adam=.false., learning_rate=1.5)
+        call trainable_network%train(mini_batches, cost, adam=.false., learning_rate=1.5)
       end do
 
       block
         real, parameter :: tolerance = 1.E-06
 #if defined _CRAYFTN || __GFORTRAN__
         type(tensor_t), allocatable :: network_outputs(:)
-        network_outputs = trainable_engine%infer(inputs)
+        network_outputs = trainable_network%infer(inputs)
 #else
-        associate(network_outputs => trainable_engine%infer(inputs))
+        associate(network_outputs => trainable_network%infer(inputs))
 #endif
           test_passes = maxval(abs([(network_outputs(i)%values() - inputs(i)%values(), i=1,num_pairs)])) < tolerance
 #if defined _CRAYFTN || __GFORTRAN__
@@ -492,7 +492,7 @@ function perturbed_identity_converges() result(test_passes)
     type(mini_batch_t), allocatable :: mini_batches(:)
     type(input_output_pair_t), allocatable :: input_output_pairs(:)
     type(tensor_t), allocatable :: inputs(:)
-    type(trainable_engine_t)  trainable_engine
+    type(trainable_network_t)  trainable_network
     type(bin_t), allocatable :: bins(:)
     real, allocatable :: cost(:)
     integer, parameter :: num_pairs = 6
@@ -500,11 +500,11 @@ function perturbed_identity_converges() result(test_passes)
     integer, parameter :: num_bins = 5 
     integer i, bin, epoch
 
-    trainable_engine = perturbed_identity_network(perturbation_magnitude=0.1)
+    trainable_network = perturbed_identity_network(perturbation_magnitude=0.1)
 
-    associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs())
+    associate(num_inputs => trainable_network%num_inputs(), num_outputs => trainable_network%num_outputs())
 
-      call assert(num_inputs == num_outputs,"trainable_engine_test_m(identity_mapping): # inputs == # outputs", &
+      call assert(num_inputs == num_outputs,"trainable_network_test_m(identity_mapping): # inputs == # outputs", &
         intrinsic_array_t([num_inputs, num_outputs]) &
       )
 #ifdef _CRAYFTN
@@ -523,16 +523,16 @@ function perturbed_identity_converges() result(test_passes)
       do epoch = 1,num_epochs
         call shuffle(input_output_pairs)
         mini_batches = [(mini_batch_t(input_output_pairs(bins(bin)%first():bins(bin)%last())), bin = 1, size(bins))]
-        call trainable_engine%train(mini_batches, cost, adam=.true., learning_rate=1.5)
+        call trainable_network%train(mini_batches, cost, adam=.true., learning_rate=1.5)
       end do
 
       block
         real, parameter :: tolerance = 1.E-06
 #if defined _CRAYFTN || __GFORTRAN__
         type(tensor_t), allocatable :: network_outputs(:)
-        network_outputs = trainable_engine%infer(inputs)
+        network_outputs = trainable_network%infer(inputs)
 #else
-        associate(network_outputs => trainable_engine%infer(inputs))
+        associate(network_outputs => trainable_network%infer(inputs))
 #endif
           test_passes = maxval(abs([(network_outputs(i)%values() - inputs(i)%values(), i=1,num_pairs)])) < tolerance
 #if defined _CRAYFTN || __GFORTRAN__
@@ -545,4 +545,4 @@ function perturbed_identity_converges() result(test_passes)
 
   end function
 
-end module trainable_engine_test_m
+end module trainable_network_test_m

From e89c7b4c6c7e6d278f4b17a52dd306bad84c2b8d Mon Sep 17 00:00:00 2001
From: Damian Rouson <rouson@lbl.gov>
Date: Mon, 14 Oct 2024 01:12:32 -0700
Subject: [PATCH 6/6] fix(train-cloud-micro):let write_lines() open file

---
 demo/app/train-cloud-microphysics.F90 | 12 +++---------
 1 file changed, 3 insertions(+), 9 deletions(-)

diff --git a/demo/app/train-cloud-microphysics.F90 b/demo/app/train-cloud-microphysics.F90
index e96038f5d..9a8c5efcc 100644
--- a/demo/app/train-cloud-microphysics.F90
+++ b/demo/app/train-cloud-microphysics.F90
@@ -428,15 +428,9 @@ subroutine read_train_write(training_configuration, args, plot_file)
                       print *, epoch, average_cost
                       write(plot_file%plot_unit,*) epoch, average_cost
 
-                      block 
-                        integer net_unit
-
-                        open(newunit=net_unit, file=network_file, form='formatted', status='unknown', iostat=io_status, action='write')
-                        associate(json_file => trainable_network%to_json())
-                          call json_file%write_lines(string_t(network_file))
-                        end associate
-                        close(net_unit)
-                      end block
+                      associate(json_file => trainable_network%to_json())
+                        call json_file%write_lines(string_t(network_file))
+                      end associate
 
                     end if image_1_maybe_writes
 
