Improve access operator (#24)

* Separate `Matrix::get` in const and non-const

* Change `Matrix::get` into `operator()`

* Update definitions

* Update tests

src/cuda/DataFormats/Matrix.h

@@ -65,10 +65,12 @@ class Matrix {
   friend __host__ Matrix<U> operator*(const Matrix<U>& m1, const Matrix<E>& m2);
 
   template <typename E>
-  friend __host__ std::vector<E> operator*(const Matrix<E>& matrix, const std::vector<E>& vec);
+  friend __host__ std::vector<E> operator*(const Matrix<E>& matrix,
+                                           const std::vector<E>& vec);
 
   template <typename U, std::convertible_to<U> E>
-  friend __host__ std::vector<U> operator*(const Matrix<U>& matrix, const std::vector<E>& vec);
+  friend __host__ std::vector<U> operator*(const Matrix<U>& matrix,
+                                           const std::vector<E>& vec);
 
   __host__ Matrix<T>& operator+=(const Matrix<T>& other);
   template <typename E>
@@ -90,7 +92,10 @@ class Matrix {
 
 template <typename T>
 Matrix<T>::Matrix(int n_rows, int n_cols)
-    : m_nrows{n_rows}, m_ncols{n_cols}, m_size{n_rows * n_cols}, m_data(n_rows * n_cols) {}
+    : m_nrows{n_rows},
+      m_ncols{n_cols},
+      m_size{n_rows * n_cols},
+      m_data(n_rows * n_cols) {}
 
 template <typename T>
 template <typename E>
@@ -270,7 +275,8 @@ Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {
       throw(0);
     }
   } catch (int num) {
-    std::cout << "The two matrices can't be multiplied because their dimensions are not compatible. \n";
+    std::cout << "The two matrices can't be multiplied because their dimensions are not "
+                 "compatible. \n";
   }
 
   Matrix<T> result(N, M);
@@ -302,7 +308,8 @@ Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {
   cudaMemcpy(m_b.data, m2.data().data(), size_b, cudaMemcpyHostToDevice);
   const size_t shared_size{2 * block_size * block_size * sizeof(T)};
   matrix_multiply<<<grid, block, shared_size>>>(m_a, m_b, m_c, block_size);
-  cudaMemcpy(const_cast<T*>(result.data().data()), m_c.data, size_c, cudaMemcpyDeviceToHost);
+  cudaMemcpy(
+      const_cast<T*>(result.data().data()), m_c.data, size_c, cudaMemcpyDeviceToHost);
 
   cudaFree(d_a);
   cudaFree(d_b);
@@ -322,7 +329,8 @@ Matrix<T> operator*(const Matrix<T>& m1, const Matrix<E>& m2) {
       throw(0);
     }
   } catch (int num) {
-    std::cout << "The two matrices can't be multiplied because their dimensions are not compatible. \n";
+    std::cout << "The two matrices can't be multiplied because their dimensions are not "
+                 "compatible. \n";
   }
 
   Matrix<T> result(N, M);
@@ -355,7 +363,8 @@ Matrix<T> operator*(const Matrix<T>& m1, const Matrix<E>& m2) {
   cudaMemcpy(m_b.data, m2.data().data(), size_b, cudaMemcpyHostToDevice);
   const size_t shared_size{2 * block_size * block_size * sizeof(T)};
   matrix_multiply<<<grid, block, shared_size>>>(m_a, m_b, m_c, block_size);
-  cudaMemcpy(const_cast<T*>(result.data().data()), m_c.data, size_c, cudaMemcpyDeviceToHost);
+  cudaMemcpy(
+      const_cast<T*>(result.data().data()), m_c.data, size_c, cudaMemcpyDeviceToHost);
 
   cudaFree(d_a);
   cudaFree(d_b);

src/cuda/DataFormats/VectorKernels.h

@@ -73,7 +73,10 @@ __global__ void vec_divide(T* a, E constant, int n) {
 }
 
 template <typename T, typename E, typename U>
-__global__ void matrix_multiply(const matrix_t<T> a, const matrix_t<E> b, matrix_t<U> c, int block_size) {
+__global__ void matrix_multiply(const matrix_t<T> a,
+                                const matrix_t<E> b,
+                                matrix_t<U> c,
+                                int block_size) {
   // Allocate memory on shared memory
   extern __shared__ int shared[];
   int* s_a{shared};

src/cuda/include/ErrorFunction.h

@@ -14,7 +14,8 @@ using shared = std::shared_ptr<T>;
 
 template <typename T, typename W, template <typename E> typename Activator>
 struct MeanSquaredError {
-  __host__ double operator()(const std::vector<T>& node_values, const std::vector<T>& expected_values) {
+  __host__ double operator()(const std::vector<T>& node_values,
+                             const std::vector<T>& expected_values) {
     double error{};
     int N{node_values.size()};
     for (int node_index{}; node_index < N; ++node_index) {
@@ -35,7 +36,8 @@ struct MeanSquaredError {
       int N{layers[layer_id]->size()};
       std::vector<double> delta(N);
       for (int node_index{}; node_index < N; ++node_index) {
-        delta[node_index] = (*layers[layer_id])[node_index] - static_cast<T>(expected_values[node_index]);
+        delta[node_index] =
+            (*layers[layer_id])[node_index] - static_cast<T>(expected_values[node_index]);
       }
 
       return delta;
@@ -45,7 +47,8 @@ struct MeanSquaredError {
       std::vector<double> delta(N);
 
       for (int node_index{}; node_index < N; ++node_index) {
-        std::vector<double> previous_delta{grad(expected_values, layer_id + 1, layers, weights)};
+        std::vector<double> previous_delta{
+            grad(expected_values, layer_id + 1, layers, weights)};
         delta[node_index] = act.grad((*layers[layer_id])[node_index]) *
                             (weights[layer_id]->transpose() * previous_delta)[node_index];
       }

src/cuda/include/Layer.h

@@ -47,7 +47,8 @@ template <typename T>
 Layer<T>::Layer(int n_nodes) : m_nodes(n_nodes), n_nodes{n_nodes} {}
 
 template <typename T>
-Layer<T>::Layer(std::vector<T> nodes) : m_nodes{std::move(nodes)}, n_nodes{m_nodes.size()} {}
+Layer<T>::Layer(std::vector<T> nodes)
+    : m_nodes{std::move(nodes)}, n_nodes{m_nodes.size()} {}
 
 template <typename T>
 Layer<T>::Layer(std::stringstream& stream) {
@@ -77,7 +78,8 @@ void Layer<T>::load(std::stringstream& stream) {
       }
       ++node_index;
     } catch (int num) {
-      std::cout << "The data provided exceedes the number of nodes expected for the layer\n";
+      std::cout
+          << "The data provided exceedes the number of nodes expected for the layer\n";
     }
   }
 }
@@ -90,7 +92,8 @@ void Layer<T>::set_node_data(int i, T value) {
     }
     m_nodes[i] = value;
   } catch (...) {
-    std::cout << "The index " << i << " is larger than the number of nodes in the layer\n";
+    std::cout << "The index " << i
+              << " is larger than the number of nodes in the layer\n";
   }
 }
 

src/cuda/include/Network.h

@@ -54,7 +54,9 @@ class Network {
   void set_bias_data(int layer_id, std::vector<W> bias_vector);
   void set_bias_data(int layer_id, shared<std::vector<W>> bias_vector_ptr);
 
-  std::vector<T> forward_propatation(shared<Layer<T>>, shared<Matrix<W>>, shared<std::vector<W>>);
+  std::vector<T> forward_propatation(shared<Layer<T>>,
+                                     shared<Matrix<W>>,
+                                     shared<std::vector<W>>);
   void forward_propatation();
 
   template <typename U>
@@ -94,7 +96,8 @@ Network<T, W, Activator, Loss>::Network(const std::vector<int>& nodes_per_layer)
       m_bias(n_layers - 1) {
   for (int i{}; i < n_layers - 1; ++i) {
     m_layers[i] = std::make_shared<Layer<T>>(nodes_per_layer[i]);
-    m_weights[i] = std::make_shared<Matrix<W>>(nodes_per_layer[i + 1], nodes_per_layer[i]);
+    m_weights[i] =
+        std::make_shared<Matrix<W>>(nodes_per_layer[i + 1], nodes_per_layer[i]);
     m_bias[i] = std::make_shared<std::vector<W>>(nodes_per_layer[i + 1]);
 
     // Generate random weight matrices
@@ -160,7 +163,8 @@ template <typename T,
           typename Activator,
           template <typename E, typename LW, template <typename K> typename Act>
           typename Loss>
-void Network<T, W, Activator, Loss>::set_matrix_data(int layer_id, Matrix<W> weight_matrix) {
+void Network<T, W, Activator, Loss>::set_matrix_data(int layer_id,
+                                                     Matrix<W> weight_matrix) {
   m_weights[layer_id] = std::make_shared<Matrix<W>>(weight_matrix);
 }
 
@@ -170,7 +174,8 @@ template <typename T,
           typename Activator,
           template <typename E, typename LW, template <typename K> typename Act>
           typename Loss>
-void Network<T, W, Activator, Loss>::set_matrix_data(int layer_id, shared<Matrix<W>> weight_matrix_ptr) {
+void Network<T, W, Activator, Loss>::set_matrix_data(
+    int layer_id, shared<Matrix<W>> weight_matrix_ptr) {
   m_weights[layer_id] = weight_matrix_ptr;
 }
 
@@ -180,7 +185,8 @@ template <typename T,
           typename Activator,
           template <typename E, typename LW, template <typename K> typename Act>
           typename Loss>
-void Network<T, W, Activator, Loss>::set_bias_data(int layer_id, std::vector<W> bias_vector) {
+void Network<T, W, Activator, Loss>::set_bias_data(int layer_id,
+                                                   std::vector<W> bias_vector) {
   m_bias[layer_id] = std::make_shared<std::vector<W>>(bias_vector);
 }
 
@@ -190,7 +196,8 @@ template <typename T,
           typename Activator,
           template <typename E, typename LW, template <typename K> typename Act>
           typename Loss>
-void Network<T, W, Activator, Loss>::set_bias_data(int layer_id, shared<std::vector<W>> bias_vector_ptr) {
+void Network<T, W, Activator, Loss>::set_bias_data(
+    int layer_id, shared<std::vector<W>> bias_vector_ptr) {
   m_bias[layer_id] = bias_vector_ptr;
 }
 
@@ -200,9 +207,10 @@ template <typename T,
           typename Activator,
           template <typename E, typename LW, template <typename K> typename Act>
           typename Loss>
-std::vector<T> Network<T, W, Activator, Loss>::forward_propatation(shared<Layer<T>> layer,
-                                                                   shared<Matrix<W>> weight_matrix,
-                                                                   shared<std::vector<W>> bias_vector) {
+std::vector<T> Network<T, W, Activator, Loss>::forward_propatation(
+    shared<Layer<T>> layer,
+    shared<Matrix<W>> weight_matrix,
+    shared<std::vector<W>> bias_vector) {
   std::vector<W> next_layer_nodes{*weight_matrix * layer->nodes() + *bias_vector};
 
   return Activator<T>()(next_layer_nodes);
@@ -216,7 +224,8 @@ template <typename T,
           typename Loss>
 void Network<T, W, Activator, Loss>::forward_propatation() {
   for (int i{}; i < n_layers - 1; ++i) {
-    std::vector<T> new_layer_data{forward_propatation(m_layers[i], m_weights[i], m_bias[i])};
+    std::vector<T> new_layer_data{
+        forward_propatation(m_layers[i], m_weights[i], m_bias[i])};
     m_layers[i + 1]->set_node_data(new_layer_data);
   }
 }
@@ -245,7 +254,8 @@ template <typename T,
           template <typename E, typename LW, template <typename K> typename Act>
           typename Loss>
 template <typename U>
-void Network<T, W, Activator, Loss>::back_propagation(double eta, const std::vector<U>& target) {
+void Network<T, W, Activator, Loss>::back_propagation(double eta,
+                                                      const std::vector<U>& target) {
   for (int layer_id{n_layers - 2}; layer_id >= 0; --layer_id) {
     back_propagation(target, layer_id, eta);
   }
@@ -304,7 +314,8 @@ void Network<T, W, Activator, Loss>::import_network(const std::string& filepath)
       bias.push_back(std::stod(value));
     }
 
-    m_weights[i] = std::make_shared<Matrix<W>>(m_weights[i]->nrows(), m_weights[i]->ncols(), weights);
+    m_weights[i] = std::make_shared<Matrix<W>>(
+        m_weights[i]->nrows(), m_weights[i]->ncols(), weights);
     m_bias[i] = std::make_shared<std::vector<W>>(bias);
   }
 }

src/nnhep/DataFormats/Matrix.hpp

@@ -88,7 +88,16 @@ namespace nnhep {
     /// @param i The row index of the element to get
     /// @param j The column index of the element to get
     /// @return The element of the matrix
-    inline const T get(int i, int j) const;
+    ///
+    /// @note returns by reference
+    inline T& operator()(int i, int j);
+    /// @brief Get an element of the matrix
+    /// @param i The row index of the element to get
+    /// @param j The column index of the element to get
+    /// @return The element of the matrix
+    ///
+    /// @note returns by const reference
+    inline const T& operator()(int i, int j) const;
 
     ///
     inline Matrix transpose();
@@ -193,7 +202,10 @@ namespace nnhep {
 
   template <typename T>
   Matrix<T>::Matrix(int n_rows, int n_cols)
-      : m_nrows{n_rows}, m_ncols{n_cols}, m_size{n_rows * n_cols}, m_data(n_rows * n_cols) {}
+      : m_nrows{n_rows},
+        m_ncols{n_cols},
+        m_size{n_rows * n_cols},
+        m_data(n_rows * n_cols) {}
 
   template <typename T>
   template <typename E>
@@ -276,7 +288,12 @@ namespace nnhep {
   }
 
   template <typename T>
-  const T Matrix<T>::get(int i, int j) const {
+  T& Matrix<T>::operator()(int i, int j) {
+    return m_data[j + m_ncols * i];
+  }
+
+  template <typename T>
+  const T& Matrix<T>::operator()(int i, int j) const {
     return m_data[j + m_ncols * i];
   }
 
@@ -286,7 +303,7 @@ namespace nnhep {
 
     for (int i{}; i < this->m_nrows; ++i) {
       for (int j{}; j < this->m_ncols; ++j) {
-        matrix.set_data(j, i, this->get(i, j));
+        matrix.set_data(j, i, (*this)(i, j));
       }
     }
 
@@ -335,14 +352,15 @@ namespace nnhep {
         throw(0);
       }
     } catch (int num) {
-      std::cout << "The two matrices can't be multiplied because their dimensions are not compatible. \n";
+      std::cout << "The two matrices can't be multiplied because their dimensions are "
+                   "not compatible. \n";
     }
 
     for (int i{}; i < m1.m_nrows; ++i) {
       for (int j{}; j < m2.m_ncols; ++j) {
         T sum{};
         for (int k{}; k < m1.m_ncols; ++k) {
-          sum += m1.get(i, k) * m2.get(k, j);
+          sum += m1(i, k) * m2(k, j);
         }
         result.set_data(i, j, sum);
       }
@@ -360,14 +378,15 @@ namespace nnhep {
         throw(0);
       }
     } catch (int num) {
-      std::cout << "The two matrices can't be multiplied because their dimensions are not compatible. \n";
+      std::cout << "The two matrices can't be multiplied because their dimensions are "
+                   "not compatible. \n";
     }
 
     for (int i{}; i < m1.m_nrows; ++i) {
       for (int j{}; j < m2.m_ncols; ++j) {
         T sum{};
         for (int k{}; k < m1.m_ncols; ++k) {
-          sum += m1.get(i, k) * m2.get(k, j);
+          sum += m1(i, k) * m2(k, j);
         }
         result.set_data(i, j, sum);
       }
@@ -447,21 +466,27 @@ namespace nnhep {
   template <typename T>
   template <typename E>
   Matrix<T>& Matrix<T>::operator*=(E constant) {
-    std::transform(this->m_data.cbegin(), this->m_data.cend(), this->m_data.begin(), [constant](auto x) {
-      x *= constant;
-      return x;
-    });
+    std::transform(this->m_data.cbegin(),
+                   this->m_data.cend(),
+                   this->m_data.begin(),
+                   [constant](auto x) {
+                     x *= constant;
+                     return x;
+                   });
 
     return *this;
   }
 
   template <typename T>
   template <typename E>
   Matrix<T>& Matrix<T>::operator/=(E constant) {
-    std::transform(this->m_data.cbegin(), this->m_data.cend(), this->m_data.begin(), [constant](auto x) {
-      x /= constant;
-      return x;
-    });
+    std::transform(this->m_data.cbegin(),
+                   this->m_data.cend(),
+                   this->m_data.begin(),
+                   [constant](auto x) {
+                     x /= constant;
+                     return x;
+                   });
 
     return *this;
   }