diff --git a/bin/abc_lab b/bin/abc_lab
index 0c1f68d..31c11c2 100755
Binary files a/bin/abc_lab and b/bin/abc_lab differ
diff --git a/include/modules/neural_networks/activation_functions/Softmax.h b/include/modules/neural_networks/activation_functions/Softmax.h
index 837b91c..a6431f8 100644
--- a/include/modules/neural_networks/activation_functions/Softmax.h
+++ b/include/modules/neural_networks/activation_functions/Softmax.h
@@ -4,7 +4,11 @@
 
 #include "./utils/vector.h"
 #include "./utils/matrix.h"
-#include "./utils/random.h"
+
+#include "./numerics/max.h"
+#include "./numerics/matsubtract.h"
+#include "./numerics/exponential.h"
+#include "./numerics/matdiv.h"
 
 
 namespace neural_networks{
@@ -12,15 +16,20 @@ namespace neural_networks{
 	template <typename T>
 	struct activation_softmax{
 
+			utils::Matrix<T> exp_values;
+			utils::Matrix<T> probabilities;
 			utils::Matrix<T> outputs;
 
-	        void forward(utils::Matrix<T> inputs){
-			    //outputs = numerics::max(inputs, T{0});
-			    //outputs.print();
+	        void forward(const utils::Matrix<T> inputs){
+	
+
+			    exp_values = numerics::exponential(numerics::matsubtract(inputs, numerics::max(inputs, "rows"), "col"));
+			    probabilities = numerics::matdiv(exp_values, numerics::matsum(exp_values, "col"), "col");
+
+			    outputs = probabilities;
+
 	        }
 
-
-
 	};
 
 
diff --git a/include/modules/neural_networks/datasets/spiral.h b/include/modules/neural_networks/datasets/spiral.h
index d0dd473..4b95a8b 100644
--- a/include/modules/neural_networks/datasets/spiral.h
+++ b/include/modules/neural_networks/datasets/spiral.h
@@ -10,11 +10,11 @@
 
 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) {
+		template <typename TX, typename Ty>
+		void create_spital_data(const uint64_t samples, const uint64_t classes, utils::Matrix<TX>& X, utils::Vector<Ty>& y) {
 
 			const uint64_t rows = samples*classes;
-			T r, t;
+			TX r, t;
 			uint64_t row_idx;
 
 
@@ -27,34 +27,15 @@ namespace neural_networks{
 
 			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);
+					r = static_cast<TX>(j)/static_cast<TX>(samples);
+					t =  static_cast<TX>(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;
+					X(row_idx, 0) = r*std::cos(t*2.5) + utils::random(TX{-0.15}, TX{0.15});
+					X(row_idx, 1) = r*std::sin(t*2.5) + utils::random(TX{-0.15}, TX{0.15});
+					y[row_idx] = static_cast<Ty>(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];
-		        }
-		    }*/
 		}
 
 
diff --git a/include/modules/neural_networks/layers/dense_layer.h b/include/modules/neural_networks/layers/dense_layer.h
index 28aa973..51a8ba7 100644
--- a/include/modules/neural_networks/layers/dense_layer.h
+++ b/include/modules/neural_networks/layers/dense_layer.h
@@ -23,7 +23,7 @@ namespace neural_networks{
 	        	
 	        	weights.random(n_inputs, n_neurons, -1, 1);
 	        	biases.resize(n_neurons, T{0});
-	        	weights.print();
+	        	//weights.print();
 	        	//outputs.resize()
 	        	
 	        }
diff --git a/include/modules/neural_networks/loss/Loss _CategoricalCrossentrophy.h b/include/modules/neural_networks/loss/Loss _CategoricalCrossentrophy.h
new file mode 100644
index 0000000..fa4c6c2
--- /dev/null
+++ b/include/modules/neural_networks/loss/Loss _CategoricalCrossentrophy.h	
@@ -0,0 +1,34 @@
+#pragma once
+
+#include "./core/omp_config.h"
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+
+
+namespace neural_networks{
+
+	template <typename Td, typename Ti>
+	struct Loss{
+
+			utils::Matrix<Td> sample_losses;
+			Td data_losses;
+
+			virtual utils::Vector<Td> forward(const utils::Matrix<Td>& output, const utils::Matrix<Ti>& y) = 0;
+
+			Td calculate(const utils::Matrix<Td>& output, const utils::Matrix<Ti>& y){
+				// Calculate sample losses
+				sample_losses = forward(output, y);
+
+				// Calculate mean loss
+				data_losses = numerics::mean(sample_losses);
+				return data_losses;
+
+			}
+
+
+	};
+
+
+
+} // end namespace neural_networks
\ No newline at end of file
diff --git a/include/modules/neural_networks/loss/loss.h b/include/modules/neural_networks/loss/loss.h
new file mode 100644
index 0000000..fa4c6c2
--- /dev/null
+++ b/include/modules/neural_networks/loss/loss.h
@@ -0,0 +1,34 @@
+#pragma once
+
+#include "./core/omp_config.h"
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+
+
+namespace neural_networks{
+
+	template <typename Td, typename Ti>
+	struct Loss{
+
+			utils::Matrix<Td> sample_losses;
+			Td data_losses;
+
+			virtual utils::Vector<Td> forward(const utils::Matrix<Td>& output, const utils::Matrix<Ti>& y) = 0;
+
+			Td calculate(const utils::Matrix<Td>& output, const utils::Matrix<Ti>& y){
+				// Calculate sample losses
+				sample_losses = forward(output, y);
+
+				// Calculate mean loss
+				data_losses = numerics::mean(sample_losses);
+				return data_losses;
+
+			}
+
+
+	};
+
+
+
+} // end namespace neural_networks
\ No newline at end of file
diff --git a/include/modules/neural_networks/neural_networks.h b/include/modules/neural_networks/neural_networks.h
index f905c20..c0d32ab 100644
--- a/include/modules/neural_networks/neural_networks.h
+++ b/include/modules/neural_networks/neural_networks.h
@@ -6,4 +6,7 @@
 #include "layers/dense_layer.h"
 
 
-#include "activation_functions/ReLU.h"
\ No newline at end of file
+#include "activation_functions/ReLU.h"
+#include "activation_functions/Softmax.h"
+
+#include "loss/loss.h"
\ No newline at end of file
diff --git a/include/numerics/exponential.h b/include/numerics/exponential.h
new file mode 100644
index 0000000..cfe59a7
--- /dev/null
+++ b/include/numerics/exponential.h
@@ -0,0 +1,39 @@
+#pragma once
+
+#include <cmath>
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	T exponential(const T a){
+		return std::exp(a);
+	}
+
+	template <typename T>
+	utils::Vector<T>  exponential(const utils::Vector<T>& a){
+		utils::Vector<T> b = a;
+		for (uint64_t i = 0; i < a.size(); ++i){
+			b[i] =  numerics::exponential(a[i]);
+		}
+		return b;
+	}
+
+	template <typename T>
+	utils::Matrix<T>  exponential(const utils::Matrix<T>& A){
+		utils::Matrix<T> B = A;
+		for (uint64_t i = 0; i < A.rows(); ++i){
+			for (uint64_t j = 0; j < A.cols(); ++j){
+				B(i,j) = numerics::exponential(A(i,j));
+			}
+		}
+		return B;
+	}
+
+
+
+} // namespace numerics
+
diff --git a/include/numerics/matadd.h b/include/numerics/matadd.h
index 58853e8..d3890c8 100644
--- a/include/numerics/matadd.h
+++ b/include/numerics/matadd.h
@@ -1,15 +1,12 @@
 #ifndef _matadd_n_
 #define _matadd_n_
 
-
+#include "./utils/vector.h"
 #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) {
 	    
diff --git a/include/numerics/matdiv.h b/include/numerics/matdiv.h
new file mode 100644
index 0000000..c8e1312
--- /dev/null
+++ b/include/numerics/matdiv.h
@@ -0,0 +1,38 @@
+#ifndef _matdiv_n_
+#define _matdiv_n_
+
+
+#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+
+
+namespace numerics{
+
+// ---------------- Serial baseline ----------------
+	template <typename T>
+	utils::Matrix<T> matdiv(const utils::Matrix<T>& A, const utils::Vector<T>& b, std::string method){
+        
+        utils::Matrix<T> C = A;
+
+        if (method == "row"){
+            for (uint64_t i = 0; i < A.rows(); ++i){
+                for (uint64_t j = 0; j < A.cols(); ++j){
+                    C(i,j) /= b[j];
+                }
+            }
+        }else if (method == "col"){
+            for (uint64_t i = 0; i < A.rows(); ++i){
+                for (uint64_t j = 0; j < A.cols(); ++j){
+                    C(i,j) /= b[i];
+                }
+            }
+        }else{
+            throw std::runtime_error("matdiv: choose div by: 'row' or 'col'");
+        }
+        return C;
+	}
+
+
+} // namespace numerics
+
+#endif // _matdiv_n_
\ No newline at end of file
diff --git a/include/numerics/matsubtract.h b/include/numerics/matsubtract.h
new file mode 100644
index 0000000..b58e881
--- /dev/null
+++ b/include/numerics/matsubtract.h
@@ -0,0 +1,102 @@
+#ifndef _matsubtract_n_
+#define _matsubtract_n_
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+
+namespace numerics{
+	
+	template <typename T>
+	void inplace_matsubtract_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_matsubtract_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_matsubtract_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_matsubtract_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> matsubtract_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_matsubtract_colvec(B, x);
+
+	    return B;
+	}
+
+	template <typename T>
+	utils::Matrix<T> matsubtract_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_matsubtract_rowvec(B, x);
+
+	    return B;
+	}
+
+	template <typename T>
+	utils::Matrix<T> matsubtract(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("matsubtract: too many options for dimensions");
+			} else if (rows == N){
+				return matsubtract_rowvec(A, x);
+			} else if (cols == N){
+				return matsubtract_colvec(A, x);
+			}else{
+				throw std::runtime_error("matsubtract: undefined fault - auto");
+			}
+	    }else if(method=="row"){
+				return matsubtract_rowvec(A, x);
+		} else if (method=="col"){
+			return matsubtract_colvec(A, x);
+		}else{
+			throw std::runtime_error("matsubtract: undefined fault - defined method");
+		}
+    }
+
+
+
+
+} // namespace numerics
+
+#endif // _matsubtract_n_
\ No newline at end of file
diff --git a/include/numerics/matsum.h b/include/numerics/matsum.h
new file mode 100644
index 0000000..9bf0ebd
--- /dev/null
+++ b/include/numerics/matsum.h
@@ -0,0 +1,39 @@
+#ifndef _matsum_n_
+#define _matsum_n_
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+
+namespace numerics{
+	
+	template <typename T>
+	utils::Vector<T> matsum(utils::Matrix<T>& A, std::string method) {
+
+		utils::Vector<T> b;
+
+		if (method == "row"){
+			b.resize(A.cols(), T{0});
+			for (uint64_t i = 0; i < A.cols(); ++i){
+				for (uint64_t j = 0; j < A.rows(); ++j){
+					b[i] += A(j, i); 
+				}
+			}
+		}else if (method == "col"){
+			b.resize(A.rows(), T{0});
+
+			for (uint64_t i = 0; i < A.cols(); ++i){
+				for (uint64_t j = 0; j < A.rows(); ++j){
+					b[j] += A(j, i); 
+				}
+			}
+		}else{
+			throw std::runtime_error("matsum: choose sum by: 'row' or 'col'");
+		}
+		return b;
+	}
+
+
+} // namespace numerics
+
+#endif // _matadd_n_
\ No newline at end of file
diff --git a/include/numerics/max.h b/include/numerics/max.h
index 3f8402f..4e0db24 100644
--- a/include/numerics/max.h
+++ b/include/numerics/max.h
@@ -39,13 +39,37 @@ namespace numerics{
 	utils::Matrix<T> max(const utils::Matrix<T>& A, const T b){
 
 		utils::Matrix<T> B = A;
-
 		inplace_max(B, b);
-		
-
 		return B;
 	}
 
+	template <typename T>
+	utils::Vector<T> max(const utils::Matrix<T>& A, std::string method){
+
+		utils::Vector<T> b;
+
+		if (method == "cols"){
+			b.resize(A.cols(), T{0});
+			for (uint64_t i = 0; i < A.cols(); ++i){
+				for (uint64_t j = 0; j < A.rows(); ++j){
+					b[i] = max(A(j, i), b[i]);
+				}
+			}
+		}else if (method == "rows"){
+			b.resize(A.rows(), T{0});
+			for (uint64_t i = 0; i < A.rows(); ++i){
+				for (uint64_t j = 0; j < A.cols(); ++j){
+					//std::cout << i << ":" << j << std::endl;
+					b[i] = max(A(i, j), b[i]);
+				}
+			}
+		}else{
+			throw std::runtime_error("max: choose 'rows or 'cols'");
+		}
+		return b;
+
+	}
+
 
 
 } // namespace numerics
