Neural Network

Started on implementing neural network from NNFS. I've done ReLU and stopped at p.104. Softmax is not ready.

include/modules/neural_networks/activation_functions/ReLU.h

@@ -0,0 +1,28 @@
+#pragma once
+
+#include "./core/omp_config.h"
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+#include "./utils/random.h"
+
+
+namespace neural_networks{
+
+	template <typename T>
+	struct activation_ReLU{
+
+			utils::Matrix<T> outputs;
+
+	        void forward(utils::Matrix<T> inputs){
+			    outputs = numerics::max(inputs, T{0});
+			    //outputs.print();
+	        }
+
+
+
+	};
+
+
+
+} // end namespace neural_networks

include/modules/neural_networks/activation_functions/Softmax.h

@@ -0,0 +1,28 @@
+#pragma once
+
+#include "./core/omp_config.h"
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+#include "./utils/random.h"
+
+
+namespace neural_networks{
+
+	template <typename T>
+	struct activation_softmax{
+
+			utils::Matrix<T> outputs;
+
+	        void forward(utils::Matrix<T> inputs){
+			    //outputs = numerics::max(inputs, T{0});
+			    //outputs.print();
+	        }
+
+
+
+	};
+
+
+
+} // end namespace neural_networks

include/modules/neural_networks/datasets/spiral.h

@@ -0,0 +1,61 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "./utils/matrix.h"
+#include "./utils/vector.h"
+#include "./utils/random.h"
+
+//#include <math.h>
+
+
+namespace neural_networks{
+
+		template <typename T>
+		void create_spital_data(const uint64_t samples, const uint64_t classes, utils::Matrix<T>& X, utils::Vector<T>& y) {
+
+			const uint64_t rows = samples*classes;
+			T r, t;
+			uint64_t row_idx;
+
+
+			if ((rows != X.rows()) || (X.cols() != 2)){
+				X.resize(samples*classes, 2);
+			}
+			if (rows != y.size()){
+				y.resize(rows);
+			}
+
+			for (uint64_t i = 0; i < classes; ++i){
+				for (uint64_t j = 0; j < samples; ++j){
+					r = static_cast<T>(j)/static_cast<T>(samples);
+					t =  static_cast<T>(i)*4.0 + (4.0+r);
+					row_idx = (i*samples) + j;
+
+					X(row_idx, 0) = r*std::cos(t*2.5) + utils::random(T{-0.15}, T{0.15});
+					X(row_idx, 1) = r*std::sin(t*2.5) + utils::random(T{-0.15}, T{0.15});
+					y[row_idx] = i;
+				}
+			}
+
+			
+
+			/*
+
+			utils::Matrix<T> X(static_cast<uint64_t>(samples*classes), 3, T{0});
+		    
+		    const uint64_t rows = A.rows();
+		    const uint64_t cols = A.cols();
+
+		    if (rows != x.size()) {
+		        throw std::runtime_error("inplace_matadd_colvec: dimension mismatch");
+		    }
+
+		    for (uint64_t i = 0; i < cols; ++i) {
+		        for (uint64_t j = 0; j < rows; ++j) {
+		            A(j, i) +=  x[j];
+		        }
+		    }*/
+		}
+
+
+} // end namesoace NN

include/modules/neural_networks/layers/dense_layer.h

@@ -0,0 +1,42 @@
+#pragma once
+
+#include "./core/omp_config.h"
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+#include "./utils/random.h"
+
+
+namespace neural_networks{
+
+	template <typename T>
+	struct dense_layer{
+			utils::Matrix<T> weights;
+			utils::Vector<T> biases;
+			utils::Matrix<T> outputs;
+			
+			// Default Constructor
+			dense_layer() = default;
+
+			// Constructor
+	    	dense_layer(const uint64_t n_inputs, const uint64_t n_neurons){
+	        	
+	        	weights.random(n_inputs, n_neurons, -1, 1);
+	        	biases.resize(n_neurons, T{0});
+	        	weights.print();
+	        	//outputs.resize()
+	        	
+	        }
+
+
+	        void forward(utils::Matrix<T> inputs){
+			    outputs = numerics::matadd(numerics::matmul_auto(inputs, (weights)), biases, "row");
+	        }
+
+
+
+	};
+
+
+
+} // end namespace neural_networks

include/modules/neural_networks/neural_networks.h

@@ -0,0 +1,9 @@
+// #include "./modules/neural_networks/neural_networks.h"
+#pragma once
+
+#include "datasets/spiral.h"
+
+#include "layers/dense_layer.h"
+
+
+#include "activation_functions/ReLU.h"

include/numerics/matadd.h

@@ -0,0 +1,229 @@
+#ifndef _matadd_n_
+#define _matadd_n_
+
+
+#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+
+namespace numerics{
+	
+// =================================================
+//   y = A * x    (Matrix–Vector product)
+// =================================================
+	template <typename T>
+	void inplace_matadd_colvec(utils::Matrix<T>& A, const utils::Vector<T>& x) {
+	    
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+	    if (rows != x.size()) {
+	        throw std::runtime_error("inplace_matadd_colvec: dimension mismatch");
+	    }
+
+	    for (uint64_t i = 0; i < cols; ++i) {
+	        for (uint64_t j = 0; j < rows; ++j) {
+	            A(j, i) +=  x[j];
+	        }
+	    }
+	}
+
+	template <typename T>
+	void inplace_matadd_rowvec(utils::Matrix<T>& A, const utils::Vector<T>& x) {
+	    
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+	    if (cols != x.size()) {
+	        throw std::runtime_error("inplace_matadd_rowvec: dimension mismatch");
+	    }
+
+	    for (uint64_t i = 0; i < cols; ++i) {
+	        for (uint64_t j = 0; j < rows; ++j) {
+	            A(j, i) +=  x[i];
+	        }
+	    }
+	}
+
+	template <typename T>
+	utils::Matrix<T> matadd_colvec(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
+	    
+	    //const uint64_t rows = A.rows();
+	    //const uint64_t cols = A.cols();
+
+	    utils::Matrix<T> B = A;
+
+	    inplace_matadd_colvec(B, x);
+
+	    return B;
+	}
+
+	template <typename T>
+	utils::Matrix<T> matadd_rowvec(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
+	    
+	    //const uint64_t rows = A.rows();
+	    //const uint64_t cols = A.cols();
+
+	    utils::Matrix<T> B = A;
+
+	    inplace_matadd_rowvec(B, x);
+
+	    return B;
+	}
+
+	template <typename T>
+	utils::Matrix<T> matadd(const utils::Matrix<T>& A, const utils::Vector<T>& x, std::string method = "auto"){
+
+		const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+	    const uint64_t N = x.size();
+
+	    if (method=="auto"){
+	    	
+			if (rows==cols){
+				throw std::runtime_error("matadd: too many options for dimensions");
+			} else if (rows == N){
+				return matadd_rowvec(A, x);
+			} else if (cols == N){
+				return matadd_colvec(A, x);
+			}else{
+				throw std::runtime_error("matadd: undefined fault - auto");
+			}
+	    }else if(method=="row"){
+				return matadd_rowvec(A, x);
+		} else if (method=="col"){
+			return matadd_colvec(A, x);
+		}else{
+			throw std::runtime_error("matadd: undefined fault - defined method");
+		}
+    }
+
+
+
+
+	/*
+	// -------------- Collapse(2) OpenMP ----------------
+	template <typename T>
+	utils::Vector<T> matvec_omp(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
+	    if (A.cols() != x.size()) {
+	        throw std::runtime_error("matvec: dimension mismatch");
+	    }
+
+	    const uint64_t m = A.rows();
+	    const uint64_t n = A.cols();
+
+	    utils::Vector<T> y(m, T{0}); // <-- y has length m (rows)
+
+
+        const T* xptr = x.data();
+        const T* Aptr = A.data();	// row-major: A(i,j) == Aptr[i*n + j]
+
+        // Each row i is an independent dot product: y[i] = dot(A[i,*], x)
+	    #pragma omp parallel for schedule(static)
+	    for (uint64_t i = 0; i < m; ++i) {
+            const T* row = Aptr + i * n;     // contiguous row i
+            T acc = T{0};
+            #pragma omp simd reduction(+:acc)
+	        for (uint64_t j = 0; j < n; ++j) {
+	            acc += row[j] * xptr[j];
+	        }
+	        y[i] = acc;
+		}
+
+	    return y;
+	}
+
+	// -------------- Auto OpenMP ----------------
+	template <typename T>
+	utils::Vector<T> matvec_auto(const utils::Matrix<T>& A,
+	                             const utils::Vector<T>& x) {
+
+	    
+	    uint64_t work = A.rows() * A.cols();
+
+	    bool can_parallel = omp_config::omp_parallel_allowed();
+		#ifdef _OPENMP
+		    int threads = omp_get_max_threads();
+		#else
+		    int threads = 1;
+		#endif
+
+	    if (can_parallel || work > static_cast<uint64_t>(threads) * 4ull) {
+	        return matvec_omp(A,x);
+	    }
+	    else{
+	    	// Safe fallback
+	    	return matvec(A,x);
+	    }
+	    
+	}
+
+// =================================================
+//   y = x * A    (Vector–Matrix product)
+// =================================================
+	template <typename T>
+	utils::Vector<T> vecmat(const utils::Vector<T>& x, const utils::Matrix<T>& A) {
+	    if (x.size() != A.rows()) {
+	        throw std::runtime_error("vecmat: dimension mismatch");
+	    }
+	    const uint64_t m = A.rows();
+	    const uint64_t n = A.cols();
+
+	    utils::Vector<T> y(n, T{0});
+
+	    for (uint64_t j = 0; j < n; ++j) {
+	        for (uint64_t i = 0; i < m; ++i) {
+	            y[j] += x[i] * A(i, j);
+	        }
+	    }
+	    return y;
+	}
+
+	// -------------- Collapse(2) OpenMP ----------------
+	template <typename T>
+	utils::Vector<T> vecmat_omp(const utils::Vector<T>& x, const utils::Matrix<T>& A) {
+	    if (x.size() != A.rows()) {
+	        throw std::runtime_error("vecmat: dimension mismatch");
+	    }
+	    const uint64_t m = A.rows();
+	    const uint64_t n = A.cols();
+
+	    utils::Vector<T> y(n, T{0});
+	    #pragma omp parallel for schedule(static)
+	    for (uint64_t j = 0; j < n; ++j) {
+	        T acc = T{0};
+	        for (uint64_t i = 0; i < m; ++i) {
+	            acc += x[i] * A(i, j);
+	        }
+	        y[j] = acc;
+	    }
+	    return y;
+	}
+
+	// -------------- Auto OpenMP ----------------
+	template <typename T>
+	utils::Vector<T> vecmat_auto(const utils::Vector<T>& x,
+								 const utils::Matrix<T>& A) {
+	    
+	    uint64_t work = A.rows() * A.cols();
+
+	    bool can_parallel = omp_config::omp_parallel_allowed();
+		#ifdef _OPENMP
+		    int threads = omp_get_max_threads();
+		#else
+		    int threads = 1;
+		#endif
+
+	    if (can_parallel || work > static_cast<uint64_t>(threads) * 4ull) {
+	        return vecmat_omp(x,A);
+	    }
+	    else{
+	    	// Safe fallback
+	    	return vecmat(x,A);
+	    }
+	    
+	}
+*/
+
+} // namespace numerics
+
+#endif // _matadd_n_