diff --git a/include/numerics/mean.h b/include/numerics/mean.h
new file mode 100644
index 0000000..b8166a4
--- /dev/null
+++ b/include/numerics/mean.h
@@ -0,0 +1,31 @@
+#ifndef _mean_n_
+#define _mean_n_
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+
+namespace numerics{
+	
+	template <typename T>
+	T mean(utils::Vector<T>& A) {
+
+		T mean(T{0});
+	    
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+
+	    for (uint64_t i = 0; i < cols; ++i) {
+	        for (uint64_t j = 0; j < rows; ++j) {
+	            mean += A(j, i);
+	        }
+	    }
+	    mean /= (static_cast<T>(rows)* static_cast<T>(cols));
+	    return mean;
+	}
+
+
+} // namespace numerics
+
+#endif // _mean_n_
\ No newline at end of file
diff --git a/include/numerics/numerics.h b/include/numerics/numerics.h
index fb5d2b7..1c0d2e3 100644
--- a/include/numerics/numerics.h
+++ b/include/numerics/numerics.h
@@ -6,11 +6,16 @@
 #include "./numerics/transpose.h"
 #include "./numerics/inverse.h"
 #include "./numerics/matmul.h"
+#include "./numerics/matdiv.h"
 #include "./numerics/matvec.h"
 #include "./numerics/matadd.h"
+#include "./numerics/matsubtract.h"
+#include "./numerics/matsum.h"
 #include "./numerics/min.h"
 #include "./numerics/max.h"
 #include "./numerics/abs.h"
+#include "./numerics/mean.h"
+#include "./numerics/exponential.h"
 #include "./numerics/interpolation1d.h"                       		// base
 
 
diff --git a/obj/main.d b/obj/main.d
index 2f70e92..c672385 100644
--- a/obj/main.d
+++ b/obj/main.d
@@ -8,9 +8,12 @@ obj/main.o: src/main.cpp include/./core/omp_config.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/matadd.h include/./numerics/min.h \
- include/./numerics/max.h include/./numerics/interpolation1d.h \
+ include/./numerics/matmul.h include/./numerics/matdiv.h \
+ include/./numerics/matvec.h include/./numerics/matadd.h \
+ include/./numerics/matsubtract.h include/./numerics/matsum.h \
+ include/./numerics/min.h include/./numerics/max.h \
+ include/./numerics/mean.h include/./numerics/exponential.h \
+ include/./numerics/interpolation1d.h \
  include/./numerics/interpolation1d/interpolation1d_barycentric.h \
  include/./numerics/interpolation1d/interpolation1d_base.h \
  include/./numerics/interpolation1d/interpolation1d_cubic_spline.h \