include/numerics/matmul.h

@@ -98,6 +98,7 @@ utils::Matrix<T> matmul_auto(const utils::Matrix<T>& A,
 
     // Tiny problems: serial is cheapest.
     if (!can_parallel || work < threads*4ull) {
+        
         return matmul(A,B);
     }
     // Plenty of (i,j) work → collapse(2) is a great default.

include/numerics/max.h

@@ -17,6 +17,35 @@ namespace numerics{
 		}
 	}
 
+	template <typename T>
+	void inplace_max(utils::Matrix<T>& A, const T b){
+
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+
+				if (b > A(i,j)){
+					//std::cout << A(i,j) << std::endl;
+					A(i,j) = b;
+					//std::cout << A(i,j) << std::endl;
+				}
+			}
+		}
+	}
+
+	template <typename T>
+	utils::Matrix<T> max(const utils::Matrix<T>& A, const T b){
+
+		utils::Matrix<T> B = A;
+
+		inplace_max(B, b);
+		
+
+		return B;
+	}
+
 
 
 } // namespace numerics

include/numerics/numerics.h

@@ -7,6 +7,7 @@
 #include "./numerics/inverse.h"
 #include "./numerics/matmul.h"
 #include "./numerics/matvec.h"
+#include "./numerics/matadd.h"
 #include "./numerics/min.h"
 #include "./numerics/max.h"
 #include "./numerics/abs.h"