@@ -24,7 +27,9 @@ obj/main.o: src/main.cpp include/./core/omp_config.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/./modules/neural_networks/activation_functions/ReLU.h \
+ include/./modules/neural_networks/activation_functions/Softmax.h \
+ include/./modules/neural_networks/loss/loss.h
 include/./core/omp_config.h:
 include/./utils/utils.h:
 include/./utils/vector.h:
@@ -43,10 +48,15 @@ include/./numerics/inverse/inverse_gauss_jordan.h:
 include/./numerics/inverse/inverse_lu.h:
 include/./decomp/lu.h:
 include/./numerics/matmul.h:
+include/./numerics/matdiv.h:
 include/./numerics/matvec.h:
 include/./numerics/matadd.h:
+include/./numerics/matsubtract.h:
+include/./numerics/matsum.h:
 include/./numerics/min.h:
 include/./numerics/max.h:
+include/./numerics/mean.h:
+include/./numerics/exponential.h:
 include/./numerics/interpolation1d.h:
 include/./numerics/interpolation1d/interpolation1d_barycentric.h:
 include/./numerics/interpolation1d/interpolation1d_base.h:
@@ -65,3 +75,5 @@ 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/./modules/neural_networks/activation_functions/Softmax.h:
+include/./modules/neural_networks/loss/loss.h:
diff --git a/obj/main.o b/obj/main.o
index debce12..d2232b6 100644
Binary files a/obj/main.o and b/obj/main.o differ
diff --git a/src/main.cpp b/src/main.cpp
index 83b8280..5b09871 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -24,24 +24,27 @@ int main(int argc, char const *argv[])
 {   
 
     utils::Mf X(10,2, 0);
-    utils::Vf y(10, 0);
+    utils::Vd y(10, 0);
+
+
     neural_networks::create_spital_data(100, 3, X, y);
 
-    neural_networks::dense_layer<float> layer(2, 3);
+    neural_networks::dense_layer<float> dense1(2, 3);
+    neural_networks::activation_ReLU<float> activation1;
+    neural_networks::dense_layer<float> dense2(3, 3);
+    neural_networks::activation_softmax<float> activation2;
 
-    neural_networks::activation_ReLU<float> RelU;
+    dense1.forward(X);
+    activation1.forward(dense1.outputs);
+    dense2.forward(activation1.outputs);
+    activation2.forward(dense2.outputs);
 
-    layer.forward(X);
-    RelU.forward(layer.outputs);
+ 
 
     for (int i = 0; i < 5; ++i){
-        std::cout << RelU.outputs.get_row(i) << std::endl;
+        std::cout << activation2.outputs.get_row(i) << std::endl;
     }
 
-    //layer1_output.print();
-    //layer2_output.print();
-
-    //std::cout << output << std::endl;