include/numerics/transpose.h

@@ -103,6 +103,7 @@ namespace numerics{
 	        uint64_t threads = 1;
 		#endif
 
+
         if (can_parallel && work > threads * 4ull) {
             inplace_transpose_square_omp(A);
         }else {
@@ -118,8 +119,10 @@ namespace numerics{
 
         uint64_t work = A.rows() * A.cols();
 
+
+
         if (rows==cols){
-        	utils::Matrix<T> B(rows, cols, T{0});
+        	utils::Matrix<T> B = A;
         	inplace_transpose_square_auto(B);
         	return B;
         }

include/utils/matrix.h

@@ -2,12 +2,16 @@
 #define _matrix_n_
 
 #include "./utils/vector.h"
+#include "./utils/random.h"
 
 #ifdef _OPENMP
 	#include <omp.h>
 #endif
 
 
+#include <initializer_list>
+
+
 #include <iomanip>
 
 namespace utils{
@@ -25,6 +29,72 @@ public:
 	    : rows_(rows), cols_(cols), data_(rows * cols, value) {}
 
 
+		// Construct from list-of-lists:
+		// utils::Mf A{{1,2,3}, {4,5,6}};
+		Matrix(std::initializer_list<std::initializer_list<T>> init) {
+		    rows_ = static_cast<uint64_t>(init.size());
+		    if (rows_ > 0) {
+			    cols_ = static_cast<uint64_t>(init.begin()->size());
+			} else {
+			    cols_ = 0;
+			}
+
+		    // Validate all rows have the same length
+		    for (const auto& row : init) {
+		        if (row.size() != cols_) {
+		            throw std::runtime_error("Matrix(list of lists): ragged rows");
+		        }
+		    }
+
+		    data_.resize(rows_ * cols_);
+			uint64_t r = 0;
+			for (uint64_t i = 0; i < init.size(); ++i, ++r){
+			    const std::initializer_list<T>& row = *(init.begin() + i);
+
+			    uint64_t c = 0;
+			    for (uint64_t j = 0; j < row.size(); ++j, ++c){
+			        const T& val = *(row.begin() + j);
+			        data_[r * cols_ + c] = val;
+			    }
+			}
+		}
+
+
+		// Assign from list-of-lists after default construction:
+		// utils::Md M; M = {{1,2},{3,4}};
+		Matrix& operator=(std::initializer_list<std::initializer_list<T>> init) {
+		    // Set sizes
+		    rows_ = static_cast<uint64_t>(init.size());
+		    if (rows_ > 0) {
+		        cols_ = static_cast<uint64_t>((init.begin())->size());
+		    } else {
+		        cols_ = 0;
+		    }
+
+		    // Validate: all rows must have same length
+		    for (uint64_t i = 0; i < rows_; ++i) {
+		        const std::initializer_list<T>& row = *(init.begin() + i);
+		        if (row.size() != cols_) {
+		            throw std::runtime_error("Matrix(list of lists): ragged rows");
+		        }
+		    }
+
+		    // Allocate storage
+		    data_.resize(rows_ * cols_);
+
+		    // Copy data row by row
+		    for (uint64_t i = 0; i < rows_; ++i) {
+		        const std::initializer_list<T>& row = *(init.begin() + i);
+		        for (uint64_t j = 0; j < cols_; ++j) {
+		            const T& val = *(row.begin() + j);
+		            data_[i * cols_ + j] = val;
+		        }
+		    }
+
+		    return *this;
+		}
+
+
 		//#            MATRIX: basic properties            #
 	    uint64_t rows() const noexcept {return rows_;}
 	    uint64_t cols() const noexcept {return cols_;}
@@ -45,6 +115,20 @@ public:
 
 		}
 
+		void random(uint64_t rows, uint64_t cols, const T& lower, const T& higher){
+			rows_ = rows;
+			cols_ = cols;
+			data_.resize(rows_*cols_, 0);
+
+		    // Copy data row by row
+		    for (uint64_t i = 0; i < rows_; ++i) {
+		        for (uint64_t j = 0; j < cols_; ++j) {
+		            data_[i * cols_ + j] = utils::random(lower, higher);
+		        }
+		    }
+
+		}
+
 
 		//#            MATRIX: row helpers (copy out)            #
 		// Read whole row as an owning Vector<T>

include/utils/random.h

@@ -0,0 +1,37 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include <random>
+
+namespace utils{
+
+    // Shared engine
+    inline std::mt19937& rng() {
+        static std::random_device rd;
+        static std::mt19937 gen(rd());
+        return gen;
+    }
+
+    // Integral overload
+    template <
+        typename T,
+        typename std::enable_if<std::is_integral<T>::value, int>::type = 0
+    >
+    T random(T low, T high) {
+        std::uniform_int_distribution<T> dist(low, high);
+        return dist(rng());
+    }
+
+    // Floating-point overload
+    template <
+        typename T,
+        typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0
+    >
+    T random(T low, T high) {
+        std::uniform_real_distribution<T> dist(low, high);
+        return dist(rng());
+    }
+
+
+} // end namespace utils
+

include/utils/utils.h

@@ -4,3 +4,4 @@
 #include "./utils/vector.h"
 #include "./utils/matrix.h"
 #include "./utils/generators.h"
+#include "./utils/random.h"

include/utils/vector.h

@@ -5,6 +5,8 @@
 #include <vector>
 #include <random>
 
+#include <initializer_list>  
+
 
 #include <cstdint>
 #include <type_traits>
@@ -31,6 +33,9 @@ public:
 	        v.resize(size, value);
 	    }
 
+		// Construct from a braced list: utils::Vf v{1,2,3};
+		Vector(std::initializer_list<T> init) : v(init) {}
+
 
 
 	    //##########################################################
@@ -47,6 +52,13 @@ public:
 		// a = vector[2];
 		const T& operator[](uint64_t idx) const { return v[idx]; }
 
+
+		// Assign from a braced list after default construction:
+		Vector& operator=(std::initializer_list<T> init) {
+		    v = init;
+		    return *this;
+		}
+
 		// vector.size();
 		uint64_t size() const noexcept { return v.size(); } 
 

makefile

@@ -51,7 +51,7 @@ OMP_MAX_LEVELS  ?= 1            # 1 = no nested teams; set 2+ to allow nesting
 OMP_THREADS     ?= 16           # e.g. "16" or "8,4" for nested (outer,inner)
 OMP_DYNAMIC     ?= TRUE         # TRUE/FALSE: let runtime adjust threads
 OMP_SCHEDULE    ?= STATIC     	# STATIC recommended for matvec/matmul
-OMP_DISPLAY_ENV ?= TRUE        # TRUE to print runtime config at startup
+OMP_DISPLAY_ENV ?= FALSE        # TRUE to print runtime config at startup
 
 # Export OMP defaults so child makes or tools see them (not strictly required)
 export OMP_PROC_BIND

obj/main.d

@@ -1,15 +1,16 @@
 obj/main.o: src/main.cpp include/./core/omp_config.h \
  include/./utils/utils.h include/./utils/vector.h \
- include/./utils/matrix.h include/./utils/generators.h \
- include/./utils/generators/linspace.h include/utils/vector.h \
- include/./numerics/numerics.h include/./numerics/initializers/eye.h \
- include/./numerics/matequal.h include/./numerics/abs.h \
- include/./numerics/transpose.h include/./numerics/inverse.h \
+ include/./utils/matrix.h include/./utils/random.h \
+ include/./utils/generators.h include/./utils/generators/linspace.h \
+ include/utils/vector.h include/./numerics/numerics.h \
+ include/./numerics/initializers/eye.h include/./numerics/matequal.h \
+ include/./numerics/abs.h include/./numerics/transpose.h \
+ include/./numerics/inverse.h \
  include/./numerics/inverse/inverse_gauss_jordan.h \
  include/./numerics/inverse/inverse_lu.h include/./decomp/lu.h \
  include/./numerics/matmul.h include/./numerics/matvec.h \
- include/./numerics/min.h include/./numerics/max.h \
- include/./numerics/interpolation1d.h \
+ include/./numerics/matadd.h include/./numerics/min.h \
+ include/./numerics/max.h include/./numerics/interpolation1d.h \
  include/./numerics/interpolation1d/interpolation1d_barycentric.h \
  include/./numerics/interpolation1d/interpolation1d_base.h \
  include/./numerics/interpolation1d/interpolation1d_cubic_spline.h \
@@ -19,11 +20,16 @@ obj/main.o: src/main.cpp include/./core/omp_config.h \
  include/./decomp/decomp.h include/./modules/mesh/mesh.h \
  include/modules/mesh/mesh1d.h include/modules/fluids/fluids.h \
  include/modules/fluids/diffusion1d.h include/core/global_config.h \
- include/utils/matrix.h
+ include/utils/matrix.h \
+ include/./modules/neural_networks/neural_networks.h \
+ include/./modules/neural_networks/datasets/spiral.h \
+ include/./modules/neural_networks/layers/dense_layer.h \
+ include/./modules/neural_networks/activation_functions/ReLU.h
 include/./core/omp_config.h:
 include/./utils/utils.h:
 include/./utils/vector.h:
 include/./utils/matrix.h:
+include/./utils/random.h:
 include/./utils/generators.h:
 include/./utils/generators/linspace.h:
 include/utils/vector.h:
@@ -38,6 +44,7 @@ include/./numerics/inverse/inverse_lu.h:
 include/./decomp/lu.h:
 include/./numerics/matmul.h:
 include/./numerics/matvec.h:
+include/./numerics/matadd.h:
 include/./numerics/min.h:
 include/./numerics/max.h:
 include/./numerics/interpolation1d.h:
@@ -54,3 +61,7 @@ include/modules/fluids/fluids.h:
 include/modules/fluids/diffusion1d.h:
 include/core/global_config.h:
 include/utils/matrix.h:
+include/./modules/neural_networks/neural_networks.h:
+include/./modules/neural_networks/datasets/spiral.h:
+include/./modules/neural_networks/layers/dense_layer.h:
+include/./modules/neural_networks/activation_functions/ReLU.h:

src/main.cpp

@@ -8,6 +8,8 @@
 #include "./modules/mesh/mesh.h"
 #include "modules/fluids/fluids.h"
 
+#include "./modules/neural_networks/neural_networks.h"
+
 
 
 //#include <iostream>
@@ -20,6 +22,39 @@
 
 int main(int argc, char const *argv[])
 {   
+
+    utils::Mf X(10,2, 0);
+    utils::Vf y(10, 0);
+    neural_networks::create_spital_data(100, 3, X, y);
+
+    neural_networks::dense_layer<float> layer(2, 3);
+
+    neural_networks::activation_ReLU<float> RelU;
+
+    layer.forward(X);
+    RelU.forward(layer.outputs);
+
+    for (int i = 0; i < 5; ++i){
+        std::cout << RelU.outputs.get_row(i) << std::endl;
+    }
+
+    //layer1_output.print();
+    //layer2_output.print();
+
+    //std::cout << output << std::endl;
+
+
+
+
+
+
+
+
+
+
+
+
+    /*
     utils::Vd a = utils::linspace<double>(1, 10, 10, true);
     a.print();
     mesh::Mesh1D<double> mesh(a);
@@ -40,7 +75,7 @@ int main(int argc, char const *argv[])
 
     fluids::Diffusion1D<double> diffusion(cfg, mesh, Gamma);
     diffusion.assemble(A, b, s);
-
+*/
     
     return 0;
 